基于SPH方法的凝膠燃料單滴微爆過程模擬

強(qiáng)洪夫， 張林濤， 陳福振， 劉 虎， 石 超

(火箭軍工程大學(xué)動力工程系, 陜西 西安 710025)

1 引 言

凝膠燃料是一種新型的高性能燃料，凝膠燃料液滴在燃燒過程中呈現(xiàn)出與液體液滴燃燒過程不同的現(xiàn)象，隨著液滴內(nèi)部溫度逐漸升高，外部低沸點(diǎn)燃料達(dá)到沸點(diǎn)首先蒸發(fā)，并逐漸形成膠凝劑膜，液滴內(nèi)部出現(xiàn)沸騰蒸汽并聚合形成氣泡，氣泡逐漸聚合膨脹引起液滴產(chǎn)生周期性的膨脹、破裂和蒸汽噴射，稱為微爆過程。微爆過程是凝膠燃料燃燒一個(gè)非常重要的環(huán)節(jié)，國內(nèi)外學(xué)者在凝膠單滴燃燒實(shí)驗(yàn)中發(fā)現(xiàn)了凝膠的微爆現(xiàn)象。

Muller[1]在1997年進(jìn)行了煤油RP-1/Al 金屬凝膠液滴燃燒實(shí)驗(yàn); Solomon等[2]在2009年使用掛滴燃燒實(shí)驗(yàn)系統(tǒng)，研究了JP-8煤油基非金屬和金屬兩種凝膠液滴在靜止環(huán)境中的蒸發(fā)燃燒特性; Mishra[3]在2012年進(jìn)行了航空燃料(ATF)非金屬凝膠液滴的燃燒實(shí)驗(yàn)。上述研究表明: 有機(jī)凝膠液滴燃燒過程主要?dú)v經(jīng)表面低沸點(diǎn)推進(jìn)劑組分蒸發(fā)、彈性膠凝劑膜形成、膠凝劑膜周期性膨脹-破裂以及最終膠凝劑分解燃燒四個(gè)過程。國內(nèi)研究方面，Liu[4]和何博等[5]開展了有機(jī)非金屬凝膠偏二甲肼(UDMH)液滴在四氧化二氮(NTO)環(huán)境中著火燃燒的實(shí)驗(yàn)研究，研究了氧化劑濃度、溫度、壓力、對流速度和液滴初始尺寸等因素對其著火燃燒特性的影響。

部分國內(nèi)外學(xué)者使用數(shù)值模擬方法研究凝膠燃料的燃燒微爆問題。Kunin[7-8]在2005年建立了初步的非金屬凝膠液滴在靜止環(huán)境中的準(zhǔn)穩(wěn)態(tài)蒸發(fā)燃燒模型，并重點(diǎn)仿真分析了凝膠液滴初始尺寸對其表面膠凝劑膜厚度、膠凝劑膜首次形成時(shí)間及其首次破裂時(shí)間的影響。Muller等[9]建立了金屬凝膠液滴蒸發(fā)燃燒模型，并重點(diǎn)分析了金屬顆粒含量、金屬顆粒尺度對RP-1煤油/Al金屬有機(jī)凝膠液滴剛性外殼形成及其發(fā)生二次霧化的影響。何博[6]通過建立非穩(wěn)態(tài)蒸發(fā)燃燒模型，仿真分析了環(huán)境溫度、壓力對UDMH有機(jī)凝膠液滴燃燒的影響。

光滑粒子流體動力學(xué)(Smoothed Particles Hydrodynamics，SPH)方法[10]是一種純拉格朗日無網(wǎng)格粒子法，能夠更有效地捕捉凝膠液滴在微爆過程中具體的液滴形態(tài)的變化、追蹤液滴和氣泡界面。本研究探索性地運(yùn)用SPH方法對凝膠燃料單滴燃燒中微爆過程進(jìn)行了數(shù)值模擬。模擬中采用完全變光滑長度SPH方法，引入修正的JCD(Johnson-Cook Damage)強(qiáng)度模型和改進(jìn)的固壁邊界算法，并考慮粘性作用、界面?zhèn)鳠岷捅砻鎻埩Φ纫蛩?，仿真模擬了凝膠燃料單個(gè)液滴的微爆過程。

