曹欽柳,封 鋒,鄧寒玉

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

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

曹欽柳,封 鋒,鄧寒玉

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

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

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

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

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

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

1 數(shù)值方法

1.1 計(jì)算模型

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

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

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

1.3 黏性模型

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

(1)

(2)

(3)

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

(4)

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

(5)

(6)

(7)

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

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

(8)

1.4 表面張力模型

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

fs=σκ(F)

(9)

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

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

Cas=|U|μ/σ

(10)

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

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

2 數(shù)值驗(yàn)證

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

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

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

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

Lafaurie等[7]基于VOF方法對(duì)牛頓流體氣泡的偽流動(dòng)進(jìn)行了研究,指出表面張力增大時(shí)偽流動(dòng)越明顯。表面張力系數(shù)是表面張力中的一個(gè)重要參數(shù)。為研究偽流動(dòng)強(qiáng)弱隨表面張力系數(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,對(duì)其進(jìn)行數(shù)值計(jì)算。

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)束語(yǔ)

表面張力系數(shù)增大時(shí),兩相界面處的偽流動(dòng)速度增大,界面扭曲程度增大,但是偽流動(dòng)強(qiáng)度值變化不大,因?yàn)閭瘟鲃?dòng)強(qiáng)度定義中同樣包含了表面張力系數(shù)的影響,從而抵消表面張力系數(shù)的影響。

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

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

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

[1] 劉晨.復(fù)雜燃燒流場(chǎng)數(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的柴油機(jī)噴霧混合過(guò)程的多維數(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] 強(qiáng)洪夫,劉虎,陳福振,等.基于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].合肥:中國(guó)科學(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] 張德良.計(jì)算流體力學(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ù)費(fèi)專項(xiàng)基金(30920140112001);江蘇省普通高校學(xué)術(shù)學(xué)位研究生科研創(chuàng)新計(jì)劃(KYLX16_0474)

曹欽柳(1993- ),男,博士研究生,研究方向?yàn)楹娇沼詈酵七M(jìn)理論與工程。E-mail:115101000190@njust.edu.cn。

V513

A

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