超聲對(duì)兩相液體中分散水滴作用的動(dòng)力學(xué)

王公正，李曉蒙，莫潤(rùn)陽，葉 明，楊厚祿

(陜西師范大學(xué)物理學(xué)與信息技術(shù)學(xué)院/陜西省超聲學(xué)重點(diǎn)實(shí)驗(yàn)室，陜西西安 710062)

超聲空化機(jī)理已經(jīng)得到了人們的廣泛關(guān) 注［1-7］，1917年，Rayleigh首先提出無限大不可壓縮液體中單個(gè)球形空化泡運(yùn)動(dòng)學(xué)模擬［5］，由于沒考慮表面張力、含氣量和壓縮性等因素對(duì)空泡運(yùn)動(dòng)的影響，其計(jì)算的結(jié)果有時(shí)與實(shí)際情況不符。Prospertti A通過引入聲輻射這一因素，對(duì)其修正得到R-P方程［8］，但這些都是研究液體中微小氣泡的運(yùn)動(dòng)情況?，F(xiàn)實(shí)生活中，超聲作用下的液-液兩相非均相系統(tǒng)或膜隔離的兩相系統(tǒng)經(jīng)常遇到如下問題:超聲催化乳化體系［9-13］，超聲對(duì)生命體中細(xì)胞膜內(nèi)、外物質(zhì)和能量的傳輸［14-16］，超聲作用下軟物質(zhì)液滴(如:凝膠、囊泡等)會(huì)怎樣變化，以及超聲吸脂減肥過程中［17-19］，超聲對(duì)脂肪粒和脂肪細(xì)胞這些類似于水滴物質(zhì)的作用機(jī)理如何等，這些問題都是令人關(guān)注的。因此關(guān)于超聲場(chǎng)中反相乳液的研究就顯得尤為重要。本文參考空化泡運(yùn)動(dòng)方程的推導(dǎo)方法［20-22］，推導(dǎo)出水滴在超聲作用下的運(yùn)動(dòng)方程，用數(shù)值模擬的方法研究了在超聲作用下油中水滴半徑的變化，為進(jìn)一步研究超聲對(duì)反相乳液作用機(jī)理提供理論依據(jù)。

1 水滴運(yùn)動(dòng)學(xué)方程的理論推導(dǎo)

1.1 模型的簡(jiǎn)化

本文以油為分散相，研究一個(gè)水滴在邊緣無限大的分散相中的運(yùn)動(dòng)，如圖1所示。液體的可壓縮性比較小，但外加高強(qiáng)度超聲作用時(shí)，水滴因外加壓力作用而引起的微小體積變化與半徑為微米級(jí)的初始水滴體積相比，其變化不容忽略，即水滴可壓縮性不能被忽略，而相同外加壓力作用在分散相液體上所引起的微小體積變化與分散相這個(gè)大環(huán)境整個(gè)體積相比是微乎其微的，因而此過程分散相液體可以被假設(shè)為不可壓縮的。假設(shè)水滴在超聲場(chǎng)作用下發(fā)生體積膨脹收縮現(xiàn)象，則肯定會(huì)伴有能量轉(zhuǎn)換存在。由能量守恒推導(dǎo)出多參數(shù)下的水滴運(yùn)動(dòng)方程，并對(duì)模型進(jìn)行以下簡(jiǎn)化:①水滴在運(yùn)動(dòng)過程中始終保持球形;②水滴壁只做徑向運(yùn)動(dòng)，避免因水滴壁運(yùn)動(dòng)的不對(duì)稱性所帶來的復(fù)雜情況;③ 假設(shè)分散相液體為不可壓縮液體，水滴為可壓縮液體;④ 忽略重力影響;⑤ 水滴在運(yùn)動(dòng)過程中無能量損失。

1.2 水滴壁受力分析

由于水滴壁界面兩邊分子性質(zhì)不同，界面上分子所處力場(chǎng)不是各向同性的，因此存在收縮界面的界面張力［23］。如果知道兩相物質(zhì)的表面張力，那么就可以利用調(diào)和平均方程求出兩相之間的界面張力，調(diào)和平均方程主要適用于低能物質(zhì)(如聚合物、有機(jī)液體、水等)之間的界面張力的求解［24］。凸液面表面的附加壓強(qiáng)為［25］

