1 引言

凝胶燃料是一种新型的高性能燃料，凝胶燃料液滴在燃烧过程中呈现出与液体液滴燃烧过程不同的现象，随着液滴内部温度逐渐升高，外部低沸点燃料达到沸点首先蒸发，并逐渐形成胶凝剂膜，液滴内部出现沸腾蒸汽并聚合形成气泡，气泡逐渐聚合膨胀引起液滴产生周期性的膨胀、破裂和蒸汽喷射，称为微爆过程。微爆过程是凝胶燃料燃烧一个非常重要的环节，国内外学者在凝胶单滴燃烧实验中发现了凝胶的微爆现象。

Muller^[1]在1997年进行了煤油RP-1/Al金属凝胶液滴燃烧实验; Solomon等^[2]在2009年使用挂滴燃烧实验系统，研究了JP-8煤油基非金属和金属两种凝胶液滴在静止环境中的蒸发燃烧特性; Mishra^[3]在2012年进行了航空燃料(ATF)非金属凝胶液滴的燃烧实验。上述研究表明:有机凝胶液滴燃烧过程主要历经表面低沸点推进剂组分蒸发、弹性胶凝剂膜形成、胶凝剂膜周期性膨胀-破裂以及最终胶凝剂分解燃烧四个过程。国内研究方面，Liu^[4]和何博等^[5]开展了有机非金属凝胶偏二甲肼(UDMH)液滴在四氧化二氮(NTO)环境中着火燃烧的实验研究，研究了氧化剂浓度、温度、压力、对流速度和液滴初始尺寸等因素对其着火燃烧特性的影响。

部分国内外学者使用数值模拟方法研究凝胶燃料的燃烧微爆问题。Kunin^[7-8]在2005年建立了初步的非金属凝胶液滴在静止环境中的准稳态蒸发燃烧模型，并重点仿真分析了凝胶液滴初始尺寸对其表面胶凝剂膜厚度、胶凝剂膜首次形成时间及其首次破裂时间的影响。Muller等^[9]建立了金属凝胶液滴蒸发燃烧模型，并重点分析了金属颗粒含量、金属颗粒尺度对RP-1煤油/Al金属有机凝胶液滴刚性外壳形成及其发生二次雾化的影响。何博^[6]通过建立非稳态蒸发燃烧模型，仿真分析了环境温度、压力对UDMH有机凝胶液滴燃烧的影响。

光滑粒子流体动力学(Smoothed Particles Hydrodynamics，SPH)方法^[10]是一种纯拉格朗日无网格粒子法，能够更有效地捕捉凝胶液滴在微爆过程中具体的液滴形态的变化、追踪液滴和气泡界面。本研究探索性地运用SPH方法对凝胶燃料单滴燃烧中微爆过程进行了数值模拟。模拟中采用完全变光滑长度SPH方法，引入修正的JCD(Johnson-Cook Damage)强度模型和改进的固壁边界算法，并考虑粘性作用、界面传热和表面张力等因素，仿真模拟了凝胶燃料单个液滴的微爆过程。

2 计算方法 2.1 控制方程

在计算中考虑热传导作用，考虑流体的粘性和表面张力作用，采用如下形式的Lagrange描述的流体力学控制方程^[10]:

$\frac{{{\rm{d}}\rho }}{{{\rm{d}}t}} = - \rho \nabla \cdot v$

(1)

$\frac{{{\rm{d}}v}}{{{\rm{d}}t}} = - \frac{1}{\rho }\nabla \rho + {\mathit{\boldsymbol{F}}^{\left( v \right)}} + {\mathit{\boldsymbol{F}}^{\left( s \right)}}$

(2)

式中，p为压力，Pa; ρ为密度, kg·m^-3; v为粒子速度, m·s^-1; F ^(v)为流体的粘性项，F ^(s)为表面张力项。

液滴在膨胀受热时会吸收外界热量，需要解决热传导问题，根据能量守恒和傅里叶定律可以得到微元中热传导方程为:

$\rho \frac{{{\rm{d}}E}}{{{\rm{d}}t}}\nabla \cdot \left( {k\nabla T} \right) + \dot Q$

(3)

式中，dE=c_p·ΔT，c_p为材料的比热容, J·kg^-1·K^-1，${\dot Q}$为材料的内热源。

为计算方程(2) 中的压力项，需要引入状态方程。因为在高压环境中，燃料液滴内部气体不再满足理想气体状态方程，Soave使用偏心因子对RK状态方程进行了非球形的常数修正，得到了SRK状态方程^[11]:

