0 引　言

鱼雷燃烧室是鱼雷热动力装置中最重要的部件之一^[1]，HAP三组元推进剂在燃烧室内燃烧，产生高温高压的燃气作为工质推动发动机做功^[2]。为保证燃烧室在工作过程中稳定可靠，燃烧室壳体需进行合理设计^[3]。现有结构下的燃烧室，在工作后其内壳体上有明显的烧蚀环槽。

经分析是因公差燃烧室内壳体与内衬在接触位置处存在缝隙，燃气经缝隙泄漏造成冲刷现象，形成烧蚀环槽。为解决这一问题，结构改进的一种方案为在内衬上增开一对压力平衡孔，以解决壳体的烧蚀现象。而压力平衡孔的位置及大小将直接影响结构改进方案的效果。

热流固耦合数值仿真的方法广泛应用于鱼雷、航空等各领域的燃烧室的仿真计算工作当中^[4-10]。本文针对具有不同压力平衡孔位置及大小的鱼雷燃烧室结构为研究对象，利用Fluent软件开展燃烧室高工况工作条件下的热流固耦合数值仿真计算，分析燃烧室内衬压力平衡孔的位置及大小对燃气域流场及燃烧室内壳体应力场的影响规律。

1 燃烧室工作过程中物化过程的数学模型

鱼雷燃烧室工作时，涉及到HAP三组元推进剂破碎雾化、液滴颗粒及内流场气相的两相流动、液滴颗粒的蒸发燃烧、燃烧室壳体及冷却水与燃气之间的耦合换热这4个过程。

1.1 液滴破碎控制方程

液体射流在空气中做相对运动，空气动力促使射流破碎；液体表面张力使射流不破碎。在2个力的作用下液膜破碎成小液滴，如图1所示。

由射流扰动机理，对三组元推进剂的破碎状态定义破碎系数 ${R}$ ，其数学计算式为：

$R = \frac{\lambda }{{{d_0}}} = {\text{π}} \sqrt 2 {\left(1 + \frac{{3\mu }}{{\sqrt {\rho \sigma {d_0}} }}\right)^{0.5}},$

(1)

${d^*} = 1.89{d_0} 。$

(2)

式中： $\lambda$ 为射流扰动波长； ${d_0}$ 为射流起始直径； $\sigma$ 为液体表面张力； ${d^*}$ 为破碎后的液滴的标称直径。

取R-R分布为液滴雾化后的直径分布依据，按中间质量分布准则作为直径的分布准则，可知液滴直径 $d$ 的分布表达式为：

$R = 1 - \exp \left\{ - 0.693\left[\frac{d}{{{d^*}{{(0.693)}^i}}}\right]\right\} ,$

(3)

式中， $i$ 为液滴直径分布指数。

1.2 液滴运动及气相两相流动控制方程

破碎后的液滴，在运动过程受到燃气粘性的摩擦阻力，无体积力。离散相燃料液滴的运动控制方程为：

$\frac{{{\rm{d}}{u_p}}}{{{\rm{d}}t}} = {F_D}(u - {u_p}) + \frac{{{g_x}({\rho _p} - \rho )}}{{{\rho _p}}},$

(4)

${F_D} = \frac{{18\mu }}{{{\rho _P}{d^2}}}\frac{{{C_D}Re}}{{24}} 。$

(5)

式中： $u$ 为流体相速度； ${u_p}$ 为颗粒速度； $\ \rho$ 为流体密度； $\ {\rho _p}$ 为颗粒密度； $d$ 为颗粒直径； $Re$ 为颗粒雷诺数； ${C_D}$ 为阻力系数。

对于气相部分，连续相流体运动控制方程为：

1）连续性方程

$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho u)={S}_{m} 。$

(6)

式中： ${S_m}$ 为质量源相，数值上等于注入的燃料在单位时间蒸发的质量流量。

2）动量方程

广义上连续相流体运动的控制方程为^[9]：

$\frac{\partial }{{\partial t}}(\rho {\boldsymbol{u}}) = - \nabla P + \mu {\nabla ^2} + \rho {\boldsymbol{g}} + {\boldsymbol{F}} + {S_M} 。$

(7)

式中： $P$ 为静压； $\mu$ 为气相动力粘性； ${\boldsymbol{F}}$ 为外部体积力； ${S_M}$ 为动量源相； $\nabla$ 为拉普拉斯算子。

