拉瓦尔喷管流固耦合换热数值模拟

龚红兰，李 凌

(1.上海理工大学能源与动力工程学院，上海 200093；2.上海市动力工程多相流动与传热重点实验室，上海 200093)

1 引言

拉瓦尔喷管广泛应用于工程中，比如超音速喷气式飞机、火箭等飞行器中的喷管都属于拉瓦尔喷管[1]。在发动机工作时常伴随着高温高压燃气的流动，内部热环境恶劣，对热防护层材料性能要求较高，为此，许多文献对拉瓦尔喷管内的流动换热以及固体壁面导热进行了研究。喷管内的传热过程非常复杂，涉及热传导、热对流等多种传热方式，而且喷管的热状态不仅与其工况有关，还与喷管的结构、材料物性等有关，所以应当将流场与结构温度场耦合起来进行计算。张晓光等将软件Fluent与ANSYS结合起来，完成了喷管喉衬流场及热结构的模拟分析[2]。孙林、郝雯等建立了二维轴对称喷管模型，采用流固耦合方法，以喷管内壁面温度分布为边界条件，研究了喷管在旋转状态下的流动及换热情况[3-4]。白俊华等采用经过实验验证的流固耦合换热模型，对影响喷管喉衬热交换的某些因素进行了研究[5]。吴川、王启凡等同样采用流固耦合的方法，借助流体力学软件Fluent对喷管的结构温度场开展了瞬态仿真[6-8]。

目前，对喷管热防护层进行数值计算时都是将喷管内的纯气相流场结果作为热边界条件，但是对于某些特殊场合，比如应用于固体火箭发动机这一领域，由于高能推进剂的使用，使得燃烧产物的温度逐渐增大，并且内部伴随着颗粒的产生，而颗粒对气相的影响不可忽略，必会导致两相流场结果与纯气相不同，此外颗粒还会对喷管壁面造成严重的冲刷及烧蚀，给喷管的热防护带来了更大的挑战，所以有必要对两相流动产生的影响进行更加深入的分析与研究，以便进一步开展喷管的热分析与结构设计优化工作。

故本文为了提高数值计算的精度，建立了含热防护层的喷管两相流动流固耦合[9]换热模型，对拉瓦尔喷管内部流动换热进行了研究，并将喷管内两相流场计算结果作为结构温度场模拟的热边界条件，以获得喷管热防护层以及壳体各个时刻的温度分布情况，为喷管热分析提供依据。

2 数值模型

2.1 物理模型

鉴于本文中喷管的对称性，建立二维轴对称数值模型，如图1所示。喷管内流场部分设为流体域，其它各部分都为固体域。固体区域材料性能参数如表1所示，气相密度按理想气体处理，R=320J/(kg·K)，γ=1.2216，导热率按kinetic-theory给定[10]，粘性按Sutherland定律处理(三系数形式)[11]；计算区域都采用四边形结构网格进行划分，为了保证结果的精确性，在流固耦合交界面和喷管喉部进行局部加密，网格划分情况如图2所示。

图1 喷管计算模型示意图

表1 材料物性参数

图2 计算模型网格划分

为了便于分析，可对模型进行一些简化，做出如下假设：

1)忽略内壁面的烧蚀与炭化；

2)忽略各层材料之间的接触热阻；

3)忽略辐射传热和颗粒接触传热；

4)忽略燃气的化学反应；

5)假设颗粒内部和颗粒之间没有物理和化学作用；

2.2 控制方程

本文模拟包含两个过程，由于流场的发展速率远大于喷管固体结构的温度传递速率，因而，在进行喷管结构的热传导计算时，喷管流场都处于稳定状态。所以在计算时先进行喷管两相流的稳态数值模拟，然后将其结果当作热边界条件与喷管固体壁面进行瞬态耦合计算。喷管两相流计算采用Eulerian方法，对于湍流多相流来说，其控制方程包括流体相与颗粒相的连续方程、动量方程及能量方程，由于方程较多，动量方程与能量方程见参考文献[12]，连续性方程如下所示

1)流体相连续方程

(1)

颗粒相连续方程

(2)