$p = \frac{{RT}}{{V - b}} - \frac{a}{{{T^{\frac{1}{2}}}V\left( {V + b} \right)}}$

(4)

式中，根据文献[11]选择参数为: a=12.95，b=0.049，T=288 K。

本研究中，液相的计算采用了弱可压缩(Weakly Compressible)状态方程^[12]。弱可压缩状态方程将所有理论上不可压缩的流体看作为弱可压缩，通过密度变化显式求解流场压力，状态方程的形式为:

$p = {p_0}\left[ {{{\left( {\frac{\rho }{{{\rho _0}}}} \right)}^\gamma } - 1} \right]$

(5)

式中，ρ₀为流体的初始密度, kg·m^-3; p₀为参考压强, Pa; γ是常数，一般取γ=7，参考压强p₀=c_s²ρ/γ，c_s为流场中的声速, m·s^-1，一般选取为流场最大速度的10倍左右。

2.2 控制方程的SPH离散

在微爆过程中，由于气泡膨胀会引起粒子间距增大，而传统的SPH方法在粒子分布密度不一致时，不能得到与计算域一致的核近似精度。完全变光滑长度SPH法^[13-15]修正了传统SPH法中由于变光滑长度效应所造成的计算偏差，该方法将光滑长度看作独立的坐标变量，并随粒子的运动而变化，即光滑长度同为空间和时间的函数h_i=h(r_i, t)。为了使计算的邻近粒子数目尽量保持不变，引入该约束条件，基于Lagrange控制方程组离散得到的完全变光滑方程组^[14]为:

$\frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}} = \sum\limits_{j = 1}^N {{m_j}} \left[ {{V_{ij}} \cdot {\nabla _i}{W_{ij}} + \frac{1}{2}\left( {\frac{{{\rm{d}}{h_i}}}{{{\rm{d}}t}} + \frac{{{\rm{d}}{h_j}}}{{{\rm{d}}t}}} \right)\frac{{\partial {W_{ij}}}}{{\partial h}}} \right]$

(6)

$\frac{{{\rm{d}}{V_i}}}{{{\rm{d}}t}} = - \sum\limits_{j = 1}^N {{m_j}\left( {{f_i}\frac{{{P_i}}}{{\rho _i^2}} + {f_j}\frac{{{P_j}}}{{\rho _j^2}} + {\prod _{ij}}} \right)} {\nabla _i}{W_{ij}}$

(7)

其中，光滑长度的动态变化为:

$\frac{{{\rm{d}}{h_i}}}{{{\rm{d}}t}} = - \frac{1}{{\rm{d}}}\frac{{{h_i}}}{{{\rho _i}}}\frac{{{\rm{d}}{\rho _i}}}{{{\rm{d}}t}}$

(8)

f_i为修正系数由下式表示:

${f_i} = {\left( {1 + \frac{{{h_i}}}{{{\rm{d}}{\rho _i}}}\sum\limits_{j = 1}^N {{m_j}} \frac{{\partial {W_{ij}}}}{{\partial {h_i}}}} \right)^{ - 1}}$

(9)

式中，Π_ij分别为人工粘度，其表达式参见文献[10]，光滑长度变化率dh_i/dt同密度变化率dρ_i/dt相互关联，文中采用Leapfrog迭代的方法求解密度方程和光滑长度。

对热传导方程(3) 式进行SPH离散得到:

$ - \nabla \cdot q = \sum\limits_{j = 1}^n {\frac{{{m_j}}}{{{\rho _i}}}} \frac{{4{k_1}{k_2}}}{{{k_1} + {k_2}}}\left( {{T_j} - {T_i}} \right)\frac{{\left( {{{\vec r}_j} - {{\vec r}_i}} \right)}}{{{{\left| {{{\vec r}_j} - {{\vec r}_i}} \right|}^2}}} \cdot \nabla {W_{ij}}$

(10)

该热传导方程能够在跨越材料界面时，自动满足热流守恒条件。

SPH方法的固壁边界条件很难像网格法一样严格实施。针对液滴与外壳间存在复杂滑移边界的情况，为有效施加边界条件，采用基于虚粒子法的新型固壁边界施加模型^[16]，可以有效地防止流体粒子对固壁边界的穿透和边界数值震荡的情况。计算中核函数采用三次样条函数^[17]，时间步长推进采用Leap-frog显式形式^[17]。

2.3 胶凝剂壳的本构模型

