基于MacCormack方法的Quasi-1D Isentropic Nozzle Flow数值模拟

物理模型

这篇文章是我学习《COMPUTATIONAL FLUID DYNAMICS: The Basics with Applications》(作者John D. Anderson, Jr.)的心得体会。

这里的Quasi-1D Nozzle Flow是书上第七章的内容:用MacCormack算法求解拉瓦尔喷管(convergent-divergent nozzle)中的流场。这种喷管的示意图如下:

这种喷管的横截面是漏斗状,先变窄再变宽(这就是’convergent-divergent’的由来),它可以把亚音速气流(从左端进入)加速到超音速(从右端流出),在最窄处恰好是音速。原因是亚音速气流和超音速气流的性质不同:亚音速流截面越小流速越快,超音速流则截面越大流速越快。

在这个问题里面,我们求解二维的问题,并假定流场的性质只和x坐标有关(也就是说在上图中,一‘竖条’上的流体属性是相同的),这就是“Quasi-1D”的含义。当然,这种假设只是对现实的一个近似,目的是计算方便。

微分方程

既然是一维求解,就需要用到一维的方程。假定流场内是理想气体,故可以从方程中消去压强场,剩下三个未知场:密度、速度(x方向)、温度。书上有相关推导,最终的微分方程是(7.42)~(7.44):

在这一章节中,我们采用无量纲量进行计算,以消除物理单位。各个参数对应的无量纲量定义在书的第297页:

其中gamma是比热容(对于空气,取1.4),下标为0代表在喷管的入口处(最左边),A*是喷管最窄处的面积。

经过推导,可以得到对应无量纲量的微分方程,即(7.46),(7.48),(7.50).

MacCormack

当然,计算机程序无法直接求解微分方程,必须设法将其转化。在这里我们采用MacCormack方法。这是一种explicit方法。所谓explicit,就像(7.42)~(7.44)一样,每一个方程里恰好含有一个关于时间的导数,这样在递推时就可以显式(explicitly)计算出来每一个场关于时间的差分,然后递推。

但直接用差分乘以时间步的精度是一阶的(它对应了泰勒展式中的头两项),这就是用到MacCormack方法的原因:它可以在只处理一阶导数的情况下将explicit方程的求解精度提升至二阶。

MacCormack方法的步骤叫做predictor-corrector step:先按照上述的一阶方法求出来下一个时刻的“预测流场”(predictor),然后用“预测流场”算出来一个“预测时间差分”,再回过头,把“预测差分”和在本时间步求得的时间差分做平均,用这个平均值去求下一个时间步的流场,作为最终结果(corrector).有一个关键点:在第一次用现有流场算时间差分时需要用forward difference(就是(f[i+1]-f[i])/dx),第二次用预测流场算时间差分时则用rearward difference(就是(f[i]-f[i-1])/dx),或者第一次用rearward,第二次用forward.

MacCormack方法是《The Effect of Viscosity in Hypervelocity Impact Cratering》,Robert W. MacCormack提出的,该论文里称用傅里叶级数进行分析可以证明,它的精度是二阶。

explicit方程本身在书中有推导.用forward difference的(predictor step)是(7.51)~(7.53):

用rearward difference的(corrector step)是(7.57)~(7.59):

变量上面加横线代表这个是“预测流场”。

时间步长

对于explicit方法,第四章的分析表示,要想获得稳定的解,时间、空间步长必须满足关系:

C是Courant number(这个理论叫Courant–Friedrichs–Lewy condition),它不应该超过1,在本例中我们取0.5.注意到当地音速a和流体速度V是随着流场改变的,对于本例,我们取计算出的所有dt的最小值。

算例1:亚音速-超音速

边界条件

书中称,由于喷管最窄处左边是亚音速流场,也就是说扰动会向左传播,直至喷管的入口,所以只能固定喷管入口处的两个变量(密度和温度),但允许速度值自由“浮动”。严谨的数学分析“超越了本书的范围”。

对于入口处的速度值和喷管出口处的所有值,由于forward difference和rearward difference总有一个在该处无法计算,所以它们的值无法直接从差分方程中获得。对于这些地方,我们运用linear extrapolation(假设每个值随x坐标线性):

 

我们规定入口处的密度值和温度值是固定的常数。由于采用了无量纲量,它们就固定为1:

参数设定

在这个算例中,我们规定喷管的长度是3L,截面函数是:

它在1.5L的地方达到最窄.

x方向的步长是0.1,这意味着一共31个格点。

在t=0的时候(也就是初始值),密度、温度、速度场如下:

实验结果

代码在此:

https://github.com/wmdcstdio/DigitalWindTunnel/blob/master/NozzleFlowQuasi1D/NozzleFlowQuasi1D/main.cpp

这份代码会先输出亚音速-超音速算例的结果。

书上的Table 7.1是初始的各个参数值,Table 7.3是进行1400次迭代后的各个参数值。可以发现,进行1400次迭代后,马赫数在格点16和17(最中间的两个格点)之间跨越了1,说明在最窄处流速确实恰为音速。

我的代码计算出的数值和书上结果完全一样,除了几处printf的舍入误差。限于篇幅(懒),我就不贴具体数据了,有兴趣者请自己下载代码运行:)

算例2:纯亚音速

在这个算例中,管子的形状稍有不同:

纯亚音速的边界情况和上面讨论的亚音速-超音速略有不同。原因是解不再唯一,实际上喷管出入口的压强比值和流场方程的解一一对应。因此我们需要规定这个比值,根据标准气体方程,实际上就是规定出口处[latex]\rho’T'[/latex]的值(为一个常数,在这个算例中等于0.93)。我们可以extrapolate出其中的一个(比如[latex]\rho[/latex]),然后用[latex]\frac{常数}{\rho}[/latex]求出另一个。

这个算例的初始值也不同:

其他和亚音速-超音速的算例相同,包括入口处的extrapolation等。

代码还是这份代码:

https://github.com/wmdcstdio/DigitalWindTunnel/blob/master/NozzleFlowQuasi1D/NozzleFlowQuasi1D/main.cpp

它应该在输出第一个算例之后,输出第二个算例迭代5000步后的结果,也就是书中的table 7.7,请自行验证:)

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注