圖1 油相中水滴簡(jiǎn)化模型Fig.1 Simplified model of a water-drop in oil phase

因?yàn)閮上嘁后w分別為油相與水，所以權(quán)衡后為了簡(jiǎn)化，下文界面張力系數(shù)統(tǒng)一取σ12=0.050 0N/m。

無超聲作用時(shí)，油和水滴組成的系統(tǒng)為靜止流體系統(tǒng)，所以其組成界面為等壓面，水滴在油中保持平衡，水滴壁內(nèi)外壓力相等，因此有

其中P0為水滴周圍油的靜壓，R0為水滴的初始半徑。

1.3 水滴動(dòng)力學(xué)方程推導(dǎo)

假設(shè)整個(gè)過程為等溫過程，那么水滴體積僅是壓力的函數(shù)。由等溫壓縮系數(shù)和體積彈性模量進(jìn)行推導(dǎo)，可以得出液體的等溫狀態(tài)方程［26-27］

其中對(duì)于恒定溫度下同種液體，A、C、n是定值，本文取水滴溫度為40℃，則A=3.688 1，C=3 653，n=0.158 11。

在無超聲作用時(shí)，水滴的等溫狀態(tài)方程為

在外加超聲時(shí)，設(shè)水滴壁內(nèi)對(duì)應(yīng)的壓力為P'in，則對(duì)應(yīng)水滴的等溫狀態(tài)方程為

聯(lián)立(1)，(2)式可得:

進(jìn)一步整理可得:

此時(shí)對(duì)應(yīng)的水滴壁外部受力為

其中PA=Pasin(ωt)為超聲波作用在水滴壁上的壓力，Pa為超聲聲壓幅值。

在外力作用下，水滴半徑從R0變化到R，合外力對(duì)水滴做功為

水滴收縮液體獲得的動(dòng)能為

根據(jù)能量守恒，水滴在合外力作用下半徑從R0變化到R，合外力對(duì)水滴做的功就等于水滴所獲得的動(dòng)能，即

此處 P=P'in-P'out。

兩邊對(duì)R微分，整理可得

若進(jìn)一步考慮水滴運(yùn)動(dòng)過程中的能量粘滯損耗，以及水滴振動(dòng)時(shí)向分散相液體輻射聲波而存在的輻射阻尼［22］，代入 P'in，P'out則公式可修正為

在不同參數(shù)影響下的水滴壁運(yùn)動(dòng)方程屬于2階非線性常微分方程，得不出解析解，需要先變形為形如R·=R(t，R)的一階微分方程組，再采用數(shù)值迭代法求出其數(shù)值解。本文采用4-5階的Runge-Kutta算法［28］，利用 MATLAB 軟件進(jìn)行計(jì)算機(jī)模擬，得出模型方程的解。

2 結(jié)果與討論

2.1 聲速及分散相液體參數(shù)對(duì)水滴半徑變化的影響

2.1.1 聲速對(duì)水滴半徑變化的影響 分散相液體不同即油的種類不同，對(duì)應(yīng)的超聲傳播速度就不同。圖2a是聲速分別取1 000 m/s，1 100 m/s，1 200 m/s，1 300 m/s，1 400 m/s的分散相液體中水滴半徑在5T(5個(gè)周期)內(nèi)的變化，可以看出不同的聲速下水滴半徑變化幾乎完全相同，增長(zhǎng)幅值大約都為初始半徑的25倍。圖2(b)為圖2(a)的二維圖，從圖中可以得出，水滴大約每0.7T完成一次增長(zhǎng)和收縮。收縮達(dá)到的最小值約為初始半徑的0.6倍，增長(zhǎng)和收縮后大約0.3T的時(shí)長(zhǎng)，水滴半徑在初始半徑附近作微小的震蕩，如圖2(c)所示。

2.1.2 密度對(duì)水滴半徑變化的影響 油的種類不同，其對(duì)應(yīng)密度就不同。圖3(a)是分散相液體密 度 分 別 為 700kg/m3，800kg/m3，900kg/m3，950kg/m3時(shí)水滴半徑變化情況，圖3(b)為圖3(a)的二維圖。從圖中可以看出，分散相液體密度越大，水滴增長(zhǎng)振幅越小。密度從700kg/m3增加到950kg/m3的過程中，水滴增長(zhǎng)幅值大約從水滴初始半徑的25.6倍降到22.1倍，水滴壓縮幅值大約為初始半徑的0.6倍。分散相液體密度不同，但水滴增長(zhǎng)壓縮所用的時(shí)長(zhǎng)及作微小震蕩所用時(shí)長(zhǎng)幾乎一樣，分別為0.7T和0.3T。