对于胶凝剂外壳，无机凝胶液滴燃烧时会产生刚性外壳^{[9, 18]}，有机金属凝胶液滴的胶凝剂外壳不同于前者，是一种弹性的连续性的固体结构，可以认为是弹塑性外壳^[2]。Solomon^[2]在实验中使用含有55%硼的JP-8煤油凝胶燃料，在燃烧时会产生具有弹塑性的凝胶壳。JCD模型是一种通用的弹塑性本构模型^[19]，采用的修正JCD本构模型可以近似表征金属有机凝胶液滴和无机凝胶液滴的胶凝剂壳，在计算胶凝剂外壳损伤及应力变化，本构模型表达式^[19]为:

${\sigma _{{\rm{eq}}}} = \left( {1 - D} \right)\left[ {A + B{r^n}} \right]\left[ {1 + C\ln \;{{\dot r}^*}} \right]\left[ {1 - {T^{*m}}} \right]$

(11)

式中，r表示累积损伤塑性应变，$r = \left( {1 - D} \right)p,{{\dot r}^*} = \left( {1 - D} \right){{\dot p}^*}$。D为损伤变量^[20]，D=0表示材料没有损伤，D=1表示材料完全失效。当材料出现宏裂纹时，损伤变量的临界值小于1，失效准则可描述为D=D_C≤1，其中D_C是损伤变量临界值^[21]。损伤变量D是累积弹性应变p的函数，可描述如下：

