固体火箭发动机喷管两相流动下的热固耦合研究

李 煜，程广益，李宗阳，窦怡彬，蒋君庭

（上海机电工程研究所，上海 201109）

0 引 言

喷管是固体火箭发动机的基本部件之一。喷管部件在高温燃气的作用下温度急剧上升，在进行火箭发动机各部件结构设计时，需要考虑各部件受热状态下的实际强度以保证发动机稳定工作。此外，燃气传热会引起各部件内部温度分布不均，从而引入热应力载荷的作用。因此，要进行合理的火箭发动机喷管结构设计就必须考虑传热的影响。在两相流动时，颗粒在喷管内会不断撞击喷管内壁，对喷管的热防护、变形量以及应力的大小和分布造成影响［1］。因此，研究喷管在两相流动下的传热特性及热固耦合有重要意义。

本文考虑喷管结构（包括烧蚀层、喉衬、绝热层以及壳体）的传热特性、喷管外壁与大气之间的传热因素，分析喷管结构对内部温度场的影响；利用颗粒轨道模型［2］，研究不同颗粒直径D和不同负载比η时两相流动对喷管结构内部的温度场的影响；进行热固耦合计算，研究纯气相和两相流动时喷管的热变形分布，分析两相流动对热变形的影响。

1 数值计算方法

计算时进行以下简化：不考虑化学反应；不考虑颗粒相的燃烧、蒸发、破碎过程，且各个工况均为稳态计算；不考虑内壁面的烧蚀；不考虑颗粒接触传热；不考虑各层材料之间的接触热阻；燃气为理想气体。

1．1 计算模型

考虑到计算量较大，本文采用了二维轴对称模型。喷管结构如图1所示。各结构材料分别为：耐烧蚀层为石墨布；喉衬为碳碳复合材料；绝热层为碳酚醛；壳体为45号钢。各材料物性参数如表1所示［3－7］。

图1 喷管模型Fig．1 The model of nozzle

表1 材料物性参数Tab．1 The physical property parameters of the materials

1．2 计算方法

对于高速可压缩流场计算，通常使用密度基求解［8］，本文采用隐式Reo格式求解方法。目前对固体火箭发动机流场进行模拟时，湍流模型大多采用标准k－ε模型，可取得较高的精度。因此，本文主要采用标准k－ε模型来进行计算。

描述气固两相流动时，采用颗粒轨道模型，颗粒相的瞬时速度vp由式（1）积分确定。

式中：mp为颗粒相的质量；vp为颗粒相速度；t为时间；Fp为拖曳力，表达式为

式中：ρ为气相密度；dp为颗粒直径；CD为拖曳系数，v为气相速度。

颗粒的位置可通过式（3）计算。

式中：Xp是t＋Δt时刻颗粒的位置；为t时刻颗粒的位置，对式（1）和式（3）进行数值积分可以获得颗粒的轨迹。

1．3 边界条件

根据实际情况给定适当的边界条件。假定壁面无滑移，流体与喷管内壁为耦合传热。喷管入口设为压力入口，给定压力入口总压7 MPa，总温3 000 K；喷管出口为压力出口；出口和入口均给定湍流强度和水力直径。喷管外壁面、壳体尾部端面、绝热层尾部端面、烧蚀层尾部端面的对流换热系数根据大空间自然对流换热计算公式得到［9］，为10 W／（m2·K），发射率为0．8。壳体前部端面、绝热层前部端面、烧蚀层前部端面均为绝热面，相邻固体材料之间为内部壁面。外部自然来流温度为300 K。颗粒相采用三氧化二铝，密度为2 970 kg／m3，比热容为897 J／（kg·K），导热率为202．4 W／（m·K），计算选择的颗粒直径分别为1!m、10!m、20!m 和50!m，负载比分别为0．01、0．05、0．1、0．2 和0．3。

1．4 网格无关性验证

不同数量的网格可能会对计算结果产生一定的影响，因而需要进行网格无关性验证，确保网格对计算结果不会产生明显影响。整体网格及壁面边界层附近的网格如图2所示。本文选取的3种网格数量及相关的纯气相计算结果如表2所示。

图2模型网格Fig．2 The mesh of model

表2 网格数量及计算结果Tab．2 Grid number of meshes and the calculation results

经对比发现，3种网格对计算结果的影响甚微。为了节省计算资源，本文选择Mesh1作为计算网格。

2 流场计算结果与分析

2．1 纯气相流场分析

图3所示为流固交界面、耐烧蚀层和绝热层交界面、绝热层和壳体交界面以及壳体外壁面的温度分布。图4为文献［3］中流固交界面上的温度分布。对比图3和图4可以发现，由于文献中的燃面不断推移，因而在图中左侧有较长的平直段，同时文献中的喷管扩张段较短，因而温度下降较为明显，但在喷管的收缩段及喉部的温度分布和文献中的喷管内壁面温度分布一致，从而验证了本文计算结果的准确性。

