基于VOF模型的凝膠推進劑的偽流動研究

曹欽柳,封 鋒,鄧寒玉

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

基于VOF模型的凝膠推進劑的偽流動研究

曹欽柳,封 鋒,鄧寒玉

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

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

凝膠推進劑;表面張力;流體體積函數;偽流動

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

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

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

1 數值方法

1.1 計算模型

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

1.2 界面重構方法

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

1.3 黏性模型

凝膠推進劑在流動上表現出非牛頓流體的性質,其流變本構方程可寫為如下形式[12]:

(1)

(2)

(3)

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

(4)

由上述公式,可以得到應力張量各分量表達式為

(5)

(6)

(7)

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

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

(8)

1.4 表面張力模型

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

fs=σκ(F)

(9)

式中:σ為表面張力系數,κ為兩相界面的曲率,F為網格中液體相的體積分數。上式求得的表面張力實際是一種壓力對坐標的變化量,并作用于網格內所有的流體,它作為源項添加到動量方程中。

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

Cas=|U|μ/σ

(10)

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

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

2 數值驗證

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

3 計算結果與分析

本文為深入探討VOF模型中凝膠液滴的偽流動現象,對不同表面張力系數,網格尺度和冪率模型本構參數下凝膠液滴的偽流動進行了數值計算。空氣密度取1.25kg/m3,凝膠推進劑的密度在不同的配比下,其值有所差異,本文直接取為1 000kg/m3。按照冪率模型,當剪切速率為0時,其表觀黏度趨向于無窮大,因此數值計算時需要對上界值進行限定,本文取最大表觀黏度值ηmax為10Pa·s。

3.1 表面張力系數的影響

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

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

3.2 網格尺度的影響

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

表1 網格參數

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

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

3.3 本構參數的影響

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

表2 凝膠推進劑的本構參數

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

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

4 結束語

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

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

偽流動強度隨著凝膠推進劑本構參數n值的增加而增強,并且當n>0.6時增加的梯度更大。偽流動強度最大值5.18比0.043要高出2個數量級。

偽流動強度隨凝膠推進劑的平均黏度度量K值的增加而增強,并且幾乎成系數為0.05的比例增加。

[1] 劉晨.復雜燃燒流場數值模擬方法研究[D].南京:南京航空航天大學,2009. LIU Chen.Numerical methods for complex combustion flowfields[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2009.(in Chinese)[2] 文華.基于CFD的柴油機噴霧混合過程的多維數值模擬[D].武漢:華中科技大學,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].推進技術,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] 馬東軍.可壓縮/不可壓縮流體交界面高精度數值方法的研究[D].合肥:中國科學技術大學,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方法中幾種界面重構技術的比較[J].熱科學與技術,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] 張德良.計算流體力學教程[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

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

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

V513

A

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