1 引言

小推力液体火箭发动机通常采用的冷却方式是边区冷却喷嘴射流撞击到壁面形成液膜继而进行冷却^[1]。在数值仿真工作中，准确地模拟出液膜冷却的效果会为发动机设计制造提供极大的便利。因此准确的液膜形态和合适的流体/结构耦合传热方法在研究该类问题的数值仿真显得尤为重要。

国内外对该类发动机的数值仿真工作有很多，有些描述了燃烧室中的流动和传热，但没有描述液膜形态^[2~6]；有些描述液膜行为，但没有描述与结构间的耦合传热方法^[7~9]。故提及一套精细的液膜方程和一种完备的考虑液膜条件下的耦合传热方法对相关数值仿真工作有较重大的意义。

目前多数研究都采用基于Euler-Lagrange的方法来对发动机的燃烧流场进行求解。对于液膜的求解，有采用简单的基于经验公式的求解方法，有采用基于Euler思想的VOF方法^[10]，少有提及基于Lagrange求解液膜的方法。由于多数喷雾模型都采用离散相的求解思路，如果采用Euler方法计算液膜的话势必要引入复杂的壁面/液滴作用机制，用于把Lagrange液滴转化为Euler液膜。故在采用离散相模型对气相场中液滴追踪模拟的条件下，Lagrange液膜模型在程序上的处理更方便。

Lagrange模型是把成片的液膜模拟成一个个离散相的液膜液滴，针对各个液膜液滴分别建立运动、蒸发、传热方程，统计各个液膜液滴的行为进而得出整片液膜的行为。目前多数Lagrange液膜模型的研究都是在内燃机的工作条件下开展的。O’rourke和Amsden^{[11, 12]}率先提出一套Lagrange液膜运动与传热模型（以下简称OA模型），通过内燃机的仿真实验对比得出该模型既可与离散相喷注模型保持较好的兼容性，又可以在壁面网格较粗糙的情况下保证计算结果有较高的精度。Han等^[13]指出对照Mathews所做的平板液膜实验，OA模型的仿真结果存在液膜前端延展过长的问题，Han通过修正碰撞液滴的速度分布得到了较理想的液膜形态。Zhang^{[14, 15]}在Han工作的基础上，考虑了液滴撞击液膜的影响，对OA模型做了进一步的完善。上述研究均用于指导内燃机的仿真工作^{[16, 17]}，在液体火箭发动机中对Lagrange模型的应用较少。

本文先搭建耦合传热方法，在此耦合传热方法的基础上建立一套适合小推力液体火箭发动机的Lagrange液膜模型，进而进行数值仿真工作。

2 物理模型和计算方法 2.1 耦合传热方法

本文采用松耦合的分区求解方法搭建流场/结构传热耦合计算平台，热流有三种传递走向：

（a）热流由气相场传给结构，结构传给空间。

（b）热流由气相场传给液膜，液膜传给结构，结构传给空间。

（c）气相场把热流传给液膜，结构把一部分热流传给液膜，把另一部分热流传到空间。

当前壁面未附着液膜，气相直接加热燃烧室内壁面，热流从内壁面传到外壁面继而向空间辐射，此种情况对应热流走向（a）。当前壁面附着液膜，气相加热液膜。如果内壁面温度低于液膜平均温度，气相场的热流从液膜传到内壁面，对应热流走向（b）。如果内壁面温度高于液膜平均温度，液膜除了从气相场吸收热流，还要从壁面吸收热流，对应热流走向（c）。流场计算收敛后会分别统计气相传给壁面的热流和液膜传给（吸收）壁面的热流，壁面接收热流信息进行结构导热计算，返给流场温度分布作为下一次迭代的流场边界。

低温液膜在相对高温的气体环境中，首先会经历一个不平衡蒸发阶段，同时液膜温度不断上升；当液膜温度上升到某一个值后，进入液膜的热量刚好用于支付液体的蒸发潜热以及将蒸汽加热到环境气体的温度，此时液膜温度不会继续升高，称为温度平衡阶段^[18]。当内壁面温度低于液膜平衡温度时，壁面接收液膜传递的热流信息。当壁面温度高于平衡温度时，壁面要把液膜平衡温度作为温度边界条件，返给液膜热流用于液膜的蒸发，以上计算流程如图 1所示。

