劉 虎，強洪夫，韓亞偉，陳福振

（第二炮兵工程大學(xué)，陜西 西安710025）

0 引言

冪律型流體是指流變行為符合冪定律的流體，這類流體在石油工業(yè)、航空航天等領(lǐng)域應(yīng)用較為廣泛[1]。突縮型管路流動是指流體在管路中流動時，由于管路截面面積突然縮小而產(chǎn)生的流動現(xiàn)象[2]。對于火箭發(fā)動機而言，研究冪律型流體的突縮型管路流動特性，獲取速度、剪切速率、粘度等的分布及變化情況，一方面可以為具有冪律型本構(gòu)關(guān)系的推進劑的管路設(shè)計提供參考；另一方面，突縮型管路也可以看作噴注器噴口的簡化形式。因此，對噴注器設(shè)計也具有一定指導(dǎo)意義。

光滑粒子流體動力學(xué)（Smoothed Particle Hydrodynamics,SPH） 方法是求解流體流動的一種無網(wǎng)格粒子法。由于流體粘性項中包含的速度的二階導(dǎo)數(shù)很難求解，因此，早期SPH方法一般只用于求解Euler方程，即局限于非粘流動。目前，對于牛頓流體，較實用的方法是將有限差分方法與SPH方法相結(jié)合，將二階導(dǎo)數(shù)降階求解[3]；對于非牛頓流體，由于其本構(gòu)關(guān)系種類較多，且關(guān)系式較為復(fù)雜，SPH方法的研究尚處于起步階段[4-6]。

本文主要進行了以下工作：推導(dǎo)了廣義流體動力學(xué)控制方程的SPH離散式，應(yīng)用人工應(yīng)力方法消除拉伸不穩(wěn)定現(xiàn)象，應(yīng)用XSPH方法規(guī)范粒子分布秩序，應(yīng)用罰方法施加邊界條件；提出了冪律型本構(gòu)關(guān)系的SPH求解方法，推導(dǎo)了其計算公式；分別用本文方法及有限體積方法（Finite Volume Method,FVM）進行了Poiseuille流流動過程仿真，通過對比，證明了本文方法的正確性；將本文方法應(yīng)用于冪律型流體的突縮型管路流動仿真研究，獲得了管路流動特性，并與牛頓型流體（水）的流動特性進行了對比，分析了冪律型流體流動特性的成因。

1 SPH基本理論

1.1 控制方程及狀態(tài)方程

本文采用的流體動力學(xué)控制方程為：

對于動量方程（2） 中的壓力項p，采用弱可壓縮狀態(tài)方程進行求解：

式中：p0為參考壓強；γ為常數(shù)，本文中取γ=7；p0和γ共同用于控制計算中流體密度在其常態(tài)密度附近的震蕩幅度。

1.2 SPH基本方法及控制方程的離散

在SPH算法中，連續(xù)的流場離散成為一系列有相互作用的粒子，通過核函數(shù)估計技術(shù)在這些粒子上離散控制方程組，得到一組描述各粒子物理量隨時間變化的常微分方程組，即SPH基本方程組，再對這組方程采用相應(yīng)的常微分方程組求解方法來推進時間進程的求解。

式（1）-（3）的 SPH 離散式為：

式中：i，j為粒子編碼；N為粒子i支持域內(nèi)的粒子數(shù)；ρi為粒子i的密度；mj，ρj為粒子j的質(zhì)量和密度；為核函數(shù)，它的選取直接影響計算的誤差和穩(wěn)定性，本文選用三次樣條核函數(shù)；表示核函數(shù)對粒子i的偏導(dǎo)數(shù)；h為光滑長度，表示核函數(shù)不顯著為零的取值范圍，控制著SPH粒子的影響域。

式（7） 被稱為XSPH方法，其基本思想是通過施加臨近粒子的影響使自身的運動速度與臨近粒子的平均速度相近，使粒子分布更加有序，消除由于分布不均勻帶來的粒子的非物理聚集問題。 ξ(0≤ξ≤1)是一個常數(shù)，本文中ξ=0.3。

1.3 邊界條件的施加

采用罰方法施加邊界條件，邊界力的計算方法為[8]：

1.4 時間積分

采用蛙跳 (Leap-frog)方法對SPH離散方程進行求解，它對時間是二階精度，具有存儲量低，計算效率高的特點。

