14.4 The MUSCL-Hancock Scheme

14.4 The MUSCL-Hancock Scheme

通过MUSCL格式能够实现二阶精度。

14.4.1 The Basic Scheme

非线性双曲守恒律：

\[\mathbf{U}_t + \mathbf{F}(\mathbf{U}_x) = 0 \]

目标是为了构建二阶精度的通量\(\mathbf{F}_{i+\frac{1}{2}}\)。一般来说MUSCL格式包含以下步骤：

1. 数据重构。假设已经选择好了变量\(\mathbf{W}\)。可以重构守恒量\(\mathbf{W} = \mathbf{U}\)，当然也可以选择原始变量，对于欧拉方程原始变量为，\(\mathbf{W} = (\rho, u, p,)^{T}\)。 以下的格式均为守恒量重构。在重构的过程中，网格平均值\(\mathbf{U}_{i}^{n}\)被替换为局部分段线性函数：

   \[\mathbf{U}_i(x) = \mathbf{U}_{i}^{n} + \frac{(x-x_i)}{\Delta x} \Delta_i, \qquad x \in [0, \Delta x] \]

   其中，\(\Delta_i\)是\(\mathbf{U}_i(x)\)在网格中的斜率向量。局部左边中的极点\(x=0\)和\(x = \Delta x\)对应全局坐标系中的网格边界\(x_{i-\frac{1}{2}}\)和\(x_{i+\frac{1}{2}}\)。极点处的值，也叫作边界外插值：

   \[\mathbf{U}_{i}^{L} = \mathbf{U}_{i}^{n} - \frac{1}{2}\Delta_i, \qquad \mathbf{U}_{i}^{R} = \mathbf{U}_{i}^{n} + \frac{1}{2}\Delta_i \]

   注：对于欧拉方程来说，\(\mathbf{U}\)和\(\Delta_i\)是有三个component的向量，故需要六个数量。

2. 演化。在每个边界上的外插值\(\mathbf{U}_{i}^{L}\)和\(\mathbf{U}_{i}^{R}\)根据以下规则演化\(\frac{1}{2}\Delta t\)时间：

   \[\bar{\mathbf{U}}_{i}^{L} = \mathbf{U}_{i}^{L} + \frac{1}{2}\frac{\Delta t}{\Delta x}[\mathbf{F}(\mathbf{U}_{i}^{L}) - \mathbf{F}(\mathbf{U}_{i}^{R})] \\ \bar{\mathbf{U}}_{i}^{R} = \mathbf{U}_{i}^{R} + \frac{1}{2}\frac{\Delta t}{\Delta x}[\mathbf{F}(\mathbf{U}_{i}^{L}) - \mathbf{F}(\mathbf{U}_{i}^{R})] \]

   注：演化过程完全在一个网格内部。在每个界面\(i+\frac{1}{2}\)处，有两个通量\(\mathbf{F(\mathbf{U}_{i}^{R})}\)和\(\mathbf{F(\mathbf{U}_{i+1}^{L})}\)。

3. 黎曼过问题。为了计算界面的通量，需要计算以下黎曼问题：

   \[\mathbf{U}_{L} = \bar{\mathbf{U}}_{i}^{R}, \qquad \mathbf{U}_{R} = \bar{\mathbf{U}}_{i+1}^{L} \]

   以得到自相似解\(\mathbf{U}_{i+\frac{1}{2}}(x/t)\)。面上的通量与一阶迎风Godunov格式相同，即：

   \[\mathbf{F}_{i+\frac{1}{2}} = \mathbf{F}(\mathbf{U}_{i+\frac{1}{2}}(0)) \]

对斜率向量\(\Delta_i\)的一个可能的选择：

\[\Delta_i = \frac{1}{2}(1+\omega)\Delta_{i-\frac{1}{2}} + \frac{1}{2}(1-\omega)\Delta_{i+\frac{1}{2}} \label{slope} \]

