无粘通量-通量重构格式 MUSCL

无粘通量格式 MUSCL 基本思想、推导过程、基本形式等

Reconstruction

求解网格单元通量需要基于单元表面两侧物理量值计算通量差值,而通过已知离散点估计物理量空间分布的过程,称为重构 (Reconstruction)

  • 有限差分:通过网格节点上的值估计半点值
  • 有限体积(单元中心格式):通过控制体质点值估计控制体表面值
  • 有限元:通过单元函数估计单元表面值

最简单的通量重构方式:

  • 有限差分:半点值与网格节点值一致
  • 有限体积:控制体内物理量均匀分布,内表面物理量值等于体积分平均值

一阶精度重构

qi+1/2L=qi,qi+1/2R=qi+1\begin{equation} q_{i+1/2}^L = q_i, \quad q_{i+1/2}^R = q_{i+1} \end{equation}

MUSCL

monotone upstream-centered schemes for conservation laws 直译为 守恒定律的单调上游中心格式

MUSCL 二阶精度重构格式

qi+1/2L=qi+14[(1κ)Δ+(1+κ)Δ+]i,qi+1/2R=qi+114[(1κ)Δ++(1+κ)Δ]i+1\begin{gather} q_{i + 1/2}^L = q_i + \frac{1}{4} \begin{bmatrix} (1 - \kappa) \Delta_{-} + (1 + \kappa) \Delta_{+} \end{bmatrix}_{i}, \\ q_{i + 1/2}^R = q_{i+1} - \frac{1}{4} \begin{bmatrix} (1 - \kappa) \Delta_{+} + (1 + \kappa) \Delta_{-} \end{bmatrix}_{i+1} \end{gather}

式中 Δ\Delta 为变量梯度

(Δ+)i=qi+1qi,(Δ)i=qiqi1\begin{gather} (\Delta_{+})_{i} = q_{i + 1} - q_{i}, \quad (\Delta_{-})_{i} = q_{i} - q_{i-1} \end{gather}

展开上式可以得到

qi+1/2L=[κ14,1κ2,κ+14][qi1qiqi+1]qi+1/2R=[κ+14,1κ2,κ14][qiqi+1qi+2]\begin{align*} q_{i + 1/2}^L &= \begin{bmatrix} \displaystyle { \frac{\kappa - 1}{4}, 1 - \frac{\kappa}{2}, \frac{\kappa + 1}{4} } \end{bmatrix} \begin{bmatrix} q_{i - 1} \\ q_i \\ q_{i + 1} \end{bmatrix} \\ q_{i + 1/2}^R &= \begin{bmatrix} \displaystyle { \frac{\kappa + 1}{4}, 1 - \frac{\kappa}{2}, \frac{\kappa - 1}{4} } \end{bmatrix} \begin{bmatrix} q_{i} \\ q_{i + 1} \\ q_{i + 2} \end{bmatrix} \\ \end{align*}

  • κ=1\kappa = -1 为二阶迎风差分格式
  • κ=1\kappa = 1 为中心差分格式(需引入人工粘性)
  • κ=1/3\kappa = 1/3 为迎风偏置格式(1D为三阶,2D/3D为二阶)

表中 [,,,][*, *, *, *] 表示 [i1,i,i+1,i+2][i-1, i, i+1, i+2] 节点权重

κ\kappa qi+1/2Lq_{i + 1/2}^L qi+1/2Rq_{i + 1/2}^R
1-1 [1,3,0,0]:2[-1, 3, 0, 0] : 2 [0,0,3,1]:2[0, 0, 3, -1] : 2
00 [1,4,1,0]:4[-1, 4, 1, 0] : 4 [0,1,4,1]:4[0, 1, 4, -1] : 4
11 [0,1,1,0]:2[0, 1, 1, 0] : 2 [0,1,1,0]:2[0, 1, 1, 0] : 2
1/31/3 [1,5,2,0]:6[-1, 5, 2, 0] : 6 [0,2,5,1]:6[0, 2, 5, -1] : 6

