5.2.2. 示例 2: 一维Euler方程组 (LU-SGS)

本示例演示双时间步 LU-SGS 隐式时间推进方法。

(1) 离散方法

  • 有限体积方法

  • 均匀结构网格

  • 半离散方法

  • 隐式时间推进 (LU-SGS, 二阶精度)

  • 一阶空间离散 (ROE 格式)

  • 无边界条件

一维 Euler 方程非定常流动的时间推进方程 (Eq.3.2.63) 可写作 Eq.3.2.76:

\[\begin{split}\left[ \left( \frac{1}{\Delta \tau} + \frac{1+\phi}{\Delta t} \right) \Delta x I + \delta_x A^{(m)} \right] \Delta \mathbf{U}^{(m)} = \text{RHS}^{(m)} \\ = \frac{\phi \Delta x}{\Delta t} ( \mathbf{U}^{n} - \mathbf{U}^{n-1} ) - \frac{(1+\phi) \Delta x}{\Delta t} ( \mathbf{U}^{(m)}-\mathbf{U}^{n} ) + \mathbf{R}(\mathbf{U}^{(m)})\end{split}\]

其中, \(\phi=0.5\) 时时间方向为二阶精度, \(\phi=0\) 时为一阶精度。 \(\sigma_A = (1+\varepsilon)|\lambda|_{\text{max}, A_i}\), \(\varepsilon \in [0,0.01]\)。 残差项 \(\mathbf{R}(\mathbf{U}^{(m)})\) 计算方法与显式时间推进中的处理方法相同 (式 Eq.3.2.38)。

L 块向前扫描运算: (Eq.3.2.87)

\[\Delta \mathbf{U}^*_i = \frac{ \text{RHS}^{(m)} + (A^{+} \Delta \mathbf{U})_{i-1}^{(m)} } { \frac{\Delta x}{\Delta \tau} + \frac{(1+\phi)\Delta x}{\Delta t} + \sigma_A }\]

U 块向后扫描运算: (Eq.3.2.89)

\[\Delta \mathbf{U}^{(m)}_i = \Delta \mathbf{U}^*_i - \frac{ (A^{-} \Delta \mathbf{U})_{i+1}^{(m)}} {\frac{\Delta x}{\Delta \tau} + \frac{(1+\phi)\Delta x}{\Delta t} + \sigma_A}\]

(2) 代码示例

以下为部分代码。 其中，调用的通量和重构格式参考 Eq.3.2.45, Eq.3.2.57 相关实现。

Tip

由于是二阶时间精度，相比于 example_01 物理时间步长必须减小。

使用结构体形式存储每个网格的数据:

 1class Cell():
 2    '''
 3    Data of each cell
 4    '''
 5    def __init__(self) -> None:
 6
 7        # conserved variables: [rho, rho u, rho E]
 8        # primitive variables: [rho, u, p]
 9
10        #* Physical step n-1
11        self.U_nm1 = np.zeros(3)    # conserved variables (physical step n-1)
12
13        #* Physical step n
14        self.U_n = np.zeros(3)      # conserved variables (physical step n)
15
16        #* pseudo step m
17        self.U_m = np.zeros(3)      # conserved variables (pseudo step m)
18        self.DU_m = np.zeros(3)     # difference of the conserved variables (pseudo step m)
19        self.DU_star = np.zeros(3)  # temporal DU* after L sweep (pseudo step m)
20        self.residual = np.zeros(3) # the right hand side of the dU/dt=Res
21
22        self.abs_eigenvalues = np.zeros(3)  # absolute eigenvalues [abs(u+a), abs(u), abs(u-a)]
23        self.lambda_max_m = 0.0             # maximum absolute eigenvalue of Jacobin A
24        self.mA_m = np.zeros([3,3])         # Jacobin A = dF/dU
25
26        self.pseudo_dt = 0.0        # pseudo time step
27        self.ratio_lu = 0.0         # ratio of (RHS + [A+]DU) or [A-]DU
28
29    @staticmethod
30    def JacobinA(u: float, tE: float) -> np.ndarray:
31        '''
32        >>> mA = Cell.JacobinA(U[1]/U[0], U[2]/U[0])
33        '''
34        u2 = u**2
35        gE = tE*GAMMA
36
37        mA = np.zeros([3,3])
38        mA[0,0] = 0
39        mA[0,1] = 1
40        mA[0,2] = 0
41        mA[1,0] = 0.5*(GAMMA-3)*u2
42        mA[1,1] = (3-GAMMA)*u
43        mA[1,2] = GAMMA - 1
44        mA[2,0] = (GAMMA-1)*u*u2 - gE*u
45        mA[2,1] = -1.5*(GAMMA-1)*u2 + gE
46        mA[2,2] = GAMMA*u
47
48        return mA
49
50    @staticmethod
51    def JacobinA_plus (mA: np.ndarray, lambda_max_m: float, ratio=1.01) -> np.ndarray:
52        return 0.5*(mA + np.eye(3)*lambda_max_m*ratio)
53
54    @staticmethod
55    def JacobinA_minus(mA: np.ndarray, lambda_max_m: float, ratio=1.01) -> np.ndarray:
56        return 0.5*(mA - np.eye(3)*lambda_max_m*ratio)