其中：

\[\Delta_{i-\frac{1}{2}} = \mathbf{U}_{i}^{n} - \mathbf{U}_{i-1}^{n}; \qquad \Delta_{i+\frac{1}{2}} = \mathbf{U}_{i+1}^{n} - \mathbf{U}_{i}^{n} \]

且\(\omega \in [-1, 1]\)

14.4.2 A Variant of the Scheme

也可以在第一步和第二步中使用原始变量\(\mathbf{W}\)。对于欧拉方程，使用\(\mathbf{W} = (\rho, u, p)^{T}\)得到：

\[\mathbf{W}_t + \mathbf{A}(\mathbf{W})\mathbf{W}_x = 0 \]

此时第一步和第二步可以合为一步。

\[\bar{\mathbf{W}}_{i}^{L,R} = \mathbf{W}_{i}^{L,R} + \frac{1}{2}\frac{\Delta t}{\Delta x}\mathbf{A}(\mathbf{W}_{i}^{n})(\mathbf{W}_{i}^{L} - \mathbf{W}_{i}^{R}) \]

其中\(\mathbf{A}(\mathbf{W}_{i}^{n})\) 代表网格均值上的矩阵\(\mathbf{A}\)，最终结果：

\[\bar{\mathbf{W}}_{i}^{L} = \mathbf{W}_{i}^{n} - \frac{1}{2}[\mathbf{I} + \frac{\Delta t}{\Delta x}\mathbf{A}(\mathbf{W}_{i}^{n})] \Delta_i \\ \bar{\mathbf{W}}_{i}^{R} = \mathbf{W}_{i}^{n} + \frac{1}{2}[\mathbf{I} - \frac{\Delta t}{\Delta x}\mathbf{A}(\mathbf{W}_{i}^{n})] \Delta_i \]

此时求通量即为以下黎曼问题：

\[\mathbf{W}_{L} = \mathbf{W}_{i}^{n} + \frac{1}{2}[\mathbf{I} - \frac{\Delta t}{\Delta x}\mathbf{A}(\mathbf{W}_{i}^{n})] \Delta_i \\ \mathbf{W}_{R} = \mathbf{W}_{i+1}^{n} - \frac{1}{2}[\mathbf{I} + \frac{\Delta t}{\Delta x}\mathbf{A}(\mathbf{W}_{i+1}^{n})] \Delta_{i+1} \]

找到\(\mathbf{W}_{i+\frac{1}{2}}(x/t)\)并计算通量：

\[\mathbf{F}_{i+\frac{1}{2}} = \mathbf{F}(\mathbf{W}_{i+\frac{1}{2}}(0)) \]

14.4.3 TVD Version of the Scheme

MUSCL格式有二阶精度，根据Godunov 定理，在梯度比较大的区域会有很大的震荡。所以需要TVD格式，即把斜率\(\Delta_i\)替换为限制后的斜率\(\bar{\Delta}_i\)，有以下两种方式能够实现这一过程。

限制斜率 (Limited Slopes)

直接限制斜率：

