5.1.2. 示例 2: 一维对流扩散方程

(1) 控制方程

考虑不可压缩流体 Eq.3.1.9, 写为无量纲形式：

(5.1.6)\[\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} - a \frac{\partial u}{\partial x}\]

其中，\(u\) 是被输运物质的浓度，流动速度 \(a=1.0\), 扩散系数 \(\alpha=0.01\), \(x \in [0,1]\), \(t \in [0,1]\)。

初始条件和边界条件为：

(5.1.7)\[\begin{split}& \left\{ \begin{array}{lr} u(x,0) = \exp \left( - \frac{(x+0.5)^2}{0.00125} \right) \\ u(0,t) = \frac{0.025}{\sqrt{0.000625+0.02t}} \exp \left( - \frac{(0.5-t)^2}{0.00125+0.04t} \right) \\ u(1,t) = \frac{0.025}{\sqrt{0.000625+0.02t}} \exp \left( - \frac{(1.5-t)^2}{0.00125+0.04t} \right) \\ \end{array} \right.\end{split}\]

精确解为：

(5.1.8)\[u(x,t) = \frac{0.025}{\sqrt{0.000625+0.02t}} \exp \left( - \frac{(x+0.5-t)^2}{0.00125+0.04t} \right)\]

(2) 离散方法

有限差分方法
均匀结构网格
全离散方法
显式/隐式时间推进
一阶空间离散

记 \(\gamma\) 为空间加权系数，\(\phi\) 为时间加权系数。

时间导数的一阶差分：

(5.1.9)\[\frac{\partial u}{\partial t} \approx \frac{\hat u_i- u_i}{\Delta t}\]

空间一阶导数的迎风格式：

(5.1.10)\[\begin{split}\frac{\partial u}{\partial x} \approx (1-\phi) \left[\frac{(1-\gamma)(u_i- u_{i-1}) +\gamma(u_{i+1}-u_i)}{\Delta x}\right] \\ + \phi \left[\frac{(1-\gamma)(\hat u_i- \hat u_{i-1}) +\gamma(\hat u_{i+1}- \hat u_i)}{\Delta x}\right]\end{split}\]

空间二阶导数的中心差分：

(5.1.11)\[\frac{\partial^2 u}{\partial x^2} \approx (1-\phi) \left[ \frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2} \right] + \phi \left[ \frac{\hat u_{i+1}-2 \hat u_i+ \hat u_{i-1}}{\Delta x^2} \right]\]

将 Eq.5.1.10 和 Eq.5.1.11 代入 Eq.5.1.6, 则显式/隐式格式统一表达为：

(5.1.12)\[\hat u_i = \frac{1}{A_0} (A_1 u_{i-1} + A_2 u_{i} + A_3 u_{i+1} + A_4 \hat u_{i-1} + A_5 \hat u_{i+1})\]

(5.1.13)\[\begin{split}A_0 &= 1 - \phi [c(2\gamma-1)-2s], & A_1 = (\phi-1)[c(\gamma-1)-s],\\ A_2 &= 1 + (\phi-1)[c(1-2\gamma)+2s], & A_3 = (1-\phi)[s-c\gamma], \\ A_4 &= \phi [s+c(1-\gamma)], & A_5 = \phi[s-\gamma c]\end{split}\]

其中，\(c=a \Delta t/\Delta x\), \(s=\alpha \Delta t/\Delta x^2\)。

(2.a) Lax-Wendroff 格式

取 \(phi=0\), \(\gamma=(1-c)/2\), Eq.5.1.12 得：

(5.1.14)\[\hat u_i = \frac{1}{2}(2s+c+c^2)u_{i-1} + (1-2s-c^2)u_i + \frac{1}{2}(2s-c+c^2)u_{i+1}\]

本格式为接近迎风格式的显式时间推进格式，数值稳定区间为 \(0 < s ≤ (1-c^2)/2\)。

(2.b) Crank-Nicolson 格式

取 \(phi=1/2\), \(\gamma=1/2\), Eq.5.1.12 得：