2.2 Lagrange液膜模型 2.2.1 控制方程

在每一步都要统计壁面网格上的液膜液滴体积进而得出液膜厚度，所以液膜的质量在Lagrange液膜模型中是隐式守恒的。

壁面上的液膜的流动和传热蒸发等行为受到上方的气相场、下方的壁面和撞击到其上的射流的影响。在对液膜方程进行讨论之前，作如下假设：（1）液膜厚度相比燃烧室曲率半径要足够小，以致沿燃烧室周向方向液膜的各项性质不变。（2）液膜流动处于层流状态并且速度沿液膜厚度呈线性变化。（3）液膜足够薄，以致液膜上的重力作用可以忽略不计。由文献证明方法可得在液体火箭发动机中，当液膜厚度小于100μm时，如上假定是合理的^[11]。

只考虑气相剪切力、牛顿内摩擦力、当前时刻撞击到液膜上液滴的冲击力与周围液膜压力梯度的影响。由于液体火箭发动机燃烧室中无旋转部件，故令壁面速度为0，动量方程为

$ \begin{array}{*{20}{l}} {{\rho _{\rm{l}}}{h_{\rm{l}}}\left\{ {\frac{{\partial {\mathit{\boldsymbol{u}}_{\rm{f}}}}}{{\partial t}} + \left[ {{\mathit{\boldsymbol{u}}_{\rm{f}}}\cdot{\nabla _{\rm{s}}}} \right]{\mathit{\boldsymbol{u}}_{\rm{f}}}} \right\} + {h_{\rm{l}}}{\nabla _{\rm{s}}}{p_{\rm{f}}} = }\\ {{\tau _{\rm{W}}}\mathit{\boldsymbol{t}} - {\mu _{\rm{l}}}\left( {\bar T} \right)\frac{{{\mathit{\boldsymbol{u}}_{\rm{f}}}}}{{{h_{\rm{l}}}/2}} + {{\dot P}_{{\rm{imp}}}} - \left( {{{\dot P}_{{\rm{imp}}}}\cdot\mathit{\boldsymbol{n}}} \right)\mathit{\boldsymbol{n}} - {{\dot M}_{{\rm{imp}}}}\cdot{\mathit{\boldsymbol{u}}_{\rm{f}}} + \delta {p_{\rm{f}}}\mathit{\boldsymbol{n}}} \end{array} $

(1)

离散差分，得

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{u}}_{\rm{p}}^{n + 1} = \left\{ {\frac{{{\rho _{\rm{l}}}{h_{\rm{l}}}}}{{\Delta t}}\mathit{\boldsymbol{u}}_{\rm{p}}^n + {\tau _{{\rm{w}}, {\rm{a}}}}{\mathit{\boldsymbol{t}}_{\rm{a}}} + {\mathit{\boldsymbol{S}}_{\rm{a}}} - \left[ {\frac{{{\rho _{\rm{l}}}{h_{\rm{l}}}}}{{\Delta t}}\mathit{\boldsymbol{u}}_{\rm{p}}^n\cdot{\mathit{\boldsymbol{n}}_{\rm{a}}} + {\mathit{\boldsymbol{S}}_{\rm{a}}}\cdot{\mathit{\boldsymbol{n}}_{\rm{a}}}} \right]{\mathit{\boldsymbol{n}}_{\rm{a}}}} \right\}/}\\ {\left[ {\frac{{{\rho _{\rm{l}}}{h_{\rm{l}}}}}{{\Delta t}} + \frac{{2{\mu _{\rm{l}}}}}{{{h_{\rm{l}}}}} + {{\dot M}_{{\rm{imp}}}}} \right]} \end{array} $

(2)

式中${\rho _{\rm{l}}}$代表液膜密度；${h_{\rm{l}}}$代表液膜厚度；${\mu _{\rm{l}}}$代表液体粘度; 其它各项表征的物理意义及求解公式参见文献[14]。

