基于SPH方法的冪律流體突縮型管路流動(dòng)過(guò)程仿真

劉 虎，強(qiáng)洪夫，韓亞偉，陳福振

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

0 引言

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

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

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

1 SPH基本理論

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

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

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

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

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

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

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

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

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

1.3 邊界條件的施加

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

1.4 時(shí)間積分

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

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

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

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

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

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

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

對(duì)式（14）在粒子i處進(jìn)行粒子估計(jì)：

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

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

當(dāng)n<1時(shí)，

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

3 測(cè)試算例

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

本算例中η0=106Pa·s，以近似表示式 （17）描述的冪律型流體本構(gòu)關(guān)系。為增大時(shí)間步長(zhǎng)，采用了較大尺寸的計(jì)算模型，其中，模型總體尺寸 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)定時(shí)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方法對(duì)同等規(guī)模的網(wǎng)格模型進(jìn)行了對(duì)比計(jì)算。圖2為模擬達(dá)到穩(wěn)定狀態(tài)時(shí)，分別用SPH方法及FVM方法計(jì)算得到的x=0處的粒子速度分布圖。由圖2可以看出，SPH計(jì)算結(jié)果與FVM計(jì)算結(jié)果在速度變化趨勢(shì)上完全一致，二者的最大相對(duì)誤差為4.03%。

4 突縮型管路流動(dòng)仿真

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

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

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

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

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

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

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

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

圖5突縮型管路流動(dòng)的速度場(chǎng)變化過(guò)程Fig.5 Developing process of velocity field of flow in sharply contractive pipe

圖5 為突縮型管路流動(dòng)速度場(chǎng)的變化過(guò)程圖，其中，上圖為冪律型流體的速度變化過(guò)程圖；下圖為水的速度變化過(guò)程圖。通過(guò)對(duì)比可以得出：

1） 各時(shí)刻的冪律型流體的流動(dòng)速度均小于水流速度；

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

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

圖6 冪律型流體經(jīng)過(guò)突縮型管路轉(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)過(guò)突縮型管路的流動(dòng)特性成因：

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

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

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

5 結(jié)論

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

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

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

[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ìn)劑管路中的流動(dòng)特性[J]．推進(jìn)技術(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.