基于VOF模型的凝膠推進(jìn)劑的偽流動研究

曹欽柳,封 鋒,鄧寒玉

(南京理工大學(xué) 機械工程學(xué)院,江蘇 南京 210094)

基于VOF模型的凝膠推進(jìn)劑的偽流動研究

曹欽柳,封 鋒,鄧寒玉

(南京理工大學(xué) 機械工程學(xué)院,江蘇 南京 210094)

利用VOF模型計算凝膠推進(jìn)劑一次霧化的兩相流動時,為研究引入表面張力后偽流動的特征,開發(fā)了基于同步重疊相加(SOLA)交錯網(wǎng)格算法的非牛頓不可壓流體求解器。以凝膠推進(jìn)劑液滴為對象,研究了VOF模型下的偽流動。結(jié)果表明:使用本方法的驗證算例中數(shù)值解和解析解吻合度較高,誤差在5%以內(nèi);偽流動強度定義中包含表面張力系數(shù),一定程度上抵消了凝膠推進(jìn)劑表面張力計算誤差增大對偽流動強度的影響;4種網(wǎng)格尺度下,偽流動強度的變化趨勢及時間平均值基本一致,說明網(wǎng)格尺度對偽流動無影響;偽流動強度隨著本構(gòu)參數(shù)n的升高而增強,并且當(dāng)n>0.6時上升得更快;偽流動強度與本構(gòu)參數(shù)K呈比例系數(shù)為0.05的線性關(guān)系。

凝膠推進(jìn)劑;表面張力;流體體積函數(shù);偽流動

凝膠推進(jìn)劑作為一種新型的推進(jìn)劑兼具了固體和液體推進(jìn)劑各自的優(yōu)點,具有高安全性、易儲藏、不易泄漏和可重復(fù)啟動的特點,因此得到迅速的發(fā)展。要使凝膠推進(jìn)劑實現(xiàn)高效的燃燒,需要通過相關(guān)措施使其充分霧化。對于液滴的二次霧化問題,已經(jīng)提出相關(guān)的振蕩和破碎模型[1-2]。但是對于凝膠推進(jìn)劑的一次霧化,目前仍無較為完備的理論,相關(guān)研究者提出了一些近似的計算模型,希望能夠從霧化到燃燒進(jìn)行完整的數(shù)值模擬,采用VOF、level set和SPH等方法都取得了較好的效果[3-5]。在液體一次霧化中表面張力是不容忽視的,在數(shù)值模擬中引入表面張力模型后,由于不能精確計算表面張力和界面的壓力跳躍條件,給計算上帶來了數(shù)值誤差。這種數(shù)值誤差的直接后果便是在模擬中出現(xiàn)偽流動現(xiàn)象,當(dāng)偽流動足夠強時,會引起數(shù)值振蕩甚至得到錯誤的結(jié)果。Gueyffier等[6]指出,在一些高密度比和高黏性比的情況下,很多VOF方法都不能按照真實的參數(shù)進(jìn)行計算。

對于偽流動現(xiàn)象國內(nèi)外研究者進(jìn)行了大量研究工作。Lafaurie等[7]基于VOF方法對氣泡的偽流動進(jìn)行了研究,提出了一種衡量偽流動的參數(shù),結(jié)果表明牛頓流體中氣泡的表面張力越強時偽流動就越明顯。Popinet和Zaleski[8]用一種front-tracking方法細(xì)化了網(wǎng)格壓力的梯度離散,較好地處理了表面張力和界面壓力跳躍平衡,使得偽流動強度大大降低。Aleinov等[9]基于核函數(shù)提出了計算表面張力的高階內(nèi)核概念,使得更精確地計算界面曲率成為可能。馬東軍等[10]對VOF方法的偽流動進(jìn)行了研究,得出其強度幾乎為常值,并且不隨網(wǎng)格加密而收斂。不過上述文獻(xiàn)都是針對牛頓流體進(jìn)行的研究,目前針對凝膠等非牛頓流體偽流動現(xiàn)象的研究還未見公開報道。

本文采用VOF模型進(jìn)行界面捕捉,利用CSF模型處理表面張力,開發(fā)了求解不可壓流體控制方程組的求解器,并以凝膠推進(jìn)劑液滴為對象,對偽流動現(xiàn)象進(jìn)行了研究。分析了不同表面張力系數(shù),不同網(wǎng)格尺度和不同冪率模型本構(gòu)參數(shù)對偽流動強度的影響。相關(guān)結(jié)論可以為采用VOF方法研究凝膠推進(jìn)劑的一次霧化提供參考。