2 計(jì)算方法

2.1 控制方程

在計(jì)算中考慮熱傳導(dǎo)作用，考慮流體的粘性和表面張力作用，采用如下形式的Lagrange描述的流體力學(xué)控制方程[10]:

(1)

(2)

式中， p為壓力，Pa; ρ為密度,kg·m-3; v為粒子速度,m·s-1; F(v)為流體的粘性項(xiàng)，F(xiàn)(s)為表面張力項(xiàng)。

液滴在膨脹受熱時(shí)會吸收外界熱量，需要解決熱傳導(dǎo)問題，根據(jù)能量守恒和傅里葉定律可以得到微元中熱傳導(dǎo)方程為:

(3)

為計(jì)算方程(2)中的壓力項(xiàng)，需要引入狀態(tài)方程。因?yàn)樵诟邏涵h(huán)境中，燃料液滴內(nèi)部氣體不再滿足理想氣體狀態(tài)方程，Soave使用偏心因子對RK狀態(tài)方程進(jìn)行了非球形的常數(shù)修正，得到了SRK狀態(tài)方程[11]:

(4)

式中，根據(jù)文獻(xiàn)[11]選擇參數(shù)為:a=12.95，b=0.049，T=288 K。

本研究中，液相的計(jì)算采用了弱可壓縮(Weakly Compressible)狀態(tài)方程[12]。弱可壓縮狀態(tài)方程將所有理論上不可壓縮的流體看作為弱可壓縮，通過密度變化顯式求解流場壓力，狀態(tài)方程的形式為:

(5)

2.2 控制方程的SPH離散

在微爆過程中，由于氣泡膨脹會引起粒子間距增大，而傳統(tǒng)的SPH方法在粒子分布密度不一致時(shí)，不能得到與計(jì)算域一致的核近似精度。完全變光滑長度SPH法[13-15]修正了傳統(tǒng)SPH法中由于變光滑長度效應(yīng)所造成的計(jì)算偏差，該方法將光滑長度看作獨(dú)立的坐標(biāo)變量，并隨粒子的運(yùn)動而變化，即光滑長度同為空間和時(shí)間的函數(shù)hi=hri,t。為了使計(jì)算的鄰近粒子數(shù)目盡量保持不變，引入該約束條件，基于Lagrange控制方程組離散得到的完全變光滑方程組[14]為:

(6)

(7)

其中，光滑長度的動態(tài)變化為:

(8)

fi為修正系數(shù)由下式表示:

(9)

式中，Πij分別為人工粘度，其表達(dá)式參見文獻(xiàn)[10]，光滑長度變化率dhi/dt同密度變化率dρi/dt相互關(guān)聯(lián)，文中采用Leapfrog迭代的方法求解密度方程和光滑長度。

對熱傳導(dǎo)方程(3)式進(jìn)行SPH離散得到:

(10)

該熱傳導(dǎo)方程能夠在跨越材料界面時(shí)，自動滿足熱流守恒條件。

SPH方法的固壁邊界條件很難像網(wǎng)格法一樣嚴(yán)格實(shí)施。針對液滴與外殼間存在復(fù)雜滑移邊界的情況，為有效施加邊界條件，采用基于虛粒子法的新型固壁邊界施加模型[16]，可以有效地防止流體粒子對固壁邊界的穿透和邊界數(shù)值震蕩的情況。計(jì)算中核函數(shù)采用三次樣條函數(shù)[17]，時(shí)間步長推進(jìn)采用Leap-frog顯式形式[17]。

2.3 膠凝劑殼的本構(gòu)模型

對于膠凝劑外殼，無機(jī)凝膠液滴燃燒時(shí)會產(chǎn)生剛性外殼[9,18]，有機(jī)金屬凝膠液滴的膠凝劑外殼不同于前者，是一種彈性的連續(xù)性的固體結(jié)構(gòu)，可以認(rèn)為是彈塑性外殼[2]。Solomon[2]在實(shí)驗(yàn)中使用含有55%硼的JP-8煤油凝膠燃料，在燃燒時(shí)會產(chǎn)生具有彈塑性的凝膠殼。JCD模型是一種通用的彈塑性本構(gòu)模型[19]，采用的修正JCD本構(gòu)模型可以近似表征金屬有機(jī)凝膠液滴和無機(jī)凝膠液滴的膠凝劑殼，在計(jì)算膠凝劑外殼損傷及應(yīng)力變化，本構(gòu)模型表達(dá)式[19]為:

(11)

(12)

式中，pd是損傷閾值，pf是斷裂彈塑性應(yīng)變。Johnson和Cook[19]提出的損傷演化模型中將pf描述如下：

(13)

式中，D1～D5為材料常數(shù)，(3)式中選用常數(shù)見表1。

表1 本構(gòu)模型參數(shù)

Table 1 Parameters of the constitutive model

A/PaB/PaCDcD1D2D3D4D54.9×1058.07×1050.17390.0120.007051.732-0.054-0.01230

3 “微爆”過程的數(shù)值模擬及分析

實(shí)驗(yàn)研究表明，在凝膠燃料液滴微爆過程中，液滴內(nèi)部產(chǎn)生氣泡的大小和個(gè)數(shù)都是事先未知的，但總會存在一個(gè)較大氣泡率先沖破膠凝劑膜，產(chǎn)生微爆現(xiàn)象。因此，首先研究等效的單氣泡液滴，驗(yàn)證算法的有效性，再進(jìn)一步探索性研究多氣泡并存的情況。

3.1 模型構(gòu)建及求解

根據(jù)Solomon[2]實(shí)驗(yàn)結(jié)果，對JP-8金屬凝膠單滴燃燒實(shí)驗(yàn)的微爆過程進(jìn)行建模，初始建模結(jié)果如圖1所示，整個(gè)液滴直徑為1 mm，粒子總數(shù)7669個(gè)，其中流體粒子6225個(gè)，外殼粒子共有四層共1444個(gè)粒子。

圖1 微爆過程中SPH模型

Fig.1 The SPH model in micro explosion progress

Solomon[2]在環(huán)境溫度288 K，壓力為10 MPa的條件下，使用掛滴燃燒方法，進(jìn)行凝膠單滴燃燒實(shí)驗(yàn)研究。本研究根據(jù)Solomon的實(shí)驗(yàn)設(shè)置仿真模型參數(shù)，見表2。

表2 凝膠液滴微爆過程數(shù)值模擬的物性參數(shù)

Table 2 Physical parameters of the gelled propellant droplet numerical simulation

T0/Kcp/J·kg-1·K-1ρgel/kg·m-3Mgel/g·mol-1k/W·m-1·K-1288.0500015000.30.5

Note:T0is the initial temperature,cpis the specific heat capacity,ρgelis the density,Mgelis the molar mass,kis the heat transfer coefficient.

3.2 單氣泡燃燒“微爆”過程仿真與分析

凝膠液滴受熱膨脹的SPH仿真計(jì)算結(jié)果(圖2a)與Solomon實(shí)驗(yàn)結(jié)果(圖2b)[2]對比如圖2所示。由圖2可知，在燃燒初期隨著時(shí)間推移，氣泡逐漸受熱并不斷生長，由于膠凝劑膜的阻礙，蒸汽不能擴(kuò)散到環(huán)境中，導(dǎo)致液滴內(nèi)部壓力逐漸增加，液滴直徑逐漸增大。實(shí)驗(yàn)過程中，由于氣泡分布的隨機(jī)性導(dǎo)致液滴不規(guī)則膨脹，且液滴的懸掛方式導(dǎo)致液滴下部熱流密度更大，所以更容易產(chǎn)生氣泡。隨著氣泡體積的不斷增大，膠凝劑層不斷擴(kuò)大增加了拉伸應(yīng)力，當(dāng)達(dá)到它的彈性閾值時(shí)，進(jìn)一步增加壓力會產(chǎn)生破裂，導(dǎo)致“微爆”和蒸汽噴射的發(fā)生。從圖2可以看出，實(shí)驗(yàn)中液滴左下部位呈現(xiàn)透明狀，膠凝膜厚度很小，即將破裂，與700 ms時(shí)的仿真結(jié)果相似。在仿真計(jì)算800 ms時(shí)刻(見圖3)，可以發(fā)現(xiàn)膠凝劑膜殼破口增大，氣體噴出。