对于燃烧室气相，因重力因素忽略，无外部体积力，故气相简化后的动量方程为：

$\frac{\partial }{{\partial t}}(\rho {\boldsymbol{u}}) = - \nabla P + \mu {\nabla ^2} + {S_M} 。$

(8)

3）湍流模型

燃烧室燃烧过程中燃气处于高湍流状态，选用 $k - \varepsilon$ 湍流模型，其控制方程为：

$\frac{\partial }{{\partial t}}(\rho k) + \frac{{\rm{d}}}{{{\rm{d}}x}}(\rho uk) = \frac{{\rm{d}}}{{{\rm{d}}x}}\left[\left(\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}\right)\frac{{{\rm{d}}k}}{{{\rm{d}}x}}\right] + {G_k} + {G_b} - \rho \varepsilon - {Y_M} 。$

(9)

$\frac{\partial }{{\partial t}}(\rho \varepsilon ) + \frac{{\rm{d}}}{{{\rm{d}}x}}(\rho u\varepsilon ) = \frac{{\rm{d}}}{{{\rm{d}}x}}\left[\left(\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}\right)\frac{{{\rm{d}}k}}{{{\rm{d}}x}}\right] + C\frac{\varepsilon }{k}({G_k} + {G_b})。$

(10)

式中： $\ \mu$ 为粘性系数； $\ {\mu _t}$ 为湍流粘性系数； ${G_k}$ 为速度梯度湍动能； ${G_b}$ 为浮力湍动能，对不可压流体该值为0； ${Y_M}$ 为湍流脉 动湍动能； $C$ 为修正系数常量； ${\sigma _k}$ ， ${\sigma _\varepsilon }$ 为湍流普朗特数。

1.3 液滴蒸发燃烧控制方程

1）液滴蒸发控制方程

液滴未达到蒸发温度的温度控制方程^[10]为：

${m_p}{c_p}\frac{{{\rm{d}}{T_p}}}{{{\rm{d}}t}} = h{A_p}({T_\infty } - {T_p}) + {\varepsilon _p}{A_p}\sigma (\theta _R^4 - T_p^4)。$

(11)

达到蒸发温度后的液滴质量控制方程^[10]为：

$0 = h{A_p}({T_\infty } - {T_p}) + {\varepsilon _p}{A_p}\sigma (\theta _R^4 - T_p^4) + \frac{{{\rm{d}}{m_p}}}{{{\rm{d}}t}}{h_{fg}} 。$

(12)

式中： ${m_p}$ ， ${c_p}$ ， ${T_p}$ ， ${\varepsilon _p}$ ， ${A_p}$ 分别为颗粒的质量、比热、温度、发射率和表面积； $h$ 为对流传热系数； ${T_\infty }$ 为气相温度； $\sigma$ 为波尔兹曼常数； ${\theta _R}$ 为辐射温度； ${h_{fg}}$ 为潜热。

2）液滴燃烧控制方程

HAP三组元推进剂完全燃烧的化学反应方程式为：

$\begin{split}&\rm C_{4.8}H_{73.41}O_{45.04}N_{4.16}Cl_{2.53} \to \\ &\rm 4.8CO_2 + 2.08N_2 + 35.44H_2O + 2.53HCl 。\end{split}$

反应释放的能量由反应物与生成物的总焓差决定。

Eddy-dissipation涡耗散模型的反应速率 ${{R}}_{{i},{k}}$ 控制方程为^[10]：

${R_{1,i,k}} = \nu _{i,k}'{M_{w,i}}AB\rho \frac{\varepsilon }{\kappa }\frac{{\displaystyle\sum {_p{m_p}} }}{{\displaystyle\sum\nolimits_j^N {\nu _{j,k}^{''}{M_{w,j}}} }} ,$

(13)

${R_{2,i,k}} = \nu _{i,k}'{M_{w,i}}A\rho \frac{\varepsilon }{\kappa }\mathop {\min }\limits_R \left(\frac{{{m_R}}}{{\nu _{R,k}'{M_{w,R}}}}\right),$

(14)

${R_{i,k}} = \min ({R_{1,i,k}},{R_{2,i,k}})。$

(15)