$\dot D = \left\{ {\begin{array}{*{20}{c}} 0&{p \le {p_{\rm{d}}}}\\ {\frac{{{D_{\rm{C}}}}}{{{p_{\rm{f}}} - {p_{\rm{d}}}}}\dot p}&{p \ge {p_{\rm{d}}}} \end{array}} \right.$

(12)

式中，p_d是损伤阈值，p_f是断裂弹塑性应变。Johnson和Cook^[19]提出的损伤演化模型中将p_f描述如下：

${p_{\rm{f}}} = \left[ {{D_1} + {D_2}\exp \left( {{D_3}{\sigma ^*}} \right)} \right]\left[ {1 + {D_4}\ln {{\dot p}^*}} \right]\left[ {1 + {D_5}{T^*}} \right]$

(13)

式中，D₁~D₅为材料常数，(3) 式中选用常数见表 1。

3 “微爆”过程的数值模拟及分析

实验研究表明，在凝胶燃料液滴微爆过程中，液滴内部产生气泡的大小和个数都是事先未知的，但总会存在一个较大气泡率先冲破胶凝剂膜，产生微爆现象。因此，首先研究等效的单气泡液滴，验证算法的有效性，再进一步探索性研究多气泡并存的情况。

3.1 模型构建及求解

根据Solomon^[2]实验结果，对JP-8金属凝胶单滴燃烧实验的微爆过程进行建模，初始建模结果如图 1所示，整个液滴直径为1 mm，粒子总数7669个，其中流体粒子6225个，外壳粒子共有四层共1444个粒子。

Solomon^[2]在环境温度288 K，压力为10 MPa的条件下，使用挂滴燃烧方法，进行凝胶单滴燃烧实验研究。本研究根据Solomon的实验设置仿真模型参数，见表 2。

3.2 单气泡燃烧“微爆”过程仿真与分析

凝胶液滴受热膨胀的SPH仿真计算结果(图 2a)与Solomon实验结果(图 2b)^[2]对比如图 2所示。由图 2可知，在燃烧初期随着时间推移，气泡逐渐受热并不断生长，由于胶凝剂膜的阻碍，蒸汽不能扩散到环境中，导致液滴内部压力逐渐增加，液滴直径逐渐增大。实验过程中，由于气泡分布的随机性导致液滴不规则膨胀，且液滴的悬挂方式导致液滴下部热流密度更大，所以更容易产生气泡。随着气泡体积的不断增大，胶凝剂层不断扩大增加了拉伸应力，当达到它的弹性阈值时，进一步增加压力会产生破裂，导致“微爆”和蒸汽喷射的发生。从图 2可以看出，实验中液滴左下部位呈现透明状，胶凝膜厚度很小，即将破裂，与700 ms时的仿真结果相似。在仿真计算800 ms时刻(见图 3)，可以发现胶凝剂膜壳破口增大，气体喷出。

图 4是凝胶液滴在200 ms时的速度矢量图，可以明显看到气泡生长膨胀过程中液滴内部的速度分布状况，由于气泡生长和膨胀，引起了压力升高，液滴膨胀。图 5是凝胶液滴燃烧过程液滴直径变化图，液滴随着时间推移不断增大，第一次膨胀时，直径膨胀比D/D₀可以达到1.4倍，实验中凝胶液滴的膨胀比大约为1.5倍。当时间大约为0.75 s时，胶凝剂膜破裂，液滴体积迅速缩小。

胶凝剂膜损伤随时间变化如图 6所示。由图 6可以看出，各时刻的胶凝剂膜的损伤系数变化情况，随着气泡生长膨胀，胶凝剂膜受到的拉伸应力逐渐增大，且在临近气泡的胶凝剂膜的应力最大，当胶凝剂膜达到自身损伤极限时，损伤系数D达到1，胶凝剂膜破裂，气体喷出，压力迅速下降。

仿真结果可以得到凝胶液滴在燃烧过程中压力变化分布和液滴内部平均压力变化(本研究所使用的压力均指实际压力与环境压力的差值)。由图 7压力分布图及图 8中的单气泡压力震荡曲线(黑色线)可以看出，由于气泡膨胀使压力升高并在液滴内部传播，压力触碰到膜壳时，反射回波，引起压力震荡上升。由于膜壳膨胀和能量耗散，压力震荡情况逐渐减弱，压力整体缓慢上升。通过计算，气泡冲破膜壳时液滴内部的平均压力约为150 kPa，气体喷出之后压力迅速减小，在100 ms内压力下降到68 kPa。

整个仿真过程，粒子秩序良好，无粒子穿透和数值不稳定情况产生，使用完全变光滑长度SPH方法能够很好地解决气泡膨胀过程中粒子分布密度不一致的问题，同时新型固壁边界方程和修正的JCD本构方程也能够有效地模拟胶凝剂膜应力场的变化。

3.3 多气泡“微爆”过程仿真

鉴于真实情况下，凝胶液滴内部可能存在多个气泡，本研究探索性设置大小不一的三个气泡进行多气泡的“微爆”过程仿真，其中模型参数和粒子设置与上节相同。

典型时刻的仿真计算的液滴变化如图 9所示。由图 9可见，大气泡相对小气泡膨胀更快，且大气泡产生对壁面的压力更大，导致胶凝剂膜在靠近大气泡附近破裂。同时，破裂时间相较于上节单一气泡的750 ms缩短至550 ms，这是由于液滴内的气液比更大，汽化核多，相变速度快，导致压力上升的更为迅速。

多气泡凝胶液滴压力分布如图 10所示。由图 10可以发现，在90 ms时液滴内压力分布已经较为均匀，由于存在多个压力中心，导致压力更快地达到较为均匀状态。同样也可以根据图 8中压力震荡曲线看到平均压力与单气泡的震荡形式相近，但震荡幅度更小，同时压力上升速度更快，破裂时压力为145 kPa，比单个气泡略低，因为在此时刻多气泡液滴的体积大于单气泡，导致胶凝剂膜更加脆弱; 同时，破裂后压力降为71.5 kPa，比单个气泡略高，因为在此情况下，中央大气泡较单一气泡略小，喷射气体减少，引起压力降低幅度小于单个气泡。

4 结论

采用完全变光滑长度SPH方法，对凝胶燃料液滴在微爆过程中的变化情况进行了数值模拟，得到了凝胶液滴膨胀、破裂以及气体喷出的具体过程及形态变化情况，分析了单气泡和多气泡凝胶液滴在燃烧过程中，液滴半径、压力分布以及速度等物理量的变化规律。数值模拟结果显示:在凝胶液滴微爆过程中，随着时间的推移，体积逐渐膨胀，压力呈现震荡上升形式，当直径膨胀比D/D₀达到1.4时，胶凝剂膜破裂，液滴体积迅速缩小，压力急剧下降; 气泡冲破膜壳时液滴内部的平均压力约为150 kPa，气体喷出之后压力迅速减小，在100 ms内压力下降到68 kPa; 对于多气泡凝胶液滴，大气泡相对小气泡膨胀更快，且大气泡产生对壁面的压力更大，导致胶凝剂膜在靠近大气泡附近破裂，同时，破裂时间相较于单一气泡的750 ms缩短至550 ms。本文的数值模拟结果与实验结果相吻合，揭示了凝胶液滴微爆过程的内在机理，验证了SPH算法在解决此类问题上的有效性。后续工作可结合化学反应变化对凝胶液滴完整的燃烧过程进行模拟计算。