圖4是凝膠液滴在200 ms時(shí)的速度矢量圖，可以明顯看到氣泡生長膨脹過程中液滴內(nèi)部的速度分布狀況，由于氣泡生長和膨脹，引起了壓力升高，液滴膨脹。圖5是凝膠液滴燃燒過程液滴直徑變化圖，液滴隨著時(shí)間推移不斷增大，第一次膨脹時(shí)，直徑膨脹比D/D0可以達(dá)到1.4倍，實(shí)驗(yàn)中凝膠液滴的膨脹比大約為1.5倍。當(dāng)時(shí)間大約為0.75 s時(shí)，膠凝劑膜破裂，液滴體積迅速縮小。

膠凝劑膜損傷隨時(shí)間變化如圖6所示。由圖6可以看出，各時(shí)刻的膠凝劑膜的損傷系數(shù)變化情況，隨著氣泡生長膨脹，膠凝劑膜受到的拉伸應(yīng)力逐漸增大，且在臨近氣泡的膠凝劑膜的應(yīng)力最大，當(dāng)膠凝劑膜達(dá)到自身損傷極限時(shí)，損傷系數(shù)D達(dá)到1，膠凝劑膜破裂，氣體噴出，壓力迅速下降。

400 ms

700 ms

b. Solomon experimental results[2]

圖2 凝膠液滴受熱膨脹的SPH仿真計(jì)算與Solomon實(shí)驗(yàn)結(jié)果對比[2]

Fig.2 Comparison of the gelled droplet thermal expansion′s SPH simulative results and the Solomon′s experimental ones[2]

圖3 SPH仿真800 ms計(jì)算結(jié)果

Fig.3 SPH simulation results at 800 ms

圖4 200 ms時(shí)凝膠液滴速度矢量圖

Fig.4 A plot of velocity vector of the gelled droplet at 200 ms

圖5 凝膠液滴直徑隨時(shí)間變化的曲線

Fig.5 The curve of change in diameters of the gelled droplet with time

仿真結(jié)果可以得到凝膠液滴在燃燒過程中壓力變化分布和液滴內(nèi)部平均壓力變化(本研究所使用的壓力均指實(shí)際壓力與環(huán)境壓力的差值)。由圖7壓力分布圖及圖8中的單氣泡壓力震蕩曲線(黑色線)可以看出，由于氣泡膨脹使壓力升高并在液滴內(nèi)部傳播，壓力觸碰到膜殼時(shí)，反射回波，引起壓力震蕩上升。由于膜殼膨脹和能量耗散，壓力震蕩情況逐漸減弱，壓力整體緩慢上升。通過計(jì)算，氣泡沖破膜殼時(shí)液滴內(nèi)部的平均壓力約為150 kPa，氣體噴出之后壓力迅速減小，在100 ms內(nèi)壓力下降到68 kPa。

整個(gè)仿真過程，粒子秩序良好，無粒子穿透和數(shù)值不穩(wěn)定情況產(chǎn)生，使用完全變光滑長度SPH方法能夠很好地解決氣泡膨脹過程中粒子分布密度不一致的問題，同時(shí)新型固壁邊界方程和修正的JCD本構(gòu)方程也能夠有效地模擬膠凝劑膜應(yīng)力場的變化。

3.3 多氣泡“微爆”過程仿真

鑒于真實(shí)情況下，凝膠液滴內(nèi)部可能存在多個(gè)氣泡，本研究探索性設(shè)置大小不一的三個(gè)氣泡進(jìn)行多氣泡的“微爆”過程仿真，其中模型參數(shù)和粒子設(shè)置與上節(jié)相同。

a. 400 ms b. 700 ms c. 750 ms

圖6 膠凝劑膜損傷隨時(shí)間變化

Fig.6 Change in damage of the gelatinizer membrane with time

a. 10 ms b. 15 ms c. 40 ms

d. 90 ms e. 250 ms f. 700 ms

圖7 凝膠液滴壓力分布圖