1 數(shù)值方法

1.1 計算模型

本文研究凝膠推進(jìn)劑液滴在VOF模型下的偽流動,液滴的直徑、計算區(qū)域大小及其邊界類型如圖1所示。凝膠液滴周圍為空氣,四周為自由邊界,速度零梯度外推。初始時壓力設(shè)為環(huán)境壓力101 325 Pa,速度均為0。

1.2 界面重構(gòu)方法

在兩相共存的網(wǎng)格中,實際上是有界面的。采用施主-受主思想計算時,需要幾何重構(gòu)方法構(gòu)造界面。研究人員提出了一些幾何重構(gòu)方法,如SLIC、FLAIR和PLIC等[11]。其中PLIC方法是由Youngs首先提出的。它通過單個網(wǎng)格的斜線來近似兩相界面,斜線的斜率需要通過周圍網(wǎng)格確定,在重構(gòu)精度上可以達(dá)到二階。因此,本文采用PLIC重構(gòu)方法進(jìn)行數(shù)值計算。如圖2給出了界面法向在第三象限的4種情況。

1.3 黏性模型

凝膠推進(jìn)劑在流動上表現(xiàn)出非牛頓流體的性質(zhì),其流變本構(gòu)方程可寫為如下形式[12]:

(1)

(2)

(3)

Δ可以根據(jù)剪切速度張量各分量求得:

(4)

由上述公式,可以得到應(yīng)力張量各分量表達(dá)式為

(5)

(6)

(7)

式中:η0為表觀黏度,它是衡量非牛頓流體黏性的重要參數(shù),可以表示為

η0=KΔ(n-1)/2

(8)

1.4 表面張力模型

計算表面張力時,本文采用Brackbill提出的連續(xù)表面張力模型(CSF),假設(shè)表面張力為一種體積力,并且連續(xù)分布于兩相界面兩側(cè)的一定厚度內(nèi)。其大小與界面的曲率成正比,則網(wǎng)格內(nèi)兩相界面受的表面張力為[13]

fs=σκ(F)

(9)

式中:σ為表面張力系數(shù),κ為兩相界面的曲率,F為網(wǎng)格中液體相的體積分?jǐn)?shù)。上式求得的表面張力實際是一種壓力對坐標(biāo)的變化量,并作用于網(wǎng)格內(nèi)所有的流體,它作為源項添加到動量方程中。

在牛頓流體的VOF偽流動研究中Popinet[8]等引入了一個衡量偽流動強度的參數(shù),其表達(dá)式為

Cas=|U|μ/σ

(10)

式中:μ為牛頓流體的黏性系數(shù),U為速度矢量。目前還未提出專門針對凝膠偽流動的衡量參數(shù),因此本文嘗試直接用表觀黏度η0代替上述黏性系數(shù)。

表面張力的引入給計算帶來了一定的數(shù)值誤差。采用重構(gòu)方法獲得的兩相界面始終是對真實界面的一種近似,導(dǎo)致表面張力的計算值與實際值存在一定誤差。目前的數(shù)值方法很難處理表面張力帶來的壓力跳躍平衡條件,二者帶來的誤差會驅(qū)動流場的非物理變化,這是產(chǎn)生偽流動的根本原因。需要指出的是,偽流動是計算模型本身的固有特性,提高重構(gòu)精度或者采用更好的數(shù)值方法可以減小偽流動的強度。

2 數(shù)值驗證

本文采用基于交錯網(wǎng)格的SOLA算法[14],壓力和密度值存儲于格心處,速度值則保存在網(wǎng)格線上。壓力項和黏性項采用中心差分格式,對流項采用不完全迎風(fēng)差分格式,即迎風(fēng)格式與中心差分格式的線性組合,本文的線性系數(shù)取為0.5。為驗證本文數(shù)值方法的有效性,選取如圖3并行兩相通道流作為驗證算例,圖中上下為兩平行平板。

3 計算結(jié)果與分析

