COMSOL 博客

如何在 COMSOL Multiphysics® 中模拟超音速流动

作者 Siva Sashank Tholeti

2017年 8月 10日

你是否曾经对如何模拟超音速流动，例如协和式飞机或战斗机周围的空气流动感到好奇？通常，这个过程需要求解流动中的激波或膨胀波。求解不连续现象（如激波）需要相当精细的网格，以及对流体流动的质量、动量和能量守恒强耦合方程采用合适的数值稳定性算法。在这篇博文中，我们将讨论如何模拟通过菱形翼型的超音速流动，其中涉及求解激波和膨胀波。

亚音速、音速和超音速流动之间有什么区别？

考虑一个在空中以恒定速度 u 运动的机翼。由机翼产生的扰动在空气中以压力波的形式传播，通常表现为以该处的音速 a 传播的声波。如果机翼的速度小于音速，即 u<a，则扰动的传播速度比机翼本身快。因此，上游的流体在机翼尚未到达该位置之前就已经受到了机翼存在的影响，由于机翼的速度低于音速，所以通常称为亚音速流。如果机翼的速度增加，甚至超过了音速，在机翼到达上游流体位置之前，上游流体不会受到机翼的影响，因为压力波还没有传播到那里，如下图所示（红色马赫锥）。扰动传播使我们能够定性地了解亚音速和超音速流动状态之间的分界，这种现象也可以用马赫数（Ma =u/a，无量纲）来表示。

比较亚音速、音速和超音速机翼速度的扰动传播的示意图。
亚音速(左)、音速(中)和超音速(右)机翼速度的扰动传播。其中圆圈代表压力波/声波，箭头代表物体的运动方向。

为了定量地求解流体流动问题，我们首先要确定流动是处于层流状态，还是已经转捩到湍流状态。这通常是基于雷诺数 Re =ρuL/μ 确定的，其中，ρ、μ 和 u 分别代表流体密度、动力黏度和速度，L 则表示特征长度(对于机翼，使用弦长)。对于层流，求解具有连续性和本构关系的纳维-斯托克斯方程；对于湍流，求解具有连续性和本构关系的雷诺平均纳维-斯托克斯方程。根据问题的不同，可以选择合适的湍流模型来求解湍流问题。

如果流动变得高度可压缩(比如超音速流动)，除了前面提到的质量和动量守恒方程外，还必须求解能量方程，无论流动是处于层流还是湍流状态。超音速流动中的黏性效应通常可以忽略不计，除非在具有急剧梯度的区域，例如激波或边界层。但是，如果感兴趣区域中的黏性效应可以忽略不计，那么删除纳维-斯托克斯方程中的黏性项就可以得到欧拉方程。对于简单的几何形状，这些无黏性流动方程可以得到解析解。

菱形翼型的超音速气流示意图。
菱形翼型周围的超音速流场，显示流场中的不同区域。

为简单起见，考虑超音速流过带菱形翼型的机翼横截面的基准情况。已有 COMSOL Multiphysics^® 用户在其论文“Navier-Stokes Equations in Supersonic Flow over a Double Wedge Airfoil using Adaptive Grids”中研究了这个主题。如前所述，上游流体（上图中的 0 区）不受机翼引起的扰动的影响。因此，当流体接近翼型时，它面临着流动面积的突然减小（类似于凹角），并且必须突然改变方向以匹配1区中的边界条件。

这个过程只有在流动不连续的情况下才会发生，这种不连续性被称为斜激波，并且非常薄，也就是说，与气体分子的平均自由程有关，在大气压力下空气分子的平均自由程约为 50 nm。当激波由于物体的形状和流动的马赫数而向流动方向倾斜时，称为斜激波。在激波过程中，静压、温度和密度增加，而马赫数和总压降低。

计算激波的激波角和马赫数

如上所述，当我们假设无黏流动（这通常是估计超音速流动中流动特性的有效假设）时，可以根据物体的倾角确定弱激波的激波角(σ)。在本例中，这是菱形翼型的半角 (δ) 和上游马赫数 (Ma₀)。

\tan\delta = 2 \cot\sigma \left[ Ma_{0}^2 \sin^2\sigma -1 \over {Ma_{0}^2 (\gamma + \cos 2\sigma) + 2}\right]