Fig.7 Pressure distribution of the gelled droplet

圖8 液滴內(nèi)部平均壓力p0隨時(shí)間t變化曲線

Fig.8 The curves of change in the average pressure in the droplet with time

典型時(shí)刻的仿真計(jì)算的液滴變化如圖9所示。由圖9可見，大氣泡相對小氣泡膨脹更快，且大氣泡產(chǎn)生對壁面的壓力更大，導(dǎo)致膠凝劑膜在靠近大氣泡附近破裂。同時(shí)，破裂時(shí)間相較于上節(jié)單一氣泡的750 ms縮短至550 ms，這是由于液滴內(nèi)的氣液比更大，汽化核多，相變速度快，導(dǎo)致壓力上升的更為迅速。

多氣泡凝膠液滴壓力分布如圖10所示。由圖10可以發(fā)現(xiàn)，在90 ms時(shí)液滴內(nèi)壓力分布已經(jīng)較為均勻，由于存在多個(gè)壓力中心，導(dǎo)致壓力更快地達(dá)到較為均勻狀態(tài)。同樣也可以根據(jù)圖8中壓力震蕩曲線看到平均壓力與單氣泡的震蕩形式相近，但震蕩幅度更小，同時(shí)壓力上升速度更快，破裂時(shí)壓力為145 kPa，比單個(gè)氣泡略低，因?yàn)樵诖藭r(shí)刻多氣泡液滴的體積大于單氣泡，導(dǎo)致膠凝劑膜更加脆弱; 同時(shí)，破裂后壓力降為71.5 kPa，比單個(gè)氣泡略高，因?yàn)樵诖饲闆r下，中央大氣泡較單一氣泡略小，噴射氣體減少，引起壓力降低幅度小于單個(gè)氣泡。

a. 200 ms b. 550 ms c. 650 ms

圖9 多氣泡液滴受熱膨脹仿真計(jì)算

Fig.9 Simulation of thermal expansion of the gelled droplet with multiple bubbles

a. 10 ms b. 90 ms c. 250 ms

圖10 多氣泡凝膠液滴壓力分布圖

Fig.10 Pressure distribution of the gelled droplet with multiple bubbles

4 結(jié) 論

采用完全變光滑長度SPH方法，對凝膠燃料液滴在微爆過程中的變化情況進(jìn)行了數(shù)值模擬，得到了凝膠液滴膨脹、破裂以及氣體噴出的具體過程及形態(tài)變化情況，分析了單氣泡和多氣泡凝膠液滴在燃燒過程中，液滴半徑、壓力分布以及速度等物理量的變化規(guī)律。數(shù)值模擬結(jié)果顯示: 在凝膠液滴微爆過程中，隨著時(shí)間的推移，體積逐漸膨脹，壓力呈現(xiàn)震蕩上升形式，當(dāng)直徑膨脹比D/D0達(dá)到1.4時(shí)，膠凝劑膜破裂，液滴體積迅速縮小，壓力急劇下降; 氣泡沖破膜殼時(shí)液滴內(nèi)部的平均壓力約為150 kPa，氣體噴出之后壓力迅速減小，在100 ms內(nèi)壓力下降到68 kPa; 對于多氣泡凝膠液滴，大氣泡相對小氣泡膨脹更快，且大氣泡產(chǎn)生對壁面的壓力更大，導(dǎo)致膠凝劑膜在靠近大氣泡附近破裂，同時(shí)，破裂時(shí)間相較于單一氣泡的750 ms縮短至550 ms。本文的數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果相吻合，揭示了凝膠液滴微爆過程的內(nèi)在機(jī)理，驗(yàn)證了SPH算法在解決此類問題上的有效性。后續(xù)工作可結(jié)合化學(xué)反應(yīng)變化對凝膠液滴完整的燃燒過程進(jìn)行模擬計(jì)算。

參考文獻(xiàn):

[1] Muller D C. Aluminum/Hydrocarbon Gel Propellants: Anexperimental and theoretical investigation of secondary atomization and predicted rocket engine performance[D]. Pennsylvania: The Pennsylvania State University, 1997.

[2] Solomon Y, Natan B, Cohen Y. Combustion of gel fuels based on organic gellants[J].CombustionandFlame, 2009, 156(1): 261-268.