本文為深入探討VOF模型中凝膠液滴的偽流動現(xiàn)象,對不同表面張力系數(shù),網(wǎng)格尺度和冪率模型本構(gòu)參數(shù)下凝膠液滴的偽流動進(jìn)行了數(shù)值計算?？諝饷芏热?.25kg/m3,凝膠推進(jìn)劑的密度在不同的配比下,其值有所差異,本文直接取為1 000kg/m3。按照冪率模型,當(dāng)剪切速率為0時,其表觀黏度趨向于無窮大,因此數(shù)值計算時需要對上界值進(jìn)行限定,本文取最大表觀黏度值ηmax為10Pa·s。

3.1 表面張力系數(shù)的影響

Lafaurie等[7]基于VOF方法對牛頓流體氣泡的偽流動進(jìn)行了研究,指出表面張力增大時偽流動越明顯。表面張力系數(shù)是表面張力中的一個重要參數(shù)。為研究偽流動強弱隨表面張力系數(shù)的變化規(guī)律,本文選取0.01,0.02,0.03,0.04,0.05,0.06和0.07 7種不同的表面張力系數(shù)值,采用冪率型煤油凝膠(K=13.5,n=0.47),網(wǎng)格數(shù)為2002,對其進(jìn)行數(shù)值計算。

圖5給出了5ms時4種張力系數(shù)下的凝膠液滴表面圖。從圖中可以看出,隨著表面張力的升高,界面扭曲程度增大,說明表面出現(xiàn)了較大的偽流動速度。為定量描述偽流動強度Cas值的變化情況,圖6給出了3個時刻偽流動強度隨表面張力系數(shù)的變化規(guī)律。從圖中可以看出,表面張力系數(shù)增大時,3個時刻的偽流動強度稍有下降或者上升,但是其變化都比較小,基本保持一致。從圖5可以看出,表面張力系數(shù)較大時表面處偽流動帶來的影響明顯強于系數(shù)較小的情況。雖然表面張力系數(shù)大時,造成的界面處偽流動速度大,但同時偽流動強度公式(10)中也包含了表面張力系數(shù)的影響,從而一定程度上抵消了表面張力誤差對偽流動強度的影響。

3.2 網(wǎng)格尺度的影響

數(shù)值上計算兩相界面曲率時存在較大誤差是造成表面張力誤差的原因。如果采用更加精細(xì)的網(wǎng)格,細(xì)化對界面的描述,就會對界面曲率的計算精度產(chǎn)生影響。為研究凝膠液滴的偽流動強度隨網(wǎng)格尺度的變化規(guī)律,選用相同的冪率型煤油凝膠,在不同的網(wǎng)格尺度下,對其進(jìn)行了數(shù)值研究。網(wǎng)格的相關(guān)參數(shù)如表1所示,表中m為網(wǎng)格數(shù)量,md為液滴直徑上網(wǎng)格點的數(shù)量。

表1 網(wǎng)格參數(shù)

圖7為4個算例10ms時的凝膠液滴表面示意圖。從圖中可以看出,網(wǎng)格數(shù)目為502的算例發(fā)生了明顯的扭曲,但是界面比較平緩且褶皺較少。隨著網(wǎng)格數(shù)目的增加,兩相界面都顯著偏離了初始的圓形界面,并且界面的褶皺逐漸增多,界面的銳利性逐漸增強。網(wǎng)格數(shù)目的增加并未有效延緩偽流動的發(fā)生,并且網(wǎng)格數(shù)目增加時界面將在更小的范圍內(nèi)發(fā)生扭曲,反而使得褶皺逐漸增加。因此,網(wǎng)格數(shù)增加可以使得界面在更小的范圍內(nèi)變化,但是這種變化終會逐漸發(fā)展到可觀察的程度。

馬東軍[10]使用VOF模型研究了牛頓流體偽流動現(xiàn)象隨網(wǎng)格尺度的變化,結(jié)果表明Cas值不隨網(wǎng)格尺度的變化而變化,說明偽流動強度不具備網(wǎng)格收斂性。圖8給出了4種網(wǎng)格尺度Cas值在10ms內(nèi)隨時間的變化圖。從圖中可以看出,4種情況Cas值都隨時間逐漸升高,偽流動都隨時間逐漸增強,并且其變化規(guī)律基本一致,其時間平均值分別為0.233,0.221 8,0.222和0.221,時間平均值的標(biāo)準(zhǔn)方差為0.004 2。所以,加密網(wǎng)格后計算量大了,但是對偽流動的時間增長率和強度都未產(chǎn)生多大影響,本文采用的數(shù)值方法在凝膠的偽流動中同樣不具備網(wǎng)格收斂性。因此,通過加密網(wǎng)格來控制偽流動,對提高含表面張力計算的數(shù)值精度是無效的。