2.1.3 粘滯系數(shù)對(duì)水滴半徑變化的影響 油的種類不同，粘滯系數(shù)就不同。圖4(a)為分散相液體粘滯系數(shù)分別取1×10-3Pa·s，50×10-3Pa·s，100 ×10-3Pa·s，150 ×10-3Pa·s時(shí)水滴半徑變化，圖4(b)為圖4(a)的二維圖。從圖中可看出，分散相液體粘滯系數(shù)越小，水滴振動(dòng)的振幅越大，粘滯系數(shù)從1×10-3Pa·s增加到150×10-3Pa·s的過程中，水滴增長(zhǎng)幅值從初始半徑的大約24倍減小到5倍左右，水滴收縮幅值大約為初始半徑的0.6倍。同時(shí)隨著粘滯系數(shù)的增加，水滴增長(zhǎng)收縮所用時(shí)長(zhǎng)逐漸減小，而作微小震蕩所用時(shí)長(zhǎng)逐漸增長(zhǎng)，例如當(dāng)粘滯系數(shù)為1×10-3Pa·s時(shí)，水滴增長(zhǎng)收縮所用時(shí)長(zhǎng)為0.7T，作微小震蕩所用時(shí)長(zhǎng)為0.3T，當(dāng)粘滯系數(shù)增長(zhǎng)到150×10-3Pa·s時(shí)，水滴增長(zhǎng)收縮所用時(shí)長(zhǎng)減小到大約0.45T，作微小振動(dòng)所用時(shí)長(zhǎng)大約增長(zhǎng)到0.55T。

圖2 不同聲速下水滴半徑的變化Fig.2 Variation of water-drop radius at the different sound velocity

圖3 不同密度下水滴半徑的變化Fig.3 Variation of water-drop radius at the different density of oil phase

2.1.4 靜壓對(duì)水滴半徑變化的影響 圖5(a)為分散相液體靜壓分別取2×1.013×105Pa，3×1.013×105Pa，4×1.013×105Pa時(shí)水滴半徑的變化情況，圖5(b)為圖5(a)的二維圖。從中可以看出，當(dāng)靜壓不小于超聲激勵(lì)聲壓幅值(2×1.013×105Pa)時(shí)，水滴半徑在初始半徑附近作小幅度的震蕩，并且隨著靜壓的增大，水滴增長(zhǎng)收縮幅值都越來越小，當(dāng)靜壓和超聲激勵(lì)聲壓幅值相等時(shí)，水滴增長(zhǎng)幅值大約只有初始半徑的1.14倍，收縮幅值大約只能達(dá)到初始半徑的0.97倍。而從圖2，3，4可以得知，當(dāng)超聲激勵(lì)聲壓幅值與靜壓分別取2×1.013×105Pa與1.013×105Pa時(shí)，即超聲激勵(lì)聲壓幅值大于靜壓時(shí)，水滴都會(huì)有相對(duì)明顯的增長(zhǎng)和收縮，水滴增長(zhǎng)幅值可達(dá)到初始半徑的25倍左右，收縮幅值幾乎為初始半徑的0.6倍?？梢娕c超聲激勵(lì)聲壓幅值大于靜壓相比，靜壓不小于超聲激勵(lì)聲壓幅值時(shí)，水滴更不容易增長(zhǎng)和壓縮，在激勵(lì)聲壓幅值不變的情況下，靜壓越大水滴體積越不容易發(fā)生變化。

圖4 不同粘滯系數(shù)下水滴半徑的變化Fig.4 Variation of water-drop radius at the different of the viscous coefficient

圖5 不同靜壓下水滴半徑的變化Fig.5 Variation of water-drop radius at the different of the static pressure

2.2 初始半徑及外界超聲激勵(lì)參量對(duì)水滴半徑變化的影響

