3.2.3. 高分辨率格式

(1) 抑制间断附近的数值振荡

二阶以上的 FDM 或 FVM 在间断附近的解可能会出现数值振荡，而一阶精度的格式通常会把激波、接触间断等 “抹平”。数值振荡是由于单调性的丧失， 保单调性 可以抑制间断附近的振荡。

一阶迎风格式和 Lax 格式是线性单调格式，因此不会出现数值振荡。根据 Godunov定理 : 线性单调格式最多能达到一阶精度。因此，保单调性的二阶或二阶以上精度的高分辨率格式必定是 非线性格式 , 即使求解的方程是线性的。其中的一种典型 高分辨率格式 是 TVD 格式。

TVD (Total Variation Diminishing) 格式：如果差分格式满足 时间推进 时总变差不增加，则称为 TVD 格式。谓总变差，是衡量函数的波动状况的一种指标，在离散点上，定义为: \(TV(u)=\sum_{j=-\infty}^{\infty} |u_{j+1}-u_j|\)

可以证明，非线性双曲型守恒方程组（系统）的解析解本身不具备 TVD 的特点。所以要求求解守恒方程的格式具有 TVD 的性质是不需要的，甚至是有害的。实际上，有人已经证明，不存在求解非线性守恒系统的 TVD 格式。所以所谓非线性守恒系统的 TVD 格式只是对标量方程TVD格式的形式上的推广。

TVD 条件是比较严的保单调条件, 不能直接用于多维问题。可以证明，对二维、三维情况，不存在二阶及以上精度的 TVD 格式（即使对于线性标量方程）。因此，将二阶 TVD 格式推广到多维问题，也只是形式上的推广。

极值原理 (Maximum Principal) 是可以直接应用于多维问题抑制振荡的准则。

(2) 一维 Euler 方程的 TVD 格式

Euler 方程的 TVD 格式，在形式上和有限体积框架下的 Roe 格式是类似的。主要区别在于重构界面左右端守恒变量 \(\mathbf{U}_{j+1/2}^L, \mathbf{U}_{j+1/2}^R\) 的算法。

(2.a) 简化算法

(3.2.51)\[\begin{split}\mathbf{D}_j &= \text{min mod} (\mathbf{U}_{j+1}-\mathbf{U}_{j}, &\mathbf{U}_{j}-\mathbf{U}_{j-1})/\Delta x \\ \mathbf{D}_{j+1} &= \text{min mod} (\mathbf{U}_{j+2}-\mathbf{U}_{j+1}, &\mathbf{U}_{j+1}-\mathbf{U}_{j})/\Delta x\end{split}\]

(3.2.52)\[\begin{split}\mathbf{U}_{j+1/2}^L &= \mathbf{U}_{j} + 0.5 \mathbf{D}_j \Delta x \\ \mathbf{U}_{j+1/2}^R &= \mathbf{U}_{j+1} - 0.5 \mathbf{D}_{j+1} \Delta x \\\end{split}\]

其中, \(\text{min mod}\) 算子为:

(3.2.53)\[\begin{split}\text{min mod}(a,b) = \left\{ \begin{array}{l} \text{sign}(a) \min(|a|, |b|) &, \text{if } ab > 0 \\ 0 &, \text{otherwise} \end{array} \right.\end{split}\]

在 TVD 格式中, 还可使用以下算子 \(B(a,b)\) 作为限制器。

(3.2.54)\[\begin{split}B _\text{Van Leer}(a,b) = \frac{ab[\text{sign}(a)+\text{sign}(b)]}{|a|+|b|+\epsilon} \\ B _\text{Van Albada}(a,b) = \frac{\max(ab,0)(a+b)}{a^2+b^2+\epsilon}\end{split}\]

TVD 格式是一类计算激波的数值方法，限制器并不限于 \(\text{min mod}\) 算子。

(2.b) 基于特征方向的 TVD 格式

完整的 TVD 格式应当根据 Riemann 问题的特征方向对 \(\mathbf{U}_{j+1/2}^L, \mathbf{U}_{j+1/2}^R\) 进行重构，而不是简单地施加限制器。

首先，计算 ROE 平均的特征值向量 \([\tilde{\lambda}]=[\tilde{\lambda}_1,\tilde{\lambda}_2,\tilde{\lambda}_3]^T\) (Eq.3.2.34) 和特征向量矩阵 \(\tilde{R}, \tilde{L}\) (Eq.3.2.35, Eq.3.2.36)。

之后，计算特征变量的梯度向量：

(3.2.55)\[\begin{split}\mathbf{D}_j &= \text{min mod} (\tilde{L}(\mathbf{U}_{j+1}-\mathbf{U}_{j}), &\tilde{L}(\mathbf{U}_{j}-\mathbf{U}_{j-1}))/\Delta x \\ \mathbf{D}_{j+1} &= \text{min mod} (\tilde{L}(\mathbf{U}_{j+2}-\mathbf{U}_{j+1}), &\tilde{L}(\mathbf{U}_{j+1}-\mathbf{U}_{j}))/\Delta x\end{split}\]

再计算 Riemann 问题的左右状态向量 (Eq.3.2.37):

(3.2.56)\[\begin{split}\mathbf{W}^L &= \tilde{L}\mathbf{U}_{j} + \frac{1}{2}(\mathbf{D}_j\Delta x - \mathbf{D}_j \odot [\tilde{\lambda}] \Delta t) \\ \mathbf{W}^R &= \tilde{L}\mathbf{U}_{j+1} - \frac{1}{2}(\mathbf{D}_{j+1}\Delta x + \mathbf{D}_{j+1} \odot [\tilde{\lambda}] \Delta t)\end{split}\]

其中, \(\odot\) 是哈达马积 (Hadamard product), 指两个维度相同矩阵对应位置相乘。是 Kronecker 积 (直积, 张量积) \(\otimes\) 在两矩阵维度相同时的简化形式。

最后，计算左右守恒向量:

(3.2.57)\[\begin{split}\mathbf{U}_{j+1/2}^L &= \tilde{R} \mathbf{W}^L \\ \mathbf{U}_{j+1/2}^R &= \tilde{R} \mathbf{W}^R \\\end{split}\]

原则上，可用任意时间离散格式，但 Shu-Osher 证明了如果用 TVD Runge-Kutta 格式 (Eq.3.1.24) 进行时间离散，得到的全离散格式也具有 TVD 性质。

以下为 TVD 格式通量在一维 Euler 问题中的实现，为了更清晰，代码的实现不考虑效率。

    @staticmethod
    def Upwind1_TVD(Ujm1: np.ndarray, Uj: np.ndarray, 
                    Ujp1: np.ndarray, Ujp2: np.ndarray,
                    limiter) -> Tuple[np.ndarray, np.ndarray]:
        '''
        For L&R of face j+1/2:
        
        limiter = min_mod, van_leer, van_albada
        '''
        uUL = Uj   + 0.5*limiter(Ujp1-Uj,   Uj-Ujm1)
        uUR = Ujp1 - 0.5*limiter(Ujp2-Ujp1, Ujp1-Uj)
        return uUL, uUR
    
    @staticmethod
    def Upwind1_TVD_eigen(Ujm1: np.ndarray, Uj: np.ndarray, Ujp1: np.ndarray, 
            Ujp2: np.ndarray, limiter, roe_average, dx, dt) \
            -> Tuple[np.ndarray, np.ndarray]:
        '''
        For L&R of face j+1/2:
        
        limiter = min_mod, van_leer, van_albada
        '''
        rhoL = Uj  [0]; uL = Uj  [1]/rhoL; pL=(GAMMA-1)*(Uj  [2]-0.5*rhoL*uL**2)
        rhoR = Ujp1[0]; uR = Ujp1[1]/rhoR; pR=(GAMMA-1)*(Ujp1[2]-0.5*rhoR*uR**2)
        tHL  = (Uj  [2] + pL)/rhoL
        tHR  = (Ujp1[2] + pR)/rhoR
        
        roe_rho, roe_u, roe_H, roe_a \
            = roe_average(rhoL, rhoR, uL, uR, tHL, tHR)

        mL = np.zeros([3,3])
        mR = np.zeros([3,3])

        g1 = GAMMA-1
        u2 = roe_u**2
        a2 = roe_a**2
        ua = roe_u*roe_a
        
        mR[0, :] = np.array([1,           1,      1          ])
        mR[1, :] = np.array([roe_u-roe_a, roe_u,  roe_u+roe_a])
        mR[2, :] = np.array([roe_H-ua,    0.5*u2, roe_H+ua   ])
        
        mL[0, :] = np.array([0.5*g1*u2+ua, -g1*roe_u-roe_a,  g1]) / (2*a2)
        mL[1, :] = np.array([a2-0.5*g1*u2,  g1*roe_u,       -g1]) / a2
        mL[2, :] = np.array([0.5*g1*u2-ua, -g1*roe_u+roe_a,  g1]) / (2*a2)
        
        aa = np.dot(mL, Uj  -Ujm1)
        bb = np.dot(mL, Ujp1-Uj  )
        cc = np.dot(mL, Ujp2-Ujp1)
        
        dL = limiter(bb, aa)/dx
        dR = limiter(cc, bb)/dx
        
        eig= np.array([roe_u-roe_a, roe_u, roe_u+roe_a])

        wL = np.dot(mL, Uj  ) + 0.5*dL*(dx-dt*eig)
        wR = np.dot(mL, Ujp1) - 0.5*dR*(dx+dt*eig)

        uUL = np.dot(mR, wL)
        uUR = np.dot(mR, wR)

        return uUL, uUR

(3) MUSCL 格式

MUSCL 格式即 Monotonic Upstream-centered Scheme for Conservation Laws, (van Leer, 1979)。