虽然Lagrange液膜动量方程求解的是液膜液滴的速度，但由于液膜厚度是当前网格中所有液膜液滴共同累积而成的，故其它液膜液滴对所计算液膜液滴的影响隐式反映到液膜厚度上，进而反映在气相剪切力和内摩擦力上。

液膜能量方程采用双线性温度模型，即把液膜内温度分为上下两段，如图 2所示。由界面守恒条件可知，气相场传给液膜的热量一部分用于液膜的蒸发，另一部分传入液膜内部。由于在液膜内部考虑双线性温度模型，则液膜能量方程为

$ \begin{array}{*{20}{l}} {{\rho _{\rm{l}}}{h_{\rm{l}}}{C_{{\rm{vl}}}}\left\{ {\frac{{\partial {{\bar T}_{\rm{l}}}}}{{\partial t}} + \left[ {\left( {{\mathit{\boldsymbol{{\bar u}}}_{\rm{l}}} - {\mathit{\boldsymbol{u}}_{\rm{w}}}} \right)\cdot{\nabla _{\rm{s}}}} \right]{{\bar T}_{\rm{l}}}} \right\} = }\\ {{\lambda _{\rm{l}}}\left( {{{\bar T}_{\rm{l}}}} \right)\left[ {\frac{{{T_{\rm{s}}} - {{\bar T}_{\rm{l}}}}}{{{h_{\rm{l}}}/2}} - \frac{{{{\bar T}_{\rm{l}}} - {T_{\rm{w}}}}}{{{h_{\rm{l}}}/2}}} \right] + {{\mathop Q\limits^\cdot }_{{\rm{imp}}}} - {I_{\rm{l}}}\left( {{{\bar T}_{\rm{l}}}} \right){{\dot M}_{{\rm{imp}}}}} \end{array} $

(3)

后两项可以化简消去，整理得

$ {\rho _{\rm{l}}}{C_{{\rm{vl}}}}{h_{\rm{l}}}\frac{{T_{\rm{p}}^{n + 1} - T_{\rm{p}}^n}}{{{\rm{\Delta }}t}} = \frac{{2{\lambda _{\rm{l}}}}}{{{h_{\rm{l}}}}}\left( {{T_{\rm{s}}} + {T_{\rm{w}}} - 2T_{\rm{p}}^n} \right) $

(4)

式中${C_{{\rm{vl}}}}$代表液体比热；${\lambda _{\rm{l}}}$代表液膜导热系数；${T_{\rm{s}}}$与$T_{\rm{w}}$分别代表液膜表面温度和壁面温度，其余各项物理意义见上式。

界面守恒方程为

$ {Q_{{\rm{gas}}}} = {\dot M_{{\rm{vap}}}}{H_{\rm{l}}} + {\lambda _{\rm{l}}}\frac{{{T_{\rm{s}}} - T_{\rm{p}}^{n + 1}}}{{{h_{\rm{l}}}/2}} $

(5)

式中${Q_{{\rm{gas}}}}$是气相场传给液膜的能量；${H_{\rm{l}}}$是蒸发潜热；${\dot M_{{\rm{vap}}}}$是液膜蒸发速率，计算公式在下文说明。

由于液膜采用的Lagrange方法计算，气相场采用Euler方法，如果液膜很薄的情况下，液膜对气相场计算的影响是可以忽略不计的。在统计发动机内流场向壁面传递热流的计算中，通常采用较薄的网格作为近壁处第一层网格从而令网格中气相温度处于线性层内，继而能准确地计算出向壁面传递的热流。无论在液膜的蒸发升温计算还是与结构的耦合传热计算，气相场传给液膜的热流都至关重要，假设液膜对气相场计算的影响是忽略不计的，则采用傅里叶导热定律计算气相场传给液膜的热流，即

$ {Q_{{\rm{gas}}}} = {\lambda _{{\rm{gas}}}}\frac{{{T_{{\rm{gas}}}} - {T_{{\rm{l}}, {\rm{s}}}}}}{{{h_{\rm{l}}}/2}} $

(6)