由图3可以明显地看出：绝热层碳酚醛的导热率较小，对壳体起到了很好的隔热作用；烧蚀层的石墨布、喉衬的碳碳复合材料和壳体的钢导热率较大，因而在烧蚀层和壳体内部温度变化较小，交界面之间的温度变化也较小。

各界面温度的最高点出现在喷管喉部的上游附近，沿喷管轴线方向，温度先增大后减小。喷管内流体温度沿轴线方向逐渐减小，由于喉衬的导热率较大、传热快，大量热量通过交界面迅速地传递给烧蚀层，因而在上游出现最高温。下游的流体温度相对较低，温度呈现下降趋势。喷管尾部温度下降的梯度较大，这是由于尾部端面的对流换热及辐射换热导致了大量热量散失到外部自然来流。

图3 各交界面的温度分布Fig．3 Temperature distribution at each interface

图4 文献［3］中流固交接面的温度分布Fig．4 Temperature distribution at the fluid－solid interface in the literature［3］

2．2 不同颗粒直径对流场的影响

计算不同颗粒直径的两相流动时，选择的颗粒直径D分别为1!m、10!m、20!m和50!m，负载比取0．3。

图5给出了纯气相及不同直径颗粒下两相流场中喷管温度沿轴向的分布。图6为轴线气相的马赫数分布。对比文献［10］发现，本文的计算结果在轴向位置0．05 m之前的曲线分布趋势同文献计算结果保持一致。这是由于两相流动时颗粒相的惯性以及颗粒相与气相之间的黏性，轴线附近的流场受到颗粒相的传热作用，使得两相流动时轴线温度较纯气相情况下的轴线温度更高。但在轴向位置0．05 m之后，曲线的分布有较大差异。这是由于本文的计算模型和文献［10］所采用的喷管模型不相同，本文所选计算模型的扩张段更长，在扩张段尾部，直径较大颗粒的随流性较差，形成的无颗粒区域较大，大量粒子在靠近轴线附近运动，因而轴线附近的温度就较高。这也说明，在喷管喉部下游的一小段区域内，轴线马赫数和温度分布具有一定的规律，但是在下游一定区域之后，其趋势并不呈现出明显的规律。由此可知，不同颗粒直径下两相流场沿轴线的温度分布与模型扩张段有重要关系。

图5 不同颗粒直径时轴线温度分布Fig．5 Temperature distribution on the axis with different particle diameters

图6 不同颗粒直径时轴线气相马赫数分布Fig．6 Mach number distribution of axial gas phase with different particle diameters

图7给出了纯气相及不同直径颗粒下流固交界面、绝热层和烧蚀层交界面、壳体与绝热层交界面以及壳体外壁的温度分布。由图7可知，两相流场对喷管结构内部温度场的影响很大。与纯气相情况相比，两相流动下，喷管各交界面在扩张段的温度明显偏高，而在喷管收缩段温度无明显差异，并且颗粒直径越大，扩张段的温度越低。

图7 不同颗粒直径时喷管各交界面温度分布Fig．7 Temperature distribution at each interface of the nozzle with different particle diameters

2．3 不同负载比对流场的影响

针对不同负载比，计算两相流动时的流场参数，选取的负载比分别为0．01、0．05、0．1、0．2和0．3，颗粒直径为1!m。

图8和图9分别为不同负载比时沿轴线的温度和马赫数分布。

图8 不同负载比时沿轴线的温度分布Fig．8 Temperature distribution of the axis at different load ratios

图9 不同负载比时沿轴线气相的马赫数分布Fig．9 Mach number distribution of axial gas phase at different load ratios

图10所示分别为沿流固交界面、绝热层和烧蚀层交界面、壳体和绝热层交界面以及壳体外壁的温度分布。由图10可以看出，随着负载比的减小，各交界面的温度在喷管喉部及扩张段呈现下降的趋势，而在收缩段温度大小和变化趋势基本一致，说明负载比对喷管喉部及扩张段结构的温度有较大的影响。

图10 不同负载比时喷管各交界面的温度分布Fig．10 Temperature distribution at each interface of the nozzle at different load ratios

3 热固耦合计算结果与分析

通过对流场参数的计算获得流固耦合交界面上的温度分布后，可以进一步计算固体域的热变形。为了对比分析两相流动对喷管热变形的影响，本文主要计算纯气相、不同颗粒直径以及不同负载比时两相流动的喷管热固耦合，工况均与第2章所述工况相同。表3给出了喷管各结构材料的力学参数。

喷管的网格划分采用多区域法，并且在前部和尾部端面以及螺栓和壳体连接处加密。总网格数量为27．3万。图11给出了喷管的整体网格。