对于比热比为 γ 的流体，1 区斜激波后的马赫数可以由下式估算：

Ma_1 = \sqrt{(\gamma+1)^2 Ma_0^4 \sin^2\sigma-4(Ma_0^2 \sin^2\sigma -1)(\gamma Ma_0^2 sin^2\sigma +1) \over [2\gamma Ma_0^2 \sin^2\sigma – (\gamma-1)][(\gamma -1)Ma_0^2 \sin^2 \sigma + 2]}

在沿着菱形翼型的激波之后，流动遇到翼型顶部周围的膨胀区域（类似于一个凸角）。流动方向的变化以匹配边界条件是通过膨胀波实现的，也称为 Prandtl-Meyer 膨胀波。在膨胀过程中，静压、温度和密度降低，而马赫数增加。同样，通过假设无黏流动，我们可以根据Ma₁ 和转向角(θ)估计膨胀波后的马赫数：

\theta = \sqrt{\gamma +1 \over \gamma -1} \left[ \tan^{-1} \left( \sqrt{{\gamma -1 \over \gamma + 1}(Ma_2^2 -1)}\right) – \tan^{-1} \left( \sqrt{{\gamma -1 \over \gamma + 1}(Ma_1^2 -1)}\right)\right]- \left[\tan^{-1}\left(\sqrt{Ma_2^2-1}\right) – \tan^{-1}\left(\sqrt{Ma_1^2-1}\right) \right]

基于上述方程，我们可以估计激波或膨胀波的冲激波角和马赫数，也可以计算压力、密度和温度等类似估计值。然而，对于冲击/膨胀波相互作用的更复杂的几何形状或高黏性流体，通过整体守恒方程的数值解来求解流动会更合适。接下来，我们使用 COMSOL 来对菱形翼型这个简单基准案例进行建模，并将结果与上述方程的估计值进行比较。

在 COMSOL Multiphysics^® 中模拟通过菱形翼型的超音速流动

如前所述，这种情况的流动在超音速下变得高度可压缩，导致所有守恒方程的强耦合。现在必须一起求解这些方程，以求解流场中出现的激波和膨胀波，这可以使用 COMSOL 软件 CFD 模块中的高马赫数流动 接口来完成。让我们看一下如何使用高马赫数流动 接口建立通过菱形翼型的超音速流模型。

菱形翼型模型的示意图，包括边界条件。
模型配置，包括边界条件。

模型设置的示意图，以及入口处的静压和温度如上所示。根据入口条件和机翼弦长，发现雷诺数大于 3 x 10⁵，表明流动处于湍流状态。因此，这里使用了 k-ε 湍流模型。

在出口边界条件下使用混合流动条件，这个条件允许该边界处的流动为超音速或亚音速。如果流动是亚音速的，那么将确保条件设置下指定的静压。在翼型上定义了一个无滑移条件，并在顶壁和底壁上使用滑移壁和出口（混合流）边界条件的组合，用于模拟解析求解的无黏模型方程。此外，由于这肯定是一个非线性问题，需要注意的是，我们可以通过指定可能接近最终解的初始条件来实现问题的更好收敛，正如之前关于非线性静态有限元问题的博客中所讨论的那样。

如前所述，激波是流场中的一种不连续形式。因此，必须在压力波处和周围使用非常精细的网格，以便在数值上解析这个区域。激波或膨胀波的空间位置无法提前得知。此外，我们可能对计算不同攻角和不同入口条件下的流场感兴趣，这将导致域中不同位置的冲击。因此，如果对这类问题使用均匀网格，则需要在任何地方都具有非常高的网格密度，从而大大增加计算成本。但是，在 COMSOL Multiphysics^® 软件中，我们可以利用自适应网格细化（求解器中提供的选项）动态调整域中的网格细化，以便在不连续处附近创建细网格，并在其余域中使用相对较粗的网格。在这里，自适应网格细化的条件是在稳态研究步骤的设置下定义。

下表列出了研究结果和上一节中给出的表达式的结果，以进行比较。COMSOL Multiphysics 的结果是在流线区域的中间部分进行评估的，如下面的流线图所示。如表格所示，仿真结果和理论估计值相差在 3% 以内。值得注意的是，数值模型考虑了问题中存在的任意黏性效应，这使得数值模型比基于无粘性理论的解析估计在物理上更准确。