式中${\lambda _{{\rm{gas}}}}$为液膜表面所在网格的混合气体热导率；${T_{{\rm{gas}}}}$是近壁处第一层网格的气体温度；$h$是近壁处第一层网格的厚度。在此方程中，把液膜看成零厚度，主要是基于用Lagrange方法计算的液膜对用Euler方法计算的气相场影响不大，且近壁处第一层网格气相温度处于线性层内，从而能得出更准确的热流的考虑。

2.2.2 蒸发速率的求解

在没有到达平衡温度之前液膜的蒸发速率为

$ {\dot M_{{\rm{vap}}}} = {H_{\rm{Y}}}{\rm{ln}}\left( {\frac{{1 - {Y_{\rm{v}}}}}{{1 - {Y_{{\rm{vs}}}}}}} \right) $

(7)

式中${Y_{\rm{v}}}$是液膜表面所在气相场的燃料蒸汽质量分数；${Y_{{\rm{vs}}}}$是在液膜表面温度下的蒸汽平衡质量分数；$H_Y$是蒸发因子，具体计算公式参见文献[15]。

$ {Y_{{\rm{vs}}}} = \frac{1}{{1 + \left( {\frac{p}{{{p_{\rm{l}}}}} - 1} \right)\frac{{{M_\infty }}}{{{M_{\rm{l}}}}}}} $

(8)

式中$p$为液膜表面气体的压力；${p_{\rm{l}}}$是在液膜表面温度下的饱和蒸汽压，可采用Riedel蒸气压方程计算；${M_\infty }$是离液膜表面混合气体摩尔质量；${M_{\rm{l}}}$是液膜蒸汽的摩尔质量。

仿照液滴的蒸发模型^[19]，定义传热数${B_{\rm{T}}}$和传质数${B_{\rm{D}}}$

$ {B_{\rm{T}}} = \frac{{{C_{p, {\rm{g}}}}\left( {{T_{{\rm{gas}}}} - {T_{\rm{s}}}} \right)}}{{{H_{\rm{l}}} + {C_{{\rm{vl}}}}\left( {{T_{\rm{s}}} - {T_{\rm{p}}}} \right)}} $

(9)

$ {B_{\rm{D}}} = \frac{{{Y_{\rm{v}}} - {Y_{{\rm{vs}}}}}}{{{Y_{{\rm{vs}}}} - 1}} $

(10)

式中${C_{p, {\rm{g}}}}$是气体定压比热，其余物理量见上述。

随着液膜表面温度${T_{\rm{s}}}$升高，传热数${B_{\rm{T}}}$减小，传质数${B_{\rm{D}}}$增大，当表面温度升高到某一值使得${B_{\rm{D}}} = Le\cdot{B_{\rm{T}}}$时，则令此温度为液膜平衡温度，液膜表面温度限制在此温度下。$Le$为路易斯数，此时，进入液膜的热量刚好用于支付液体的蒸发潜热以及将蒸汽加热到环境气体的温度，蒸发速率为

$ {\dot M_{{\rm{vap}}}} = {Q_{{\rm{gas}}}}/\left( {{H_{\rm{l}}} + {C_{p, {\rm{g}}}}\left( {{T_{{\rm{gas}}}} - {T_{\rm{s}}}} \right)} \right) $

(11)

2.2.3 气液耦合作用

气相场中液滴对气相场的影响是通过向气相场加入源项体现出来的，由于液膜也采用Lagrange方法计算，故其对气相场的影响也可以仿照液滴的方法。则以下各式分别代表质量源项、动量源项和能量源项^[20]。液膜计算流程图如图 3所示。

$ {{S}_{\text{p}}}=\overset{{}}{\mathop{\mathop{\sum }^{}}}\,\left[ {{\left( {{m}_{\text{d}}} \right)}_{\text{in}}}-{{\left( {{m}_{\text{d}}} \right)}_{\text{out}}} \right]/\text{ }\!\!\Delta\!\!\text{ }V $

(12)

$ {{S}_{\text{p}}}=\overset{{}}{\mathop{\mathop{\sum }^{}}}\,\left[ {{\left( {{m}_{\text{d}}}{{\mathit{\boldsymbol{u}}}_{\text{d}}} \right)}_{\text{in}}}-{{\left( {{m}_{\text{d}}}{{\mathit{\boldsymbol{u}}}_{\text{d}}} \right)}_{\text{out}}} \right]/\text{ }\!\!\Delta\!\!\text{ }V $

