物理模型
这篇文章是我学习《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的时候（也就是初始值），密度、温度、速度场如下：
实验结果
代码在此：
这份代码会先输出亚音速-超音速算例的结果。
书上的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等。
代码还是这份代码：
它应该在输出第一个算例之后，输出第二个算例迭代5000步后的结果，也就是书中的table 7.7，请自行验证：）