残差项 \(\mathbf{R}(\mathbf{U}^{(m)})\) 计算方法与显式时间推进中的处理方法相同:

 1def explicit_residual(Um2, Um1, U, Up1, Up2) -> np.ndarray:
 2    '''
 3    Calculate the right hand side of the (1/J) dU/dt = Res
 4    '''
 5    uUL, uUR = Reconstruction.Upwind1_TVD(Um2, Um1, U, Up1,
 6                limiter=Reconstruction.min_mod)
 7
 8    fFaceL = Roe.flux_face(uUL, uUR)
 9
10    uUL, uUR = Reconstruction.Upwind1_TVD(Um1, U, Up1, Up2,
11                limiter=Reconstruction.min_mod)
12
13    fFaceR = Roe.flux_face(uUL, uUR)
14
15    res = - (fFaceR - fFaceL)
16
17    return res

LU-SGS 的两次扫描过程:

  1def LU_SGS(physical_dt: float, dx: float, cfl: float, cells: List[Cell], n_pseudo_steps: int, phi: float):
  2    '''
  3    Update cell.U_n (need initial U_n & U_nm1)
  4    '''
  5    #* Initialization
  6    for i in range(len(cells)):
  7        cells[i].U_m = cells[i].U_n.copy()
  8
  9    #* Pseudo time iteration: L-U sweep
 10    for i_pseudo in range(n_pseudo_steps):
 11
 12        #* Preparation of all cells (m)
 13        for i in range(2, N_POINTS-1):
 14
 15            #* primitive variables
 16            rho = cells[i].U_m[0]
 17            u   = cells[i].U_m[1]/rho
 18            tE  = cells[i].U_m[2]/rho
 19            p   = (GAMMA-1)*(tE-0.5*rho*u**2)
 20            a   = np.sqrt(GAMMA*p/rho)
 21
 22            #* Calculate eigenvalues
 23            cells[i].abs_eigenvalues = np.array([abs(u+a), abs(u-a), abs(u)])
 24            cells[i].lambda_max_m = np.max(cells[i].abs_eigenvalues)
 25
 26            #* Calculate pseudo_dt
 27            cells[i].pseudo_dt = cfl*dx/cells[i].lambda_max_m
 28
 29            #* Calculate residual R(U(m))
 30            cells[i].residual = explicit_residual(
 31                cells[i-2].U_m, cells[i-1].U_m, cells[i].U_m, cells[i+1].U_m, cells[i+2].U_m)
 32
 33            #* Calculate [A]
 34            cells[i].mA_m = Cell.JacobinA(u, tE)
 35
 36            #* Calculate ratio of (RHS+[A+] DU) or [A-] DU
 37            denominator = dx/cells[i].pseudo_dt + (1+phi)*dx/physical_dt \
 38                        + RATIO_LAMBDA*cells[i].lambda_max_m
 39            cells[i].ratio_lu = 1.0/denominator
 40
 41
 42        #* L sweep
 43        for i in range(2, N_POINTS-1):
 44
 45            II = i-1
 46
 47            #* Calculate [A+]
 48            mAp = Cell.JacobinA_plus(cells[II].mA_m, cells[II].lambda_max_m, ratio=RATIO_LAMBDA)
 49
 50            #* Calculate RHS
 51            rhs  =    phi *dx/physical_dt*(cells[i].U_n-cells[i].U_nm1)
 52            rhs -= (1+phi)*dx/physical_dt*(cells[i].U_m-cells[i].U_n)
 53            rhs += cells[i].residual
 54
 55            #* Calculate dU*
 56            numerator = rhs + np.dot(mAp, cells[II].DU_m)
 57            cells[i].DU_star = cells[i].ratio_lu * numerator
 58
 59
 60        #* U sweep
 61        for i in range(N_POINTS-1, 2, -1):
 62
 63            II = i+1
 64
 65            #* Calculate [A-]
 66            mAm = Cell.JacobinA_minus(cells[II].mA_m, cells[II].lambda_max_m, ratio=RATIO_LAMBDA)
 67
 68            #* Calculate dU(m)
 69            numerator = np.dot(mAm, cells[II].DU_m)
 70            cells[i].DU_m = cells[i].DU_star + cells[i].ratio_lu * numerator
 71
 72
 73        #* Update conserved variables
 74        for i in range(2, N_POINTS-1):
 75            cells[i].U_m = cells[i].U_m + cells[i].DU_m
 76
 77
 78    #* Update conserved variables
 79    for i in range(2, N_POINTS-1):
 80        cells[i].U_nm1 = cells[i].U_n.copy()
 81        cells[i].U_n   = cells[i].U_m.copy()
 82
 83
 84    #* Residual
 85    residual = {}
 86    residual['physical-density'] = np.max([abs(cell.U_n[0] - cell.U_nm1[0]) for cell in cells[2:N_POINTS-1]])
 87    residual['pseudo-density']   = np.max([abs(cell.DU_m[0]) for cell in cells[2:N_POINTS-1]])
 88
 89    return residual
 90
 91if __name__ == "__main__":
 92
 93    phi = 0.0
 94    for i_physical in range(N_PHYSICAL_STEP):
 95
 96        if i_physical >= 2:
 97            phi = DS_PHI
 98
 99        residual = LU_SGS(physical_dt, DX, CFL, cells, N_PSEUDO_STEP, phi)
100
101    solution = np.zeros_like(initial_solution)
102    for i in range(mesh.shape[0]):
103        solution[i,:] = cells[i].U_n

(3) 结果展示

../../_images/example_021.jpg

图 5.2.2 ROE 格式隐式时间推进与精确解的对比