3.3 本構(gòu)參數(shù)的影響

在針對牛頓流體的偽流動研究中,一些無量綱參數(shù)是影響偽流動強度的重要因素,如凝膠液滴與空氣的密度比ρl/ρg和黏性比μl/μg,以及衡量表面張力和黏性力相互關(guān)系的Ohnesorge數(shù)[8],但是,對于凝膠推進(jìn)劑,其表觀黏度隨剪切速率而變化,流場中每個網(wǎng)格的表觀黏度都有所差別。本文算例中空氣和凝膠液滴的密度均為常數(shù),密度比為定值,而黏性比和Ohnesorge數(shù)沒有統(tǒng)一的數(shù)值,所以很難用這兩個參數(shù)來研究黏度對偽流動的影響。式(1)的本構(gòu)參數(shù)K和n體現(xiàn)了黏度的變化,因此本文選取不同的K和n值,如表2,網(wǎng)格數(shù)為2002,對凝膠液滴的偽流動進(jìn)行了計算。

表2 凝膠推進(jìn)劑的本構(gòu)參數(shù)

為定量描述偽流動強度的變化,圖9給出了Cas值隨n的變化規(guī)律。從圖中可以看出,隨著n的增大,偽流動的強度逐漸增大,最大值5.18比最小值0.043要高出2個數(shù)量級,并且當(dāng)n>0.6時,偽流動強度的變化對n很敏感。這說明,雖然在目前的算法中偽流動是不可避免的數(shù)值誤差,但是對于不同的流體,其強度是有差別的。非牛頓流體稀化能力越強,相同時間內(nèi)由偽流動帶來的誤差就越小。所以,對于稀化能力較弱的凝膠推進(jìn)劑,采用包含表面張力的VOF模型進(jìn)行計算時,其結(jié)果需要謹(jǐn)慎對待。

Lafaurie等[7]采用基于VOF的新型算法,得出牛頓流體偽流動強度隨黏性的增大而降低。文獻(xiàn)[10]采用VOF和levelset模型進(jìn)行計算,結(jié)果表明Cas值幾乎為常數(shù),與黏性無關(guān)。K值體現(xiàn)了凝膠推進(jìn)劑的黏度平均水平,圖10為Cas值隨K值的變化規(guī)律。Cas隨K值幾乎成系數(shù)為0.05的比例增加。凝膠推進(jìn)劑具有較高的黏度,對Cas值貢獻(xiàn)較高。流場速度小,其差異未能體現(xiàn)出來。所以Cas值與凝膠推進(jìn)劑的平均黏度呈現(xiàn)出了明顯的相關(guān)性。

4 結(jié)束語

表面張力系數(shù)增大時,兩相界面處的偽流動速度增大,界面扭曲程度增大,但是偽流動強度值變化不大,因為偽流動強度定義中同樣包含了表面張力系數(shù)的影響,從而抵消表面張力系數(shù)的影響。

VOF模型下4種網(wǎng)格精度的Cas值變化趨勢和時間平均值基本一致,其時間平均值的標(biāo)準(zhǔn)方差僅為0.004 2,網(wǎng)格精度對偽流動的強度無影響。采用加密網(wǎng)格不能減小小偽流動帶來的數(shù)值誤差。

偽流動強度隨著凝膠推進(jìn)劑本構(gòu)參數(shù)n值的增加而增強,并且當(dāng)n>0.6時增加的梯度更大。偽流動強度最大值5.18比0.043要高出2個數(shù)量級。

偽流動強度隨凝膠推進(jìn)劑的平均黏度度量K值的增加而增強,并且?guī)缀醭上禂?shù)為0.05的比例增加。