2.2.1 初始半徑對(duì)水滴半徑變化的影響 圖6為水滴初始半徑分別取 3μm，6μm，9μm 時(shí)水滴半徑在超聲激勵(lì)下的變化。從圖中可以看出，水滴初始半徑越大，水滴半徑變化幅度越小，在初始半徑為3μm時(shí)，水滴的增長(zhǎng)幅度可達(dá)初始半徑的45倍左右，當(dāng)初始半徑為9μm時(shí)，水滴的增長(zhǎng)幅度減小為初始半徑的15倍左右。假設(shè)水滴初始半徑一直增大，那么對(duì)應(yīng)的水滴增長(zhǎng)幅度就會(huì)一直減小，最終減小到不再增長(zhǎng)，也就是當(dāng)水滴初始半徑足夠大時(shí)，即使在超聲激勵(lì)下水滴也很難發(fā)生體積變化。

2.2.2 超聲激勵(lì)聲壓幅值對(duì)水滴半徑變化的影響 圖7(a)為超聲激勵(lì)聲壓幅值分別取2×1.013×105Pa，4 ×1.013 ×105Pa，8 ×1.013×105Pa時(shí)水滴半徑的變化情況。從圖中可以看出，隨著超聲激勵(lì)聲壓幅值的小幅度增大，水滴增長(zhǎng)振幅逐漸變大，水滴在一個(gè)周期內(nèi)增長(zhǎng)收縮所用時(shí)長(zhǎng)增加，在初始半徑附近做微小振動(dòng)時(shí)長(zhǎng)減短。當(dāng)超聲激勵(lì)聲壓幅值為2×1.013×105Pa時(shí)，水滴增長(zhǎng)幅值為初始半徑的24倍左右，增長(zhǎng)收縮所用時(shí)長(zhǎng)大約為0.7T，在初始半徑附近做微小振動(dòng)所用時(shí)長(zhǎng)大約為0.3T，當(dāng)超聲激勵(lì)聲壓幅值為8×1.013×105Pa時(shí)，水滴增長(zhǎng)幅值增長(zhǎng)到初始半徑的92倍左右，增長(zhǎng)收縮所用時(shí)長(zhǎng)大約為1T，水滴半徑幾乎不再在初始半徑附近做微小振動(dòng)。

圖6 不同初始半徑下水滴半徑的變化Fig.6 Variation of water-drop radius at the different of the initial radius

圖7(b)為超聲激勵(lì)聲壓幅值分別為10×1.013×105Pa，100×1.013×105Pa，200×1.013×105Pa，500×1.013×105Pa時(shí)水滴半徑變化情況，圖7(c)為圖7(b)的二維圖。從圖中可以看出，隨著超聲激勵(lì)聲壓幅值大幅度增長(zhǎng)，水滴增長(zhǎng)振幅在逐漸變大的同時(shí)，各個(gè)周期水滴半徑的增長(zhǎng)幅值并不相同，但2T，4T時(shí)間段水滴增長(zhǎng)振幅會(huì)周期性地逐漸趕上1T，3T，5T時(shí)間段的增長(zhǎng)振幅?？梢娂?lì)聲壓越大，水滴半徑增長(zhǎng)越大，當(dāng)激勵(lì)聲壓足夠大時(shí)，奇偶周期段水滴增長(zhǎng)幅度不同，但偶數(shù)周期段水滴半徑的變化會(huì)隨著激勵(lì)聲壓幅值的大幅度增大而周期性地趕上奇數(shù)周期段水滴半徑的變化幅度。

圖7 不同激勵(lì)聲壓幅值下水滴半徑的變化Fig.7 Variation of water-drop radius at the different of the pressure magnitude

2.2.3 超聲激勵(lì)頻率對(duì)水滴半徑變化的影響圖8(a)為超聲頻率分別為2×104Hz，4×104Hz，1×105Hz，2×105Hz時(shí)水滴半徑的變化。從中可以看出，隨著超聲頻率的增大，水滴增長(zhǎng)幅值逐漸變小，當(dāng)頻率為2×104Hz時(shí)，水滴增長(zhǎng)幅值大約為初始半徑的106倍，當(dāng)頻率為2×105Hz時(shí)，水滴增長(zhǎng)幅值減小到初始半徑的11倍左右。