式中：φ表示密度ρ及速度v； xi為粒子i的位置坐標(biāo)。為了使計算過程穩(wěn)定，采用考慮具有粘性耗散的時間步長表達式：

對于本文研究的冪律型流體而言，η為流體表觀粘度。

2 冪律型本構(gòu)關(guān)系的SPH處理

冪律型流體的本構(gòu)關(guān)系式為：

流體的剪切速率可以表示為：

二維情況下，式（13）可寫作：

對式（14）在粒子i處進行粒子估計：

式（15） 即為粒子i處剪切速率張量的SPH離散求解式。

借鑒牛頓流體中粘性的處理方法，本文應(yīng)用下式來計算冪律型流體的粘性項：

當(dāng)n<1時，

當(dāng)η0取較大值時，式 （18） 可以代替式（17）近似描述冪律型流體的本構(gòu)關(guān)系。

3 測試算例

應(yīng)用Poiseuille算例驗證第3節(jié)冪律型本構(gòu)關(guān)系的SPH求解方法的正確性。Poiseuille流的模型為：流體在分別位于y=0和y=l的兩塊平行且無限大的固定平板間流動。初始靜止的流體由于受到體力F（例如壓力梯度或外力）的作用而在兩平板間逐漸流動，最終達到穩(wěn)定狀態(tài)。

本算例中η0=106Pa·s，以近似表示式 （17）描述的冪律型流體本構(gòu)關(guān)系。為增大時間步長，采用了較大尺寸的計算模型，其中，模型總體尺寸 0.4 m×0.4 m，粒子間距 0.01 m，粒子數(shù) 1 600，其中，流體粒子數(shù)1 520；邊界粒子采用虛粒子，粒子數(shù)80。Poiseuille流模型的初始粒子分布如圖1所示。

圖1 Poiseuille流初始粒子分布圖Fig.1 Initial particle distribution of Poiseuille flow

圖2 Poiseuille流穩(wěn)定時x=0處的速度分布圖Fig.2 Velocity distribution of particles at x=0 when Poiseuille flow is stable

流體密度 ρ=1000 kg/m3。體力 F=2×10-4m/s2。流變指數(shù)n=0.8，稠度系數(shù)k=31.48 Pa·s。應(yīng)用FVM方法對同等規(guī)模的網(wǎng)格模型進行了對比計算。圖2為模擬達到穩(wěn)定狀態(tài)時，分別用SPH方法及FVM方法計算得到的x=0處的粒子速度分布圖。由圖2可以看出，SPH計算結(jié)果與FVM計算結(jié)果在速度變化趨勢上完全一致，二者的最大相對誤差為4.03%。

4 突縮型管路流動仿真

采用的突縮型管路模型如圖3所示。推板以恒定的速度V=0.1 m/s推動管路內(nèi)的流體向外流出，管路直徑D=5.5×10-3m，管路截面面積與突縮口截面面積比Spipe/Sout=2，計算中不考慮重力作用。

圖3 突縮型管路流動示意圖Fig.3 Schematic of flow in sharply contractive pipe

模型粒子數(shù)5 552，其中，流體粒子5 000個，邊界粒子552個；流體粒子間距1.1×10-4m，邊界粒子間距5.5×10-5m。

采用牛頓型流體進行了突縮管路流動特性模擬，與冪律型流體的流動特性進行對比。牛頓型流體為水，冪律型流體的物質(zhì)參數(shù)取自文獻 [3]，其中，稠度系數(shù)k=41.65 Pa·s，流變指數(shù) n=0.224，密度ρ=1000 kg/m3。結(jié)合文獻 [3]中的表觀粘度與剪切速率實驗數(shù)據(jù)，取η0=0.3 Pa·s。

圖4為突縮型管路流動仿真結(jié)果，通過對比可以發(fā)現(xiàn)，對于本算例，冪律型流體的流動特性與水有很大區(qū)別，表現(xiàn)為：

1） 相同時間內(nèi)，冪律型流體的流動距離較小，流體液柱直徑較大；

2）冪律型流體的液柱頭部呈凸面形。

圖4 突縮型管路流動仿真結(jié)果Fig.4 Simulation results of flow in sharply contractive pipe

圖5突縮型管路流動的速度場變化過程Fig.5 Developing process of velocity field of flow in sharply contractive pipe

圖5 為突縮型管路流動速度場的變化過程圖，其中，上圖為冪律型流體的速度變化過程圖；下圖為水的速度變化過程圖。通過對比可以得出：