[3] Mishra D P, Patyal A. Effects of initial droplet diameter and pressure on burning of ATF gel propellant droplets[J].Fuel, 2012, 95: 226-233.

[4] Liu Z, Hu X, He Z, et al. Experimental study on the combustion and micro explosion of freely falling gelled unsymmetrical dimethylhydrazine(UDMH) fuel droplets[J].Energies, 2012, 5(12): 3126-3136.

[5] He B, Nie W, Feng S, et al. Effects of NTO oxidizer temperature and pressure on hypergolic ignition delay and life time of UDMH organic gel droplet[J].Propellants,Explosives,Pyrotechnics, 2013, 38(5): 665-684.

[6] He B, Nie W, He H. Unsteady combustion model of nonmetalized organic gel fuel droplet[J].Energy&Fuels, 2012: 1145274908.

[7] Kunin A, Greenberg J, Natan B. A simple phenomenological model of an organic gel spray diffusion flame[J].CombustionScienceandTechnology, 2007, 180(1): 27-44.

[8] Kunin A, Natan B, Greenberg B. Theoretical model of the transient combustion of organic-gellant-based gel fuel droplets[J].JournalofPropulsionandPower, 2010, 26: 765-771.

[9] Muller D C,Turns S R. Some aspects secondary atomization of aluminum/hydrocarbon slurry propellant[J].JournalofHeatTransfer, 2003, 125: 535-537.

[10] Monaghan J J. Smoothed particle hydrodynamics[R]. Progress in Physics, 2005.

[11] 朱自強(qiáng). 化工熱力學(xué)[M]. 北京: 化學(xué)工業(yè)出版社, 1987: 15-111.

[12] Monaghan J J. Simulating free surface flows with SPH[J].JournalofComputationalPhysics, 1994, 110(2): 399-406.

[13] QIANG Hong-fu, GAO Wei-ran. A new SPH equation including variable smoothing lengths aspects and its implementation[C]∥ISCM, 2007, Beijing.

[14] 強(qiáng)洪夫, 高巍然. 修正變光滑長度SPH方法及其應(yīng)用[J]. 解放軍理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2007, 8(5): 419-424.

QIANG Hong-fu, GAO Wei-ran. Modified SPH method considering full variable smoothing lengths effects and its applications[J].JournalofPLAUniversityofScienceandTechnology, 2007, 8(5): 419-424.

[15] Monaghan J J. SPH without a tensile instability[J].JournalofComputationalPhysics, 2000, 159(2): 290-311.

[16] 劉虎, 強(qiáng)洪夫, 陳福振, 等. 一種新型光滑粒子動力學(xué)固壁邊界施加模型[J]. 物理學(xué)報(bào), 2015, 64(9): 094701.

LIU Hu, QIANG Hong-fu, CHEN Fu-zhen, et al. A new boundary treatment method in smoothed particle hydrodynamics[J].ActaPhysSin, 2015, 64(9): 094701.

[17] 劉桂榮, 劉謀斌. 光滑粒子流體動力學(xué)-一種無網(wǎng)格粒子法[M]. 韓旭, 楊剛, 強(qiáng)洪夫, 譯. 長沙: 湖南大學(xué)出版社，2005.

LIU Giu-rong, LIU Mou-bin.Smoothed particle hydrodynamics: a meshfree particle method. HAN Xu, YANG Gang, QIANG Hong-fu, Translated. Changsha: Hunan University Press,2005.

[18] Mueller D C, Turns S R. Ignition and combustion characteristics of metallized propellants-phase Ⅱ[R]. Annual Report, 1993.

[19] Johnson G R, Cook W H. Fracture characteristics of three metals subjected to various strains, strain rates, temperatures, and pressures[J].EngineeringFractureMechanics, 1985, 21(1): 31-48.

[20] Lemaitre J. A course on damage mechanics[M]. Berlin: Springer, 1996.

[21] Camacho G T, Ortiz M. Adaptive lagrangian modelling of ballistic penetration of metallic targets[J].ComputeMethodsinAppliedMechanicsandEngineering, 1997, 142: 269-301.