圖8(b)為超聲頻率，分別為1×106Hz，2×106Hz，1×107Hz時(shí)水滴半徑的變化。從中可以看出，隨著超聲頻率的增大水滴振幅逐漸減小的同時(shí)，增長(zhǎng)出現(xiàn)次峰值的個(gè)數(shù)卻在增多。頻率為1×106Hz時(shí)，水滴每做一次劇烈收縮，半徑增長(zhǎng)出現(xiàn)2個(gè)峰值，頻率為2×106Hz時(shí)，水滴做一次劇烈收縮，半徑增長(zhǎng)有時(shí)會(huì)出現(xiàn)3，4個(gè)峰值，當(dāng)頻率為1×107Hz時(shí)，水滴半徑幾乎不再發(fā)生變化。

圖8 不同激勵(lì)頻率下水滴半徑的變化Fig.8 Variation of water-drop radius at the different ultrasonic frequency

3 結(jié)語

兩相液體中分散水滴在超聲場(chǎng)作用下所發(fā)生的變化，不僅僅取決于水滴本身和分散相液體的參數(shù)，即“內(nèi)因”，也取決于外界所加超聲激勵(lì)參數(shù)即“外因”。對(duì)于不同的分散相液體所對(duì)應(yīng)的不同聲速、密度、粘滯系數(shù)而言，在激勵(lì)頻率為2×104Hz、激勵(lì)聲壓幅值為2×1.013×105Pa的超聲作用下，聲速對(duì)水滴半徑變化幾乎沒影響，密度和粘滯系數(shù)越小，水滴半徑變化越明顯，靜壓越大，水滴越不容易發(fā)生體積變化，初始半徑越小，水滴變化幅度越大;在分散相液體確定的情況下，超聲激勵(lì)聲壓幅值越大，水滴半徑增長(zhǎng)幅度越大，超聲激勵(lì)頻率越高，水滴半徑增長(zhǎng)越小，當(dāng)激勵(lì)頻率為1×107Hz時(shí)，水滴半徑幾乎處于初始半徑的穩(wěn)定狀態(tài)。在超聲激勵(lì)下，與水滴壓縮相比，水滴增長(zhǎng)更明顯些。通過分析油相中水滴在超聲場(chǎng)中半徑變化的模擬計(jì)算結(jié)果，可以推測(cè)出反相乳液系統(tǒng)或類似于反相乳液系統(tǒng)，其內(nèi)部近似水滴部分在超聲作用下將發(fā)生增長(zhǎng)和微小收縮現(xiàn)象，并且增長(zhǎng)收縮程度不只取決于外加超聲場(chǎng)的各個(gè)參量大小，而且取決于所研究系統(tǒng)本身的一些參量。這些結(jié)果會(huì)對(duì)超聲在生物醫(yī)學(xué)工程研究方面提供指導(dǎo)。

本文在水滴動(dòng)力學(xué)方程推導(dǎo)的過程中有的因素并沒有考慮進(jìn)去，如分散相液體相對(duì)更小的可壓縮性等，這有待于進(jìn)一步深化研究。

［1］ 吳曉霞，張華余，馬空軍.超聲空化泡運(yùn)動(dòng)特性的研究進(jìn)展［J］.應(yīng)用聲學(xué)，2012，31(6):416-422.

［2］ CUNHA F R，ALBERNAZ D L.Oscillatory motion of a spherical bubble in a non-Newtonian fluid［J］.Journal of Non-Newtonian Fluid Mechanics，2013，191:35-44.

［3］ 張鵬利，林書玉.聲場(chǎng)作用下兩空化泡相互作用的研究［J］.物理學(xué)報(bào)，2009，58(11):7797-7801.

［4］ 沈壯志，林書玉.聲場(chǎng)中水力空化泡的動(dòng)力學(xué)特性［J］.物理學(xué)報(bào)，2011，60(8):385-394.

［5］ HEGEDUS F，HOS C，KULLMANN L.Stable period 1，2 and 3 structures of the harmonically excited Rayleigh-Plesset equation applying low ambient pressure［J］.Journal of Applied Mathematics，2013，78(6):1179-1195.

［6］ BADER K B，MOBLEY J，CHURCH C C，et al.The effect of static pressure on the strength of inertial cavitation events［J］.Journal of The Acoustical Society of America，2012，132(4):2286-2291.

［7］ DABIRI S，SIRIGNANO W A，JOSEPH D D.Interaction between a cavitation bubble and shear flow［J］.Journal of Fluid Mechanics，2010，651:93-116.