式中： ${m_p}$ 为所有产物组分和某一特定反应物的质量比重； ${M_{w,j}}$ 为物质的分子量； $\nu _{i,k}'$ ， $\nu _{j,k}^{''}$ 分别为反应物和生成物的化学当量系数； $A$ 和 $B$ 为经验常数，分别取4和0.5。

1.4 耦合换热过程能量控制方程

流体微元体微分形式的能量方程为：

$\frac{\partial \rho e}{\partial t}+h(\rho ,\mu ,\lambda ,{c}_{p})\cdot {\displaystyle \sum \nabla T=\nabla \cdot (\lambda \nabla T)}-P\nabla \cdot u+\varPhi +Q ,$

(16)

$\varPhi = \frac{1}{2}\mu {\left(\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}}\right)^2} + \left({\mu '} - \frac{2}{3}\mu \right){\left(\frac{{\partial {u_i}}}{{\partial {x_i}}}\right)^2}。$

(17)

式中： $\dfrac{{\partial \rho e}}{{\partial t}}$ 为单位时间内流体微元体的能量变化；对流换热系数 $h$ 为流体的密度 $\rho$ 、动力粘度 $\mu$ 、导热系数 $\lambda$ 、定压比热容 ${c_p}$ 的状态函数； $\nabla \cdot (\lambda \nabla T)$ 为微元体通过控制体边界经热传导导入控制体的能量； $P\nabla \cdot u$ 为微元体边界逆压力做功消耗的能量； $\varPhi$ 为流体微元体粘性耗散产生的能量； $Q$ 为单位时间里外界产生的热量。

2 数值计算几何模型及网格划分 2.1 几何模型预处理及域提取

对燃烧室的几何结构进行简化，以降低网格划分难度、降低网格数量、提高网格质量、帮助模型更快更好的收敛。

利用提取功能处理简化完成的燃烧室固体域，得到研究对象的全体计算域，如图2所示。

保温层域的缝隙厚度由燃烧室内壳体与内衬的最大加工公差决定，缝隙厚度为0.166 mm，上下两处缝隙位置的几何结构局部放大图如图3所示。

2.2 计算域网格划分

采用混合网格的划分方式完成数值计算模型的网格划分，其中使用ICEM软件的block功能进行结构化网格的划分，使用ICEM软件中的八叉树算法进行非结构化网格的划分。在Fluent中使用make polyhedral功能将四面体网格转为多面体网格，最终划分的网格数量在千万级左右，网格划分结果如图4所示。

2.3 网格无关性验证

验证网格尺寸不会对数值计算结果造成显著影响，开展网格无关性验证。进行网格无关性验证的网格划分方案如表1所示。

本文针对燃烧室燃气域的流场进行研究，故选取燃烧室燃气域出口段的温度变化情况作为衡量网格无关性的指标，网格无关性验证的结果如图5所示。

结合表1的网格划分方案及图5的网格计算结果，可知网格划分方案2为最优的网格划分方式，可作为之后数值计算中的网格划分方案。

3 边界条件及压力平衡孔几何参数设置 3.1 边界条件设置

针对燃烧室热流固耦合数值仿真计算，以燃烧室在高工况下的燃料、冷却水、内压等工况参数作为数值仿真计算的边界参数。具体的边界条件参数设置如表2所示。

3.2 压力平衡孔几何参数

设置9组不同压力平衡孔位置及尺寸的数值仿真计算case，用于比较压力平衡孔位置及大小对燃烧室燃气域的流场及壳体的应力场的影响，具体结构参数设置如表3所示。

4 数值计算结果及分析 4.1 压力平衡孔位置对燃烧室流场及壳体应力场的影响

1）上方缝隙处

不同压力平衡孔位置下，燃气域上方缝隙处的压力、速度、热通量变化图如图6～图8所示。

由图6可知，压力平衡孔与内衬上表面距离越小，缝隙两端的压差越小。经分析原因在于燃气腔内的燃气越接近内衬上表面，燃气的压力值与内衬上表面处的压力值越接近，开通压力平衡孔后压力联通，使得缝隙两端的压差减小。

由图7可知，压力平衡孔位置的改变不对缝隙处燃气的速度变化规律造成影响，但对速度的大小会产生显著影响，且压力平衡孔与内衬上表面越近，缝隙内的燃气流速越低，这是由于缝隙两端的压差减小导致的。平衡孔位置的改变对缝隙内燃气流速的最大影响值约为3 m/s。

缝隙内燃气流速的降低，会导致燃气对流换热量的减小，所以压力平衡孔与内衬上表面的距离越小，缝隙处的热通量大小越小，图8可以印证这一分析。

2）下方缝隙处

不同压力平衡孔位置下，燃气域下方缝隙处的压力、速度、热通量变化图如图9～图11所示。

可知，改变压力平衡孔的位置对下方缝隙处的燃气的压强、流速、缝隙处的热通量等参数均不产生较大影响。分析认为原因是燃气在外腔内流动，当流至下方缝隙的上方时燃气的流动状态已经趋于稳定，故压力平衡孔对下方缝隙处的燃气流动已不会造成显著影响。

3）壳体应力场

图12为不同压力平衡孔位置燃烧室壳体的应力场云图。可知，随着压力平衡孔与内衬上表面的距离的减小，壳体缝隙处的最大应力值减小。结合图6～图8，可知原因在于缝隙两端压差降低，降低了燃气流速从而减少了缝隙处的热通量。降低了壳体在缝隙处的温度，从而减小了热应力的大小。

4.2 压力平衡孔大小对燃烧室流场及壳体应力场的影响

1）上方缝隙处

不同压力平衡孔直径下，燃气域上方缝隙处的压力、速度、热通量变化图如图13～图15所示。

由图13可知，随着压力平衡孔的直径增大，缝隙两端的压差先变小后增大。分析认为原因在于随着孔径的增大，燃气域外腔内的燃气压力会越来越接近内腔，使得缝隙两端压差减小，但孔径过大会使燃气在孔的末端发展出湍流，在湍流中心处形成低压区，反而使得缝隙两端压差增大。

由图14可知，压力平衡孔直径的改变不对缝隙处燃气的速度变化规律造成影响，但对速度的大小会产生显著影响。随着压力平衡孔直径的增大，缝隙内的燃气流速先降低后增高，这是由于缝隙两端的压差先变小后变大所导致的。且在孔径大小为4 mm时，缝隙内的燃气流速最低为8.5 m/s。随着孔径的改变，缝隙内燃气流速的最大差值约为2 m/s。

由图15可知，上方缝隙处燃气的对流换热量与燃气的流速具有相同的变化规律，即先变小后变大。

2）下方缝隙处

不同压力平衡孔直径下，燃气域下方缝隙处的压力、速度、热通量变化如图16～图18所示。

可以看出，改变压力平衡孔的直径对燃烧域下方缝隙处的缝隙两端压强、缝隙内的燃气流速、缝隙处的热通量等参数均不产生较大影响，经分析认为原因在于燃气在外腔内流动，当流至下方缝隙上方时燃气的流动状态已经趋于稳定，故压力平衡孔对下方缝隙处的燃气流动不会造成影响。

3）壳体应力场

图19为不同压力平衡孔直径燃烧室壳体的应力场云图。可知，随着压力平衡孔直径的增大，壳体缝隙处的最大应力值先减小后增大，且当孔径为4 mm时，壳体所受的最大应力值最小。结合图13～图15，可知原因在于缝隙两端压差降低，降低了燃气流速从而减少了缝隙处的热通量。降低了壳体在缝隙处的温度，从而减小了热应力的大小。

5 结　语

本文针对压力平衡孔的位置及大小对鱼雷燃烧室的燃气域流场及壳体域应力场的影响规律进行数值计算，经分析可以得出如下结论：

1）对于燃气域上方缝隙处，压力平衡孔的位置距离内衬上表面越近，燃气域缝隙两端的压差越小、燃气的流速越低、缝隙处的热通量越小；对于下方缝隙处，压力平衡孔的位置变化对燃气域无明显影响；对于壳体的最大应力值，压力平衡孔距离内衬上表面越近，最大应力值越小。

2）对于燃气域上方缝隙处，随着压力平衡孔的直径变大，燃气域缝隙两端的压差先变小后变大、燃气的流速先降低后增高、缝隙处的热通量先减小后增大；对于下方缝隙处，压力平衡孔的直径变化对燃气域无明显影响；对于壳体的最大应力值，随着压力平衡孔直径的增大，最大应力值先减小后增大，且当直径为4 mm时最大应力值最小。