(13)

$ {{S}_{\text{p}}}=\overset{{}}{\mathop{\mathop{\sum }^{}}}\,\left[ {{m}_{\text{d}}}{{\left( {{h}_{\text{d}}}+\frac{\mathit{\boldsymbol{u}}_{\text{d}}^{2}}{2} \right)}_{\text{in}}}-{{m}_{\text{d}}}{{\left( {{h}_{\text{d}}}+\frac{\mathit{\boldsymbol{u}}_{\text{d}}^{2}}{2} \right)}_{\text{out}}} \right]/\text{ }\!\!\Delta\!\!\text{ }V $

(14)

2.3 计算方法

采用分区求解的思路，各个计算区域采用独立的程序，只在耦合界面处进行数据传递。共分为四大模块：

（1）采用Euler-Lagrange-Lagrange法计算流场-离散相液滴-液膜液滴的燃烧流动模块。

（2）采用有限体积法计算强辐射能力气体（H₂O与CO₂）的燃气辐射模块。

（3）采用有限元法计算导热微分方程的结构热传导模块^[21]。

（4）采用斯忒藩-玻尔兹曼定律计算的辐射冷却模块。

四个模块在界面处传递数据的框图如图 4所示，流场在有液膜处采用液膜温度边界条件，在无液膜处采用壁面温度边界条件；流场计算收敛后把压力、温度和气体的摩尔分数传给燃气辐射模块；燃气辐射模块计算完成后，把辐射热流密度与流场的对流热流密度相加，作为结构热侧的热流密度边界；辐射冷却模块利用上一时间步结构冷侧的温度分布得出热流密度之后，作为结构冷侧的热流密度边界；结构两侧默认采用热流密度边界，当热侧某处有液膜并且结构温度高于液膜平衡温度时，结构在该处采用温度边界条件；当结构热侧吸热量和冷侧放热量相等时，耦合传热迭代计算收敛。

3 算例验证 3.1 仿真模型与计算条件

拟采用490N液体火箭发动机作为仿真模型，采用四氧化二氮和一甲基肼作为推进剂，推进剂物性、发动机结构与喷注方式详见文献[20]，燃烧室长度为78mm，以发动机喉部直径d作为特征长度，推力室工作参数由表 1所示。由于计算量过大，故采用1/8网格，图 5展示了流场和结构的网格划分图，整个计算域总的网格数量为31.2万个单元。其中，流场计算域网格单元个数为21.6万，固体结构域网格单元数是9.6万。对流场的头部、喉部和近壁处进行了加密处理，保证流场近壁处第一层网格较细以使网格中气相温度处于线性层内。为了实现可靠的热防护效果和形成较为均匀的液膜，冷却液膜喷注单元设计为与喷注面板呈22°和43°的两组喷射角度。流场计算的边界条件主要有：等温固壁边界、等温液膜边界、压力出口边界、周期边界、完全匹配对接边界；结构计算的边界条件主要有：内壁热流边界、内壁温度边界与外壁辐射热流边界；辐射计算的边界条件主要是外壁等温边界条件。

流场计算采用三维可压缩N-S方程，对控制方程中对流项的离散采用HLLC格式，通过采用MUSCL近似方法得到二阶精度格式。粘性项用中心差分离散，关于时间的离散本文采用LU-SGS隐式离散方法^[22]。液滴的粒径分布采用R-R分布，液滴运动模型详见文献[18]，蒸发模型采用自燃推进剂高压蒸发模型。燃烧模型采用利用Arrehenius公式计算化学反应源项的有限速率模型，一甲基肼和四氧化二氮的化学反应动力学采用5步反应机理，由表 2所示。

气体辐射模块采用有限体积法^[23]，把整个计算域进行空间离散和角度离散，保证控制体的每个立体角内辐射能量守恒。通过立体角内能量守恒的原理，计算气体辐射能量。

结构采用了对不规则计算区域适应性更强的有限元方法离散求解三维傅里叶导热微分方程^[24]。

辐射冷却模块把结构冷侧各单元面分别用斯忒藩-玻尔兹曼求解热流密度，再把其返还给结构。