［8］ PROSPERTTI A，LEZZI A.Bubble dynamics in a compressible liquid，part 1:first-order theory［J］.Journal of Fluid Mechanics，1986，168:457-478.

［9］ LEONG T S H，WOOSTER T J，KENTISH S E，et al.Minimising oil droplet size using ultrasonic emulsification［J］.Ultrasonics Sonochemistry，2009，16(6):721-727.

［10］ CUCHEVAL A，CHOW R C Y.A study on the emulsification of oil by power ultrasound［J］.Ultrasonics Sonochemistryy，2008，15(5):916-920.

［11］GAIKWAD S G，PANDIT A B.Ultrasound emulsification:Effect of ultrasonic and physicochemical properties on dispersed phase volume and droplet size［J］.Ultrasonics Sonochemistry，2008，15(4):554-563.

［12］ MAA Y F，HSU C C.Performance of sonication and micro fluidization for liquid-liquid emulsification［J］.Pharmaceutical Development and Technology，1999，4(2):233-240.

［13］ KALSA O，MICHON C，YANNIOTIS S，et al.Ultrasonic energy input influence on the production of submicron o/w emulsions containing whey protein and common stabilizers［J］.Ultrasonics Sonochemidtry，2013，20(3):881-891.

［14］張德俊.超聲空化及其生物醫(yī)學(xué)效應(yīng)［J］.中國超聲醫(yī)學(xué)雜志，1995，11(7):510-512.

［15］ ZHOU Yun，CUI Jian-min，CHENRI X D.Dynamics of Sonoporation Correlated with Acoustic Cavitation Activities［J］.2008 Biophysical Journal，2008，94(7):L51-L53.

［16］ PARK Ju-young，F(xiàn)AN Zhen-zhen，CHERI X D.Effects of shear stress cultivation on cell membrane disruption and intracellular calcium concentration in sonoporation of endothelial cells［J］.Journal of Biomechanics.2011，44(1):164-169.

［17］ PALUMBO P，CINQUE B，MICONI G，et al.Biological effects of low frequency high intensity ultrasound application on ex vivo human adipose tissue［J］.International Journal of Immunopathology and Pharmacology，2011，24(2):411-422.

［18］ PINE J L，SMITH L J，H.MELINDA J，et al.Ultrasound-assisted lipoplasty［J］.Plastic Surgical Nursing:Official Journal of the American Society of Plastic and Reconstructive Surgical Nurses，2003，23(3):101-108.

［19］ NAVARRO-VIANA F.Rhytidectomy assisted with ultrasound techniques:the ultra-lipo-lift technique［J］.Aesthetic Plastic Surgery，2001，25(3):175-180.

［20］BEHNIA S，ZAHIR H，YAHUAVI M，et al.Observations on the dynamics of bubble cluster in an ultrasonic field［J］.Nonlinear Dynamics，2013，72(3):561-574.

［21］HEGEDUS F，KOCH S，GAREN W，et al.The effect of high viscosity on compressible and incompressible Rayleigh-Plesset-type bubble models［J］.International Journal of Heat and Fluid Flow，2013，42:200-208.

［22］ PROSPERETTI A，LEZZI A.Bubble dynamics in a compressible liquid，part 1:first-order theory［J］.Journal of Fluid Mechanics，1986，168:457-478.

［23］ WANG P，ANDERKO A.Modeling interfacial tension in liquid-liquid systems containing electrolytes［J］.Industrial＆ Engineering Chemistry Research，2013，52(20):6822-6840.

［24］杜永娟，陳國榮，胡福增.材料表界面［M］.上海:華東理工大學(xué)出版社，2001.

［25］ PETERS F，ARABALI D.Interfacial tension between oil and water measured with a modified contour method［J］.Colloid and Surfaces A-Physicochemical and Engineering Aspects，2013，426:1-5.

［26］李志旺，嚴(yán)家騄.一個(gè)新的液體狀態(tài)方程［J］.工程熱物理學(xué)報(bào)，1990，11(2):126-128.

［27］王克強(qiáng).一個(gè)新的液體等溫狀態(tài)方程［J］.荊州師專學(xué)報(bào)(自然科學(xué)版)，1992，15(5):64-67.

［28］史振方.聲空化優(yōu)化的數(shù)值模擬研究［D］.西安:陜西師范大學(xué)，2009.