Limiter

MUSCL 方法在流场光滑区结果较好,对高马赫数激波区,会由于重构值过大或过小引起上冲或下冲,为减小由此引起的振荡,需引入限制器限制重构梯度

qi+1/2L=qi+14[(1κ)Δˉ+(1+κ)Δˉ+]i,qi+1/2R=qi+114[(1κ)Δˉ++(1+κ)Δˉ]i+1\begin{gather} q_{i + 1/2}^L = q_i + \frac{1}{4} \begin{bmatrix} (1 - \kappa) \bar{\Delta}_{-} + (1 + \kappa) \bar{\Delta}_{+} \end{bmatrix}_{i}, \\ q_{i + 1/2}^R = q_{i+1} - \frac{1}{4} \begin{bmatrix} (1 - \kappa) \bar{\Delta}_{+} + (1 + \kappa) \bar{\Delta}_{-} \end{bmatrix}_{i+1} \end{gather}

限制器的对称性会影响求解结果

Limiter κ\kappa Symmetry
minmod
Van Albada \checkmark
Van Leer \checkmark
Hemker & Koren
Venkatakrishnan
Nishikawa
Barth

minmod

Δˉ=minmod(Δ,bΔ+),Δˉ+=minmod(Δ+,bΔ)\begin{equation} \bar{\Delta}_{-} = \text{minmod}(\Delta_{-}, b\Delta_{+}), \quad \bar{\Delta}_{+} = \text{minmod}(\Delta_{+}, b\Delta_{-}) \end{equation}

式中

1b3κ1κ\begin{equation} 1 \leq b \leq \frac{3 - \kappa}{1 - \kappa} \end{equation}

κ\kappa min\min max\max
-1 1 2
0 1 3
1/3 1 4
1 1 \infty

12(1+3κ1κ)=1+11κ\begin{equation*} \frac{1}{2} \cdot \left( 1 + \frac{3 - \kappa}{1 - \kappa} \right) = 1 + \frac{1}{1 - \kappa} \end{equation*}

minmod 限制器包含两个方面

  1. 抑制振荡 sign(x)+sign(y)=0ifxy<0\text{sign}(x) + \text{sign}(y) = 0 \quad \text{if} \quad x y < 0
  2. 限制梯度 min(x,y)min(|x|, |y|)

minmod(x,y)=12[sign(x)+sign(y)]min(x,y)\begin{equation} \text{minmod}(x, y) = \frac{1}{2}\begin{bmatrix} \text{sign}(x) + \text{sign}(y) \end{bmatrix} \min(|x|, |y|) \end{equation}

sign(x)={1x01x<0\begin{equation} \text{sign} (x) = \begin{cases} 1 \quad & x \ge 0 \\ -1 \quad & x < 0 \end{cases} \end{equation}

Van Leer

κ=1\kappa = -1 时限制器形式为

Δˉ=Δ+ˉ=vanLeer(Δ,Δ+)vanLeer(x,y)=[sign(x)+sign(y)]xyx+y+ε\begin{gather} \bar{\Delta_-} = \bar{\Delta_+} = \text{vanLeer} (\Delta_-, \Delta_+) \\ \text{vanLeer}(x, y) = \frac{\begin{bmatrix} \text{sign}(x) + \text{sign}(y) \end{bmatrix} xy}{|x| + |y| + \varepsilon} \end{gather}

其中 ε\varepsilon 是一个小量,可取 10610^{-6}

Van Albada

κ=1\kappa = -1 时限制器形式为

Δˉ=Δ+ˉ=vanAlbada(Δ,Δ+)vanAlbada(x,y)=x(y2+ε)+y(x2+ε)x2+y2+2ε\begin{gather} \bar{\Delta_-} = \bar{\Delta_+} = \text{vanAlbada} (\Delta_-, \Delta_+) \\ \text{vanAlbada} (x, y) = \frac{ x (y^2 + \varepsilon) + y (x^2 + \varepsilon) }{ x^2 + y^2 + 2 \varepsilon } \end{gather}

0%