1） 各時刻的冪律型流體的流動速度均小于水流速度；

2） 流體自突縮口流出前（t＝0.00025 s），冪律型流體液速度呈凹面型，中間速度低，突縮口處速度高；流出后（t＝0.00225 s以后），呈凸面型，中間速度高，兩側(cè)速度低；水流速度始終表現(xiàn)為兩側(cè)速度高，中間速度低。

圖6為冪律型流體液經(jīng)過突縮型管路轉(zhuǎn)角處的剪切速率分布等值線圖?？梢钥闯?，流體的剪切速率值圍繞突縮型管路轉(zhuǎn)角呈環(huán)狀分布，越靠近轉(zhuǎn)角處，剪切速率值越大。

圖6 冪律型流體經(jīng)過突縮型管路轉(zhuǎn)角處的剪切速率等值線圖Fig.6 Isolines of shear rates while power-law fluid is flowing through the sharp edge of sharply contractive pipe

結(jié)合圖4～圖6，可以得到冪律型流體經(jīng)過突縮型管路的流動特性成因：

1）由于冪律型流體的粘性遠比水的粘性大，因此，冪律型流體在管路流動時受到的粘性阻力大，流動速度慢，在流量相同的情況下形成的液柱直徑比水的液柱直徑大；

2） 初始時刻，流體在推板的推動下向前運動，但并未通過突縮口，在突縮口轉(zhuǎn)角處，兩側(cè)流體必須轉(zhuǎn)向通過突縮口，導(dǎo)致速度急劇增大，此外，對于冪律型流體，速度的增大引起剪切速率變大，表觀粘度減小，轉(zhuǎn)角處受到的粘性阻力減小，同樣引起速度的增大，因此，水和冪律型流體在噴注初期，突縮口截面的速度分布均呈凹面形，中間速度低，兩側(cè)速度高；

3） 流體從突縮口流出后，冪律型流體的粘度大，產(chǎn)生的粘性阻力大于牛頓型流體，導(dǎo)致冪律型流體的運動速度迅速減小，由于速度差主要存在于流體兩側(cè)，因此，兩側(cè)的減速作用尤為明顯；此外，流體中間部分的流量比兩側(cè)流量大，也造成中間部分的流速比兩側(cè)流速高，從而導(dǎo)致冪律型流體的液柱頭部形狀及速度分布均呈凸面形。

5 結(jié)論

本文的方法可有效處理冪律型流體的本構(gòu)關(guān)系。

冪律型流體的突縮型管路流動特性與水有很大差別，主要表現(xiàn)為流動阻力大、流速低、流動困難。

在本文研究范圍內(nèi)，粘性較大以及粘性隨剪切速率的變化是形成冪律型流體流動特性的主要原因。

[1]NATAN B,RAHIMI S.The status of gel propellants in year 2000[C]//Proceedings of 5th International Symposium on Special Topics in Chemical Propulsion:Combustion of EnergeticMaterials.BocaRaton:BegelHouse,2000:1-24.

[2]張蒙正,左博．冪律型凝膠推進劑管路中的流動特性[J]．推進技術(shù),2009,30(2):246-250．

[3]BASA M,QUINLAN N J,LASTIWKA M.Robustness and accuracy of SPH formulations for viscous flow[J].International Journal for Numerical Methods in Fluids,2009,60:1127-1148.

[4]FANG Jian-nong,OWENS R G,TACHER L,et al.A numericalstudyoftheSPH method forsimulating transient viscoelastic free surface flows[J].Journal of Non-Newtonian Fluid Mech,2006,139:68-84.

[5]RFIEE A,MANZARI M T,HOSSEINI M.An incom pressible SPH method for simulation of unsteady viscoelastic free-surface flows[J].International Journal of Non-Linear Mechanics,2007,42(10):1210-1223.

[6]FAN X J,TANNER R I,ZHENG R.Smoothed particle hydrodynamics simulation of non-Newtonian moulding flow [J].Journal of Non-Newtonian Fluid Mechanics,2010,165:219-226.

[7]MONAGHAN J J.SPH without a tensile instability[J].Journal of Computational Physics,2000,159:290-311.

[8]王坤鵬.基于接觸邊界條件SPH新方法及其工程應(yīng)用[D].西安:第二炮兵工程大學(xué)201教研室,2008:25-29.