[1] 劉晨.復(fù)雜燃燒流場數(shù)值模擬方法研究[D].南京:南京航空航天大學(xué),2009. LIU Chen.Numerical methods for complex combustion flowfields[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2009.(in Chinese)[2] 文華.基于CFD的柴油機噴霧混合過程的多維數(shù)值模擬[D].武漢:華中科技大學(xué),2004. WEN Hua.Multi-dimensional numerical modeling of spray mixing process in diesel engines based on CFD[D].Wuhan:Huazhong University of Science and Technology,2004.(in Chinese)

[3] MA D J.Atomization patterns and breakup characteristics of liquid sheets formed by two impinging jets:AIAA 2011-97[R].USA:AIAA,2011.

[4] ARIENTI M,LI X,SOTEROIU M C,et al.Coupled level-set/volume-of-fluid method for simulation of injector atomization[J].Journal of Propulsion & Power,2013,29(1):147-157.

[5] 強洪夫,劉虎,陳福振,等.基于SPH方法的射流撞擊仿真[J].推進(jìn)技術(shù),2012,33(3):424-429. QIANG Hong-fu,LIU Hu,CHEN Fu-zhen,et al.Simulation on jet impingement based on SPH method[J].Journal of Propulsion Technology,2012,33(3):424-429.(in Chinese)

[6] GUEYFFIER D,LI J,NADIM A,et al.Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[J].Journal of Computational Physics,1999,152(2):423-456.

[7] LAFAURIE B,NARDONE C,SCARDOVELLI R,et al.Modelling merging and fragmentation in multiphase flows with surfer[J].Journal of Computational Physics,1994,113(1):134-147.

[8] POPINET S,ZALESKI S.A front-tracking algorithm for accurate representation of surface tension[J].International Journal for Numerical Methods in Fluids,1999,30(6):775-793.

[9] ALEINOV I,PUCKETT E G.Computing surface tension with high-order kernels[C]//Proceeding of the sixth international symposium on computational fluid dynamics.Lake Tahoe:NV,1995:13-18.

[10] 馬東軍.可壓縮/不可壓縮流體交界面高精度數(shù)值方法的研究[D].合肥:中國科學(xué)技術(shù)大學(xué),2002. MA Dong-jun.Numerical investigation of the high resolution interfacial simulation methods applied to compressible and incompressible flows[D].Hefei:University of Science and Technology of China,2002.(in Chinese)

[11] 趙大勇,李維仲.VOF方法中幾種界面重構(gòu)技術(shù)的比較[J].熱科學(xué)與技術(shù),2003,2(4):318-323. ZHAO Da-yong,LI Wei-zhong.Comparison of three kinds of interface reconstruction methods applied to VOF[J].Journal of Thermal Science and Technology,2003,2(4):318-323.(in Chinese)

[12] RAHIMI S,NATAN B.Numerical solution of the flow of power-law gel propellants in converging injectors[J].Propellants Explosives Pyrotechnics,2000,25(4):203-212.

[13] RAESSI M.On modeling surface tension-dominant,large density ratio,two-phase flows[D].Toronto:University of Toronto,2008.

[14] 張德良.計算流體力學(xué)教程[M].北京:高等教育出版社,2010:367-374. ZHANG De-liang.Computational fluid dynamics tutorial[M].Beijing:High Education Press,2010:367-374.(in Chinese)

Research on Spurious Currents of Gel Propellant Based on VOF Model

CAO Qin-liu,FENG Feng,DENG Han-yu

(School of Mechanical Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

In order to study the characteristics of spurious currents considering surface tension in the two-phase flow for gel propellant primary-atomization,an incompressible fluid solver was developed based on synchronous overlap-add(SOLA)staggered mesh.Based on the gel propellant droplet,the spurious currents of VOF method was investigated.The results indicate that the numerical solution using this method agrees well with the analytical solution,and the relevant error is within 5%.Surface tension coefficient is taken into account in the definition of the strength of spurious currents,which partly eliminates the effect of the errors of surface tension of gel propellant on spurious currents values.The four mesh scales has nearly similar trends and time-average values of spurious currents,which reveals little effect of mesh scale on spurious currents.Spurious currents improves with the increase of the indexn,and the growth rate is greater whenn>0.6.The values of spurious currents linearly varies with the indexK,and the proportionality coefficient was 0.05.

gel propellant;surface tension;VOF;spurious currents

2016-12-22

航天科技創(chuàng)新基金(CASC03-02);中央高?；究蒲袠I(yè)務(wù)費專項基金(30920140112001);江蘇省普通高校學(xué)術(shù)學(xué)位研究生科研創(chuàng)新計劃(KYLX16_0474)

曹欽柳(1993- ),男,博士研究生,研究方向為航空宇航推進(jìn)理論與工程。E-mail:115101000190@njust.edu.cn。

V513

A

1004-499X(2017)02-0065-05