流场边界网格单元与结构边界网格单元并不是完全贴合的，由于是三维到三维之间的数据传递，故采用几何意义明确，在数据足够密集的情况下插值精度足够高的三线性插值方法进行界面处插值。

3.2 计算结果及分析

流场和结构都是采用1/8发动机网格进行计算，但在统计热流的时候，按照发动机整体结构进行考虑。本文采用稳态计算方法，由能量守恒得知，当推力室内壁面的吸热量与推力室外壁面的放热量一致并且吸（放）热量保持不变时，耦合传热迭代收敛。图 5显示的是耦合迭代过程中推力室内壁面吸热量和外壁面放热量的变化情况，由于吸热量与放热量逐渐一致而且沿水平方向波动，可认为耦合迭代收敛。可见，发动机流场传给结构的总热流在11.8kW左右。

由图 6（a）可见耦合界面处温度连续，所选取的插值方法精度较高。仿真结果选取三个测量点，点a对应喷射夹角22°，点b对应1/8计算域中间，点c对应喷射夹角43°。可见，传热耦合计算得到的液膜覆盖段温度为425K，是液膜在此环境条件下的平衡温度，头部结构温度平均为530K左右，比实验实测值433K高100K左右。发动机喉部结构平均温度为1320K左右，比实验实测值1623K低300K左右。扩展段仿真温度相比实验结果偏高。整体来说，液膜冷却效果较好，喉部温度计算较低，扩张段温度计算稍高。

本文在射流撞壁行为准则中只考虑射流撞击壁面后铺展成膜的情况，液膜干涸后会形成气膜向下游移动。如图 7（a）所示，壁面处的CH₃NHNH₂极易分解产生CH₄，过多CH₄气膜积聚在壁面附近未参与燃烧致使燃烧不够充分，而过量的液膜和气膜又起到较强冷却作用，这两种原因导致喉部温度过低。气膜不光阻碍燃烧室的气流向壁面传热，还阻碍扩张段的气流冷却壁面，这就造成了扩张段结构温度下降较慢。

图 7（b）显示的是在Z=0处内流场温度分布云图，最高温度为3000K左右。燃烧室压力为0.75MPa，相比于实验值0.85MPa较低，是由于CH₃NHNH₂分解产生的部分CH₄积聚在壁面导致燃烧不充分。

图 8展示了结构热侧和结构切片温度分布图，可见液膜冷却效果很明显，图 9展示液膜分布形态及液膜液滴分布形态，液膜厚度是通过统计每个网格中液膜液滴的体积之后得出的，最大液膜厚度为60μm，小于100μm，薄膜假设正确。可见与面板夹角较小的喷雾形成的液膜离头部较近，形成的液膜更集中，生存时间更长，因此液膜延展的距离更长，约为30mm，到达收敛段，阻隔燃气热流效果更好。与面板夹角较大的喷注方式形成的液膜周向分布较广，有效阻止了结构热流流向头部。

4 结论

通过本文研究，得到如下结论：

（1）Euler-Lagrange两相燃烧流动的框架下，Lagrange液膜模型计算得到的最大液膜厚度为60μm，最大液膜长度为30mm，能到达收敛段，头部温度为433K，液膜覆盖段壁温为425K，收敛段温度为700~1000K，喉部温度为1320K。喉部温度计算稍低与射流撞壁行为准则的局限性有关，但所提出的考虑液膜条件下的耦合传热方法是正确的。

（2）与面板的小喷注夹角形成的液膜离头部较近且延展距离较长，有效阻隔了燃气传给壁面热流，大喷注夹角形成的液膜周向扩展范围更大，有效阻止了结构热流向头部传递。

本文令射流撞击壁面后全部形成液膜，没有考虑反弹和飞溅等行为，由此导致积聚气膜过多以致燃烧不够充分，致使室压偏低，过多气膜既阻碍了收敛段气流加热室壁又阻碍扩张段气流冷却室壁，不充分的燃烧和低温气膜隔热作用导致喉部温度偏低。构建适用液体火箭发动机工作条件下的射流撞壁行为准则，将是今后工作的重点。