(5.1.15)\[\begin{split}- (c+2s) \hat u_{i-1} + (4+4s) \hat u_i + (c-2s) \hat u_{i+1} = \\ (c+2s) u_{i-1} + (4-4s) u_{i} - (c-2s) u_{i+1}\end{split}\]

本格式为无条件稳定隐式时间推进格式，可以写为线性方程组形式 \(A\mathbf{u}=\mathbf{b}\)。其中 \(A\) 是一个近似三对角形式的方阵， \(\mathbf{u}\) 是下一时间步的未知解向量， \(\mathbf{b}\) 是由当前时间步的解计算得到的向量。

边界条件的处理：

(5.1.16)\[\begin{split}& A[0,0]=1, A[0,1]=0, \\ & A[n,n]=1, A[n,n-1]=0, \\ & b[0]=u(0,t), b[n]=u(1,t), n = \text{rank}(b)\end{split}\]

(3) 代码示例

以下为部分代码。

def lax_wendroff(mesh: np.array):
    
    print('Lax-Wendroff')
    
    def step(un: np.array, t: float):
        
        c1 = 0.5*(2*SS+CC+CC**2)
        c2 = 1-2*SS-CC**2
        c3 = 0.5*(2*SS-CC+CC**2)
        uu = np.zeros_like(un)
        
        uu[0] = boundary_condition_0(t)
        uu[-1]= boundary_condition_1(t)
        
        for i in range(1, N_POINTS-1):
            uu[i] = c1*un[i-1] + c2*un[i] +c3*un[i+1]

        return uu
    
    current_solution = initial_condition(mesh)
    
    for iter in range(N_TIME_STEPS):
        
        t = (iter+1)*DT
        current_solution = step(current_solution, t)
    
    return current_solution

def crank_nicolson(mesh: np.ndarray):

    def ndarray_b(un: np.ndarray):
        
        bb = np.zeros_like(un)

        _c1 = CC+2*SS
        _c2 = 4-4*SS
        _c3 = -(CC-2*SS)
    
        bb[0] = un[0]
        bb[-1]= un[-1]
    
        for i in range(1, N_POINTS-1):
            bb[i] = _c1*un[i-1] + _c2*un[i] + _c3*un[i+1]
        
        return bb

    c1 = -(CC+2*SS)
    c2 = 4+4*SS
    c3 = CC-2*SS

    # Create a new sparse PETSc matrix, fill it and then assemble it
    A = PETSc.Mat().createAIJ([N_POINTS, N_POINTS])
    A.setUp()
    
    A.setValue(0, 0, 1.0)
    A.setValue(N_POINTS-1, N_POINTS-1, 1.0)
    
    for i in range(1, N_POINTS-1):
        A.setValue(i, i-1, c1)
        A.setValue(i, i,   c2)
        A.setValue(i, i+1, c3)
    A.assemble()

    # Assemble the initial rhs to the linear system
    b = PETSc.Vec().createSeq(N_POINTS)
    b.setArray(ndarray_b(initial_condition(mesh)))

    # Allocate a PETSc vector storing the solution to the linear system
    x = PETSc.Vec().createSeq(N_POINTS)

    # Instantiate a linear solver: Krylow subspace linear iterative solver
    ksp = PETSc.KSP().create()
    ksp.setOperators(A)
    ksp.setFromOptions()
    
    chosen_solver = ksp.getType()
    print('Crank-Nicolson [solver: %s]'%(chosen_solver))
    
    for iter in range(N_TIME_STEPS):
        
        t = (iter+1)*DT
        ksp.solve(b, x)

        # Re-assemble the rhs to move forward in time
        current_solution = x.getArray()
        current_solution[0] = boundary_condition_0(t)
        current_solution[-1]= boundary_condition_1(t)

        b.setArray(ndarray_b(current_solution))
    
    return current_solution

(4) 结果展示

../../_images/example_02.jpg — 图 5.1.2 Lax-Wendroff 和 Crank-Nicolson 与精确解的对比