\[\bar{\Delta}_i = \left\{\begin{matrix} \max [0, \min(\beta \Delta_{i-\frac{1}{2}}, \Delta_{i+\frac{1}{2}}), \min(\Delta_{i-\frac{1}{2}}, \beta \Delta_{i+\frac{1}{2}})], & \Delta_{i+\frac{1}{2}} > 0 \\ \min [0, \max(\beta \Delta_{i-\frac{1}{2}}, \Delta_{i+\frac{1}{2}}), \max(\Delta_{i-\frac{1}{2}}, \beta \Delta_{i+\frac{1}{2}})], & \Delta_{i+\frac{1}{2}} < 0 \end{matrix}\right. \]

当\(\beta = 1\)的时候得到MINBEE (MINMOD)格式， \(\beta = 2\)的时候得到SUPERBEE格式

对上述方法的一个改进是特征限制，也叫wave-by-wave limiting。首先，网格内的斜率\(\Delta_i\)是边界上跳跃\(\Delta_{i-\frac{1}{2}}, \Delta_{i+\frac{1}{2}}\)的函数：

\[\bar{\Delta}_i = \Delta_i(\Delta_{i-\frac{1}{2}}, \Delta_{i+\frac{1}{2}}) \]

在这个方法中，在特征场上分解上述的两个跳跃：

\[\Delta_{i-\frac{1}{2}} = \sum_{k=1}^{N} \Delta_{i-\frac{1}{2}}^{k}; \qquad \Delta_{i+\frac{1}{2}} = \sum_{k=1}^{N} \Delta_{i+\frac{1}{2}}^{k} \]

在分解后的每个子项中使用限制斜率的方式处理。

使用斜率限制器(Use of Slope Limiters)

另一种方式是找一个斜率限制器\(\xi_i\)使得：

\[\bar{\Delta}_i = \xi_i \Delta_i \]

其中斜率通过公式\(\eqref{slope}\)计算。这种方法得到斜率限制器\(\xi_i\)的TVD区域：

\[\xi(r) = 0; \qquad r <0 \\ 0 \le \xi(r) \le \min\{ \xi_L(r), \xi_R(r) \} ; \qquad r >0 \]

其中：

\[\xi_L (r) = \frac{2\beta_{i-\frac{1}{2}}r}{1 - \omega + (1+\omega)r} \\ \xi_R (r) = \frac{2\beta_{i+\frac{1}{2}}}{1 - \omega + (1+\omega)r} \\ r = \frac{\Delta_{i-\frac{1}{2}}}{\Delta_{i+\frac{1}{2}}} \]

且：

\[\beta_{i-\frac{1}{2}} = \frac{2}{1+c}; \qquad \beta_{i+\frac{1}{2}} = \frac{2}{1-c} \]

在单波方程中，\(c\)代表Courant number。考虑\(\beta_{i-\frac{1}{2}}\)和\(\beta_{i+\frac{1}{2}}\)的极值，简化后可以直接令\(\beta_{i-\frac{1}{2}} = \beta_{i+\frac{1}{2}} = 1\)。一种改进的方法可以采用特征极值。这种方法需要求解额外的黎曼问题来得到每道波的间断和相应的Courant numbers。

几种限制器：

superbee flux limiter：

\[\xi_{sb}(r) = \left\{\begin{matrix} 0 & \text{if } r \le 0\\ 2r & \text{if } 0 \le r \le \frac{1}{2} \\ 1 & \text{if } \frac{1}{2} \le r \le 1 \\ \min\{ r, \xi_R(r), 2 \} & \text{if } r \ge 0 \end{matrix}\right. \]

van Leer–type slope limiter:

\[\xi_{vl}(r) = \left\{ \begin{matrix} 0 & \text{if } r \le 0 \\ \min\{\frac{2r}{1+r}, \xi_R(r)\} & \text{if } r \ge 0 \end{matrix} \right. \]

van Albada–type slope limiter:

\[\xi_{va}(r) = \left\{ \begin{matrix} 0 & \text{if } r \le 0 \\ \min\{\frac{r(1+r)}{1+r^2}, \xi_R(r)\} & \text{if } r \ge 0 \end{matrix} \right. \]

minbee–type slope limiter：

\[\xi_{va}(r) = \left\{ \begin{matrix} 0 & \text{if } r \le 0 \\ r & 0 \le r \le 1 \\ \min\{1, \xi_R(r)\} & \text{if } r \ge 0 \end{matrix} \right. \]

14.4.4 Summary of the MUSCL–Hancock Method

总结：

1. 边界条件。
2. 计算时间步长。根据CFL数计算。
3. 边界插值，并演化。
4. 求解黎曼问题，计算通量。
5. 更新解。
6. 下一个时间步，回到第一步。