三维超音速反应流高效隐式算法研究

三维超音速反应流高效隐式算法研究

王光利 关喜峰 王康

航空工业西飞 陕西西安 710089

摘要：发展了一种模拟超音速反应流的新的计算方法。采用一种解耦的隐式Gauss­-Seidel方法求解雷诺平均Navier-Stokes方程和组分方程，使用点隐式方法处理化学源项刚性问题。该隐式方法仅需对Navier-Stokes方程的矩阵求逆，组分方程矩阵则化为对角阵。该方法能够很容易应用于现有求解完全气体或者平衡理想气体的求解器上。应用本文方法对 “double appllo disk” 高超音速绕流进行了数值模拟并和实验结果进行了比较。

关键词：解耦隐式方法；超音速；三维Navier-Stokes方程；组分方程；化学非平衡

1 引言

气体化学非平衡效应对飞行器的影响是高超声速飞行器设计中的重要一环，这使得发展高效化学非平衡数值模拟算法成为需要。化学非平衡流的计算需要注意两个方面：数值稳定性和计算效率。稳定性问题主要由组分连续方程中的化学源项刚性引起；随着组分的增加，控制方程变得十分庞大，要求高效率的计算方法。国内外多数化学非平衡数值模拟方法多数计算速度低下，程序编写和修改工作十分繁杂，不利于程序的工程应用。

本文在参照了国内外相关工作经验基础上，结合早期解耦算法和目前常用的点隐式算法提出了一种新的隐式解耦算法^[1]。该数值模拟结果表明本文方法可以有效解决化学非平衡流动中的刚性问题，公式简单，程序编写和修改方便，易于实现和改进。采用本文程序对“双阿波罗盘”高超音速流场进行了数值模拟，并与实验进行了比较。

2 控制方程

本文应用的控制方程是三维可压缩非定常N-S方程及组分连续方程，其直角坐标系下的积分形式为：

（1）

假设所有组分气体为理想气体，热力学混合气体状态方程为：

（2）

（3）

公式所涉及到的参数以及化学源项的求解详见文献^[2]。

3 隐式高斯-赛德尔算法

控制方程可以离散成如下的形式：

（4）

可以分裂为：

； （5）

这里, 为雅可比矩阵的特征值。

令

（6）

则离散控制方程可以改写为：

（7）

上式中的上标 表示高斯–赛德尔迭代的扫描次数。一般说来，只需运用高斯–赛德尔迭代在全场网格上进行预估-校正两次扫描就可以了。

本文还采用了当地时间步长、隐式残值光顺^[3]等加速收敛措施，结合化学源项处理方法，本文算法可用较大的时间步长推进，较快地收敛到定常解。

4 点隐式方法

对于组分连续方程，还需要采用点隐式方法^[4]对化学源项进行处理，方程（7）可以改写为：

（8）

其中 是源项的雅可比矩阵，0≤ ≤1，由于本文针对定常流动进行数值模拟，因此对源项采用一阶精度进行计算， 取1。

为了简化计算，对源项雅可比矩阵进行对角化^[5]：

（9）

其中

（10）

（11）

式中 和 为化学反应的正反应速率， 和 为反应物和生成物的化学当量计算系数。

5 解耦方法

本文分别求解组分连续方程和N-S方程两个方程组，然后通过能量关系式将它们耦合起来。具体步骤如下：

（1）对流场进行初始化。

（2）采用本文中介绍的空间离散和时间推进方法求解N-S方程组，每推进一个时间步可以得到控制体中流体的密度，速度以及总内能。与此同时组分的质量分数，比热以及混合气体常数使用上一时间步的值。

（3）求解组分连续方程组，采用与求解N-S方程组相同的空间，时间格式和统一时间步长推进，每推进一个时间步，可以计算出各组分密度。

（4）由于组分连续方程和N-S方程是分别求解，因此需要将两个方程耦合起来。本文将求解组分连续方程所得到的各组分密度以质量分数的形式代入内能表达式：

（12）

再将状态方程代入内能表达式得：

（13）

通过以下的牛顿迭代法^[6]求解方程（13）可以得可以到温度，进而由状态方程求得压强，从而完成两个方程组的耦合。

6 算例和结果分析