式(1)(2)中符号和含义见参考文献[12]。

为了使方程封闭，还需加上湍流输运方程，如下所示：

3)流体相湍动能方程k

(3)

颗粒湍能输运方程kk

(4)

流体湍流耗散方程ε

(5)

式(3)(4)中符号和含义见参见文献[12]，式(5)中符号和含义见参考文献[13]。

固体区域仅考虑导热，其控制方程如下

(6)

式(6)中符号和含义见参考文献[5]。

2.3 计算方法与边界条件

为了更好的体现以两相流计算结果和纯气相计算结果作为热边界条件的不同，本文采用拟流体模型(PFM)来模拟喷管内的两相流动，对颗粒相与气相均采用Eulerian处理，该模型考虑了颗粒相的粘性、导热及扩散，相比于离散相(DPM)模型，可以全面考虑颗粒相的湍流输运。湍流模型选用RNG k-ε两方程模型，采用标准壁面函数(standard wall function)求解近壁区域物理量。喷管进口为压力进口边界，总压为7MPa，总温为3000K，颗粒体积分数为0.001，气相按理想可压流处理；出口取压力出口条件；喷管内流场与固体域接触表面设为耦合边界条件；外壁为对流换热边界，h=10 W·m-2·K-1。计算时间5s，时间步长取0.001s。

3 计算结果分析

3.1 模型验证

为了验证本文模型的正确性，对文献[14]中的JPL(Jet Propulsion Laboratory)喷管进行了纯气相模拟计算，喷管参数及边界条件见文献[15]，将本文结果与文献[14]中实验值进行了对照，如图3所示，能看出本文结果和实验值较吻合。

图3 模拟结果与实验结果对比

3.2 两相流场计算结果分析

本文先对喷管内的两相流场进行了模拟，图4是喷管内两相流场与纯气相流场计算结果的对照情况，其中图4(a)是两种工况下喷管内气相压强沿轴向变化曲线，从图中可以看出，颗粒相的存在对于喷管内流场压强分布的影响并不明显，这是由于颗粒不会像气体一样膨胀做功导致的。图4(b)为两种工况下喷管内气相温度沿轴向变化曲线，从图中可以看出，加入颗粒相之后，同一位置处燃气温度要大于纯气相情况下的温度，而且在喷管喉部和扩张段位置处颗粒相对气相温度的影响比较明显，这是因为颗粒的比热容较大，容易保持原来的温度，而且在喷管喉部之后燃气由于膨胀其温度减小，与此同时，颗粒相与气相发生热交换而使得气相温度增加。但在喷管喉部之前颗粒相对气相温度的影响并不明显，是因为这些位置燃气温度较高，颗粒相与气相的热交换较小。图4(c)(d)分别是喷管内气相速度与马赫数沿轴向的变化曲线，由这两个图可以看出，两相流情况下相同位置处燃气速度与马赫数都比纯气相时要小，这是因为颗粒具有较大的惯性，存在速度滞后对气相产生了阻力导致气流流速减小。从以上可以看出，两相流情况下的计算结果与纯气相时有较大不同，所以在进行喷管结构温度场数值计算时应该以两相流流场计算结果作为瞬态热传导模拟的初始值，而本文正是如此。

图4 两相流计算结果与纯气相计算结果对比

3.3 结构温度场结果分析

3.3.1 瞬时温度分布

图5是工作初期喷管喉部附近的温度云图，从图中可以看出，由于喉衬材料导热系数较大，炭/炭喉衬的温度传导速率明显高于烧蚀层和绝热层，喉衬与烧蚀层和绝热层的温差较大。而且时间越往后，可以发现炭/炭喉衬处热量向两侧烧蚀层传递，它们的温差逐步减小。图6是t=2，3，4，5 s时喷管的温度场分布情况，如图6所示，随着时间的推移，喷管壁面温度逐渐增加，热量在不断地往固体域内部传递，气体与内壁的对流换热和固体域内的热传导，使得喷管热防护层在径向呈现出清晰的温度梯度，由于烧蚀层、绝热层等的材料不同，其热力学性能有较大差异，使得各固体域的温度分布不同，炭/炭材料喉衬导热性能好，而绝热层和烧蚀层材料导热系数较小，温度上升缓慢，由于绝热层对热传递的阻隔，壳体温度也上升缓慢，靠近喷管内壁面的地方温度较高，喷管壁面的热量传递呈二维特性。

图5 工作初期喷管喉部附近温度云图

为了进一步分析喷管径向温度分布，分别在图7所示位置取截面，5s时截面上的热防护层和壳体温度分布如图8所示，横坐标均已无量纲化，横坐标为0的点是流体域与固体域的交界面位置处。由图8可以看出，每个截面的温度均沿径向逐渐降低，在炭/炭喉衬与绝热层以及烧蚀层与绝热层交界位置处曲线发生了转折，曲线斜率不同，但由于烧蚀层与绝热层材料导热系数相差较小导致图8(a)中效果不是太明显，而炭、炭喉衬与绝热层材料导热系数相差比较大，温度曲线斜率相差很大。壳体材质为钢，从表1可以知道其导热系数比喷管其它部分大得多，所以壳体温度保持在大于300K的小范围温度区间，图8中壳体位置的温度曲线斜率接近为零。另外，从图中可以看出，喷管内壁面处收敛段截面温度最高，然后是喉部的三个截面，左边截面温度最大，右边截面最小，且扩张段截面温度小于喉部截面温度，入口段截面温度小于收敛段截面温度。说明喷管内壁面最高温度位于收敛段中部至喉部前半段这一范围内，从图6可以看出，这是因为该区域容易积聚大量颗粒，受热严重，且颗粒的冲刷会对壁面造成不好的影响，加剧了该区域的烧蚀。

图6 不同时刻喷管温度场分布云图

图7 喷管各部位截面位置

图8 5s时喷管各部位截面温度曲线

3.3.2 对比分析

喉部是喷管的关键部位，且所处热环境恶劣，所以有必要考虑喉部温度分布及变化情况。图9是t=1，3，5 s时以两相流和纯气相流场计算结果为初始值的喷管喉部截面2上热防护层和壳体的温度分布曲线，图中横坐标原点为喷管内流场与固体域的交界面位置，横坐标为1的点为喷管外壁面位置，横坐标都已无量纲化。从图中可以看出，无论两相流情况下还是纯气相情况下，随着时间的推进，热防护层温度逐渐升高，在各个时刻，两相流情况下，喷管喉部截面2上沿径向温度分布比纯气相条件下要高，但壳体位置效果还不太明显。出现这种现象是因为颗粒比热容较大，容易保持原来的温度，而气相在喉部时已开始膨胀降温，由于温差，喷管内颗粒与燃气之间必会发生热交换，使得同一位置处气相的温度提高，同时粒子、燃气与壁面之间也存在着大量的热传递，进而热量以导热方式向热防护层传递。颗粒的存在对喷管热防护层、壳体的温度分布产生了较大影响，使得其温度值明显提高，这对喷管的安全是不利的。因此，在对喷管进行热分析时，不能忽略颗粒相的影响。

图9 不同时刻喷管喉部温度曲线对比

4 结论

相较于以往文献采用纯气相流场结果进行后续数值计算，本文不仅采用了流固耦合的方法，并考虑了颗粒的影响，以喷管两相流场结果作为固体壁面热传导初始条件对拉瓦尔喷管内的流动换热情况进行了模拟研究，得到了以下结论：

1)发动机工作过程中，燃气热量由喷管内壁向外壁传递，固体壁面温度逐渐升高，形成明显的径向温度梯度，同时喉衬导热系数大，不断往两侧传递热量，喷管壁面的热量传递呈二维特性。

2)喷管固体壁面最高温度位于内壁面收敛段中部至喉部前半部分这一区间内，此处容易积聚大量颗粒，在喷管中也最为脆弱。

3)颗粒的存在使得喷管内流场结构发生显著变化，导致喷管热防护层温度升高，对热防护层材料性能要求更高，因此，在进行喷管热分析时应考虑颗粒的作用，从而更加精确地对喷管热防护层进行设计。