设置喷管外壁面及尾部端面为对流换热和辐射换热，对流换热系数和发射率以及外部温度的设置均与流场计算时采用的设置相同。设定螺栓的3个面为固定面，相邻结构之间采用绑定接触，并且采用对称接触。

图12所示为纯气相时喷管的温度分布，可以看出，温度的最高点出现在喷管喉部上游附近处，并且沿轴线方向先增大后减小，与流场计算结果的分布规律以及数值大小基本一致。

图11 喷管整体网格Fig．11 The grids of the nozzle

图12 喷管温度分布Fig．12 Temperature distribution of nozzle

3．1 纯气相热变形分析

图13所示为纯气相情况下喷管的变形量分布。喷管的最大变形出现在喷管尾部，变形量为7．836 6 mm，最小值在螺栓处。由于螺栓的固定作用，热膨胀产生的变形只能向自由端传播，变形量沿轴线不断增大，在尾部达到最大；同理，轴向变形的变形量分布与总变形量分布相同。径向变形的最大值也出现在喷管的尾部端面，变形量为2．056 6 mm；最小值出现在烧蚀层前部，变形量为0．352 7 mm。由于壳体前端面固定，壳体在该处的变形量为零，而烧蚀层在壳体内侧，其径向变形指向内侧，即指向圆心，因而变形量在数值上为负值。喷管外部的径向变形量的分布与喷管内壁温度的变化趋势一致，先增大后减小再增大，在喷管喉部上游附近出现高温。

图13 纯气相情况下喷管的变形量分布Fig．13 Displacement distribution of nozzle in pure gas phase

3．2 两相流动对喷管变形影响分析

图14为颗粒直径D＝1μm时喷管的各变形量分布。由图14可以看出，喷管变形量的分布趋势同纯气相情况下的分布趋势一致，但两相流动下的喷管的总变形量、轴向变形量以及径向变形量比纯气相情况下的量值要大，这是因为两相流动下的喷管结构内部温度相对较高。

图14 颗粒直径1μm情况下喷管的变形量分布Fig．14 Displacement distribution of nozzle with 1μm partical diameter

图15为不同颗粒直径时喷管各变形量的最大值。比较发现，两相流动时的总变形量、轴向变形量和径向变形量的最大值比纯气相情况下的要大，并且随着颗粒直径的增大，各变形量的最大值逐渐减小。表4为不同颗粒直径时各变形量最大值的变化率，相比于纯气相时的各变形量最大值，D＝1μm时的总变形量最大值增大了12．9%，D＝50μm时总变形量最大值增大了2%。随着颗粒直径的增大，各变形量最大值的变化率也不断减小。

图15 不同颗粒直径时各变形量最大值Fig．15 Maximum deformation at different particle diameters

表4 不同颗粒直径时各变形量最大值变化率Tab．4 The change rate of maximum deformation at different particle diameters

图16为不同负载比时喷管各变形量的最大值。两相流动情况下，喷管的各变形量最大值均大于纯气相情况下的值，并且随着负载比的增大，各变形量的最大值逐渐增大。表5为不同负载比时各变形量最大值的变化率，相比于纯气相时的各变形量最大值，η＝0．01时总变形量最大值增加了0．4%，η＝0．3时总变形量最大值增加了12．9%。随着负载比的增大，变化率不断增大。

图16 不同负载比时的变形量Fig．16 Maximum deformation at different load ratios

表5 不同负载比时各变形量最大值变化率Tab．5 The rate of change of maximum deformation at different load ratios

4 结束语

本文计算了纯气相条件下喷管的流场并利用颗粒轨道模型计算了两相流场；分析了喷管结构内部温度场的变化以及两相流动对喷管内部温度的影响；通过热固耦合计算，分析了变形量的分布以及两相流动对喷管热变形的影响，研究表明：

1）在考虑喷管结构传热性能时，流固交界面及喷管内部的温度分布与绝热壁面的温度分布有较大差异；流固交界面及喷管的各交界面沿轴向温度先增大后减小，在喉部上游附近达到最大值，在喷管尾部端面处温度会发生骤降；喷管的最低温度出现在绝热层尾部端面径向位置最大处，在相邻结构交界处温度会急剧变化；绝热层内的温度梯度很大，降温幅度很大，而烧蚀层和壳体内部的温度梯度很小，降温幅度很小；

2）不同颗粒直径时的两相流动对喷管轴线处的温度影响仅在喉部下游的一小段区域内呈现规律性的变化，在此段区域之外，颗粒直径的影响不遵循相应的规律；

3）在温度场的作用下，喷管总变形量和轴向变形量沿轴线方向不断增大，在尾部变形量最大，径向变形量沿轴向先增大后减小再增大，在喷管尾部达到最大，最小值出现在烧蚀层前部；

4）两相流动对喷管的热变形有着较大影响，喷管的总变形量、轴向变形量及径向变形量均增大。