采用本文程序对“双阿波罗盘”的高超音速流场进行了数值模拟，并与相关文献和实验结果^[7]进行了对比。

“双阿波罗盘”的上下是对称的形状，其轮廓线直接测绘自阿波罗登月指令舱的下层热防护罩。模型的直径为6英寸(1英寸等于2.54厘米），在直径24英寸的RPI高超音速激波风洞内进行了实验，模型表面安装了八个压电式压力传感器^[7]。

本文计算条件为：马赫数Ma=10，雷诺数Re=1.835×105，来流温度T∞=226.65K，绝热非催化壁。

下面给出了二维和三维化学非平衡数值模拟结果，来流假设为化学反应完全气体，采用5组分17反应化学动力学模型。

图1是二维数值模拟得到的流场速度矢量图，可以看到，在弓形激波前后流动速度急剧减小，而流动绕过物体以后又开始加速（死水区除外），在背风区形成了两个基本对称的涡。

图2给出了化学组分沿驻点轴线的分布。可以看到，在波后区域氧分子和氮分子开始分裂，与此同时氧原子，氮原子以及一氧化氮开始增长，从激波前的冻结流状态（几乎为零）变化一直达到平衡状态。我们注意到各组分到达驻点处基本达到平衡，这可能是因为流体单元沿着驻点流线速度逐渐减小到零，这给各组分达到化学平衡提供了时间。

图3是物面压力和实验值的比较，可以看到，二维和三维计算结果和实验值都吻合得较好。图4展示了通过纹影法得到的实验照片和数值模拟得到的压力云图。左图是二维计算结果和实验的比较，右图是三维计算结果和实验的比较，图的上半部分是压力云图，下半部分是实验照片。通过比较可以看出，三维计算结果和实验符合较好，二维计算结果激波厚度明显偏大。这可能是因为在二维数值模拟中，流体微粒绕过物体时只能沿着一个平面运动而不能流向空间其它方向所致^[7]。

7 结论

发展了一套新的数值模拟高超音速化学非平衡流动计算程序。通过求解三维可压缩非定常N-S方程及组分扩散方程来数值模拟绕“双阿波罗盘”超音速化学非平衡流场。计算采用一种解耦的隐式Gauss­-Seidel方法求解，使用点隐式方法处理化学源项的刚性问题。本文方法只需对N-S方程雅可比矩阵求逆，组分方程矩阵则化为对角阵。该方法能够很容易应用于现有求解完全气体或者平衡理想气体的求解器上。计算结果与实验结果符合较好。

图1 速度矢量图 图2 化学组分沿驻点轴线的分布

图3 物面压力分布 图4 激波脱体距离比较

参考文献

1 张黎,叶正寅,王刚, COIL化学非平衡流动的一种解耦方法研究,空气动力学学报, 2008,26(2).

2 Robert F. Tomaro, William Z. Strang, An Implicit Algorithm for Solving Time Dependent Flows on Unstructured Grids[R], AIAA-97-0333, Reno: AIAA 1997.

3 Wang G. Strategy for the generation of tetrahedral unstructured viscous grid on three-dimensional configurations and solving the N-S equations. Master dissertation. Xian: Northwestern Polytechnical University, 2001.

4 H.C. Yee,Judy L. Shinn, Semi-implicit and fully implicit shock-capturing methods for hyperbolic conservation laws with stiff source terms , AIAA paper, 1987, 1987-1116.

5 Eberhardt, S., Imlay, S., Diagonal Implicit Scheme for Computing Flows with Finite Rate Chemistry, J. Themphysics and Heat Transfer, Vo1.6, 1992.

6 Miller J H , Shang J S , Paschkewitz J S . A Parallelj Unstructured Flow Solver for Chemically Reacting COIL Flowfields. AIAA paper, 2001, 2001-1089.

7 Toro, P.G.P, Korzenowski H, Minucci M.A.S, Nagamatsu H.T and Myabo L.N, Numerical and Experimental Investigations of "Double Apollo Disc" Inlet at Mach 10, AIAA 2000-3490, 36th A1AA/ASME/SAE/ASEE Joint Propulsion Conference, AL.

作者简介：王光利(1978-)，男，四川成都人，工程师。

──────────────────