西南大学学报 (自然科学版)  2018, Vol. 40 Issue (9): 91-95.  DOI: 10.13718/j.cnki.xdzk.2018.09.014
0
Article Options
  • PDF
  • Abstract
  • Figures
  • References
  • 扩展功能
    Email Alert
    RSS
    本文作者相关文章
    罗传胜
    李春光
    董建强
    景何仿
    欢迎关注西南大学期刊社
     

  • 求解对流扩散方程的一种高精度紧致差分格式    [PDF全文]
    罗传胜1, 李春光2, 董建强1, 景何仿3     
    1. 北方民族大学 数学与信息科学学院, 银川 750021;
    2. 北方民族大学 土木工程学院, 银川 750021;
    3. 北方民族大学 数值计算与工程应用研究所, 银川 750021
    摘要:在指数变换的基础上,将对流扩散方程变为扩散方程,消除了数值求解中较难处理的对流项,采用四阶紧致差分方法离散扩散方程的空间变量,采用扩展的$\frac{1}{3} $ -Simpson公式离散时间变量,格式的截断误差为Oτ4+h4).理论分析证明该格式是无条件稳定的.通过数值算例验证了本文方法的有效性.
    关键词高精度紧致格式    无条件稳定    指数变换    数值计算    扩展的$\frac{1}{3} $-Simpson公式    

    许多自然科学现象和工程领域均涉及对流传热问题.对于给定流场的对流传热现象,其模型方程为对流扩散方程.由于传热介质结构的复杂性和实际问题中边界条件的多样性,我们一般无法得到对流传热问题的解析解,因此研究高效率、高精度数值求解对流扩散方程的新算法具有重要的应用价值.

    对于对流扩散方程数值求解常见的差分格式有很多,如修正中心显示格式(CD)、迎风格式(FUD)、Crank-Nicolson格式(C-N)等,虽然都是绝对稳定的,但它们的截断误差均较低.目前已有许多研究对此进行了改进:文献[1-3]使用Hermite插值思路给出了求解空间的四阶差分格式,但只适用于稳态的对流扩散方程;文献[4-6]使用综合变换建立了求解对流扩散方程的一种两层四阶差分格式,但增加了离散方程节点数,计算比较复杂;文献[7]提出了求解一维定常对流扩散问题非均匀网格上的多项式型高阶紧致差分格式,但对于边界条件的处理较困难;文献[8]利用四阶精度的三次样条公式提出了时间二阶空间四阶精度的两层紧致隐格式;文献[9]提出了构造高阶精度的待定系数法,虽然很好地解决了对流项的耗散问题,但精度较低,且计算较复杂.

    本文利用指数变换将对流扩散方程变为扩散方程,然后采用三点四阶紧致差分方法离散扩散方程的空间变量,利用扩展的$ \frac{1}{3}$ -Simpson公式离散时间变量,构造出一种求解对流扩散方程的新的高精度紧致差分格式.

    1 差分格式的构造

    在计算流体力学中,天然河道中的水沙运动、水中污染物的扩散是最常见的一种流动现象,伴随着计算机技术的快速发展,建立数学模型已经成为模拟这些流动现象的一种重要手段,然而污染物的扩散和水沙运动数学模型可以写成统一的对流扩散方程形式,考虑如下简单形式的对流扩散方程

    $ \begin{array}{*{20}{c}} {\frac{{\partial u}}{{\partial t}} + p\frac{{\partial u}}{{\partial x}} = a\frac{{{\partial ^2}u}}{{{\partial ^2}x}}}&{\left( {x,t} \right) \in \mathit{\Omega } \times \left[ {0,T} \right]} \end{array} $ (1)

    的初边值问题

    $ \begin{array}{*{20}{c}} {u\left( {0,t} \right) = u\left( {1,t} \right) = 0}&t \end{array} \in \left[ {0,T} \right] $
    $ \begin{array}{*{20}{c}} {u\left( {x,0} \right) = d\left( x \right)}&{x \in \mathit{\Omega }} \end{array} $

    其中:u(xt)表示未知量,a为扩散系数,p为对流系数,Ω=[0, 1],d(x)是足够光滑的函数.

    先做指数变换

    $ u = {{\rm{e}}^{\frac{p}{{2a}}x - \frac{{p2}}{{4a}}t}}v $ (2)

    将(2)式带入(1)式后可得

    $ \frac{{\partial v}}{{\partial t}} = a\frac{{{\partial ^2}v}}{{{\partial ^2}x}} $ (3)

    先对计算区域进行离散,空间方向进行N等分网格划分,步长$h=\frac{1}{N} $,时间步长为τ,网格节点为(xjtn),其中xj=jhtn=j=0,1,2,…,N.对于任意固定的t,设vj(t)是v(xjt)的近似值,(vxx)j(t)是vxx(xjt)的近似值,(vt)j(t)是vt(xjt)的近似值.

    对空间内部节点采用四阶紧致差分公式[10]来离散空间变量

    $ {\left( {\frac{{{\partial ^2}v}}{{{\partial ^2}x}}} \right)_j} = {\left( {1 + \frac{{{h^2}}}{{12}}\delta _x^2} \right)^{ - 1}}\delta _x^2{v_j} + O\left( {{h^4}} \right) $ (4)

    其中δx2是二阶中心差分算子.将(4)式带入(3)式后可得到

    $ {\mathit{\boldsymbol{V}}_1}\left( t \right) = \mathit{\boldsymbol{CV}}\left( t \right) $ (5)

    其中:

    $ \mathit{\boldsymbol{C}} = {\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{B}} $
    $ \mathit{\boldsymbol{V}}\left( t \right) = {\left[ {\begin{array}{*{20}{c}} {v\left( {{x_1},t} \right)}\\ {v\left( {{x_2},t} \right)}\\ \vdots \\ {v\left( {{x_{N - 2}},t} \right)}\\ {v\left( {{x_{N - 1}},t} \right)} \end{array}} \right]_{\left( {N - 1} \right) \times 1}} $
    $ \mathit{\boldsymbol{A}} = {\left[ {\begin{array}{*{20}{c}} {\frac{5}{6}}&{\frac{1}{{12}}}&{}&{}\\ {\frac{1}{{12}}}&{\frac{5}{6}}& \ddots &{}\\ {}& \ddots&\ddots &{\frac{1}{{12}}}\\ {}&{}&{\frac{1}{{12}}}&{\frac{5}{6}} \end{array}} \right]_{\left( {N - 1} \right) \times \left( {N - 1} \right)}} $
    $ \mathit{\boldsymbol{B}} = \frac{a}{{{h^2}}}{\left[ {\begin{array}{*{20}{c}} { - 2}&1&{}&{}\\ 1&{ - 2}& \ddots &{}\\ {}& \ddots&\ddots &1\\ {}&{}&1&{ - 2} \end{array}} \right]_{\left( {N - 1} \right) \times \left( {N - 1} \right)}} $

    采用扩展的$ \frac{1}{3}$ -Simpson公式[11]对时间方向进行离散,设$\tilde v\left( t \right) $v(t)的近似值,在区间$[{t_{n - 1}}, {\rm{ }}{t_n}] $$\tilde v\left( t \right) $满足(5)式即

    $ {{\mathit{\boldsymbol{\tilde V}}}_t}\left( {{t_\alpha }} \right) = \mathit{\boldsymbol{C\tilde V}}\left( {{t_\alpha }} \right) $ (6)

    其中$\alpha = n - 1, \frac{{n - 1}}{2}, n $,则有

    $ \mathit{\boldsymbol{\tilde V}}\left( t \right) = {{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} + \tau D\left( s \right)\mathit{\boldsymbol{\tilde V}}_t^{n - 1} + \tau M\left( s \right)\mathit{\boldsymbol{\tilde V}}_t^{\frac{{n - 1}}{2}} + \tau H\left( s \right)\mathit{\boldsymbol{\tilde V}}_t^n $ (7)

    其中

    $ D\left( s \right) = s - \frac{3}{2}{s^2} + \frac{2}{3}{s^3},M\left( s \right) = 2{s^2} - \frac{4}{3}{s^3},H\left( s \right) = - \frac{1}{2}{s^2} + \frac{2}{3}{s^3},t = {t_{n - 1}} + \tau s,s \in \left[ {0,1} \right] $

    令(7)式中的s=1,$\mathit{\boldsymbol{\tilde V}}\left( t \right) $tn上满足(7)式,则有

    $ {{\mathit{\boldsymbol{\tilde V}}}^n} = {{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} + \frac{\tau }{6}\left( {\mathit{\boldsymbol{\tilde V}}_t^{n - 1} + 4\mathit{\boldsymbol{\tilde V}}_t^{n - \frac{1}{2}} + \mathit{\boldsymbol{\tilde V}}_t^n} \right) $ (8)

    令(7)式中的s=12,$\mathit{\boldsymbol{\tilde V}}\left( t \right) $${t_{\frac{{n - 1}}{2}}} $上满足(7)式,则

    $ {{\mathit{\boldsymbol{\tilde V}}}^{n - \frac{1}{2}}} = {{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} + \frac{\tau }{{24}}\left( {\mathit{\boldsymbol{5\tilde V}}_t^{n - 1} + 8\mathit{\boldsymbol{\tilde V}}_t^{n - \frac{1}{2}} - \mathit{\boldsymbol{\tilde V}}_t^n} \right) $ (9)

    将(6)式带入(8)式和(9)式可得到

    $ \left( {\mathit{\boldsymbol{I}} - \frac{\tau }{6}\mathit{\boldsymbol{C}}} \right){{\mathit{\boldsymbol{\tilde V}}}^n} = \left( {\mathit{\boldsymbol{I}} + \frac{\tau }{6}\mathit{\boldsymbol{C}}} \right){{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} + \frac{{2\tau }}{3}\mathit{\boldsymbol{C}}{{\mathit{\boldsymbol{\tilde V}}}^{n - \frac{1}{2}}} $ (10)
    $ \left( {\frac{\tau }{{24}}\mathit{\boldsymbol{C}}} \right){{\mathit{\boldsymbol{\tilde V}}}^n} = \left( {\mathit{\boldsymbol{I + }}\frac{{5\tau }}{{24}}\mathit{\boldsymbol{C}}} \right){{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} + \left( {\frac{\tau }{3}\mathit{\boldsymbol{C}} - \mathit{\boldsymbol{I}}} \right){{\mathit{\boldsymbol{\tilde V}}}^{n - \frac{1}{2}}} $ (11)

    其中I为单位矩阵.

    (10)式两边同时乘以$ \left( {\mathit{\boldsymbol{I}} - \frac{\tau }{3}\mathit{\boldsymbol{C}}} \right)$,(11)式两边同时乘以$\left( {\frac{{2\tau }}{3}\mathit{\boldsymbol{C}}} \right) $,两者相加可得:

    $ \left( {\mathit{\boldsymbol{I + }}\frac{\tau }{2}\mathit{\boldsymbol{C + }}\frac{{{\tau ^2}}}{{12}}{\mathit{\boldsymbol{C}}^2}} \right){{\mathit{\boldsymbol{\tilde V}}}^n} = \left( {\mathit{\boldsymbol{I}} - \frac{\tau }{2}\mathit{\boldsymbol{C + }}\frac{{{\tau ^2}}}{{12}}{\mathit{\boldsymbol{C}}^2}} \right){{\mathit{\boldsymbol{\tilde V}}}^{n - 1}} $ (12)

    通过文献[10]可知空间方向的离散具有四阶精度,通过文献[11]中的研究结果可知,时间方向也具有四阶精度,因此本文所构造的格式(12)的截断误差为O(τ4+h4).

    2 稳定性分析

    引理1   对于矩阵AB,假设λ$ \mathit{\boldsymbol{X}} \in {\mathbb{R}^{N - 1}}\left( {\mathit{\boldsymbol{X}} \ne 0} \right)$分别是矩阵A-1B的特征值和特征向量,则λ是实数且λ>0.

      由于λX分别是矩阵A-1B的特征值和特征向量,那么

    $ {\mathit{\boldsymbol{A}}^{ - 1}}\mathit{\boldsymbol{BX}} = \lambda \mathit{\boldsymbol{X}} \Rightarrow \lambda {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{AX}} = {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{BX}} $

    $ \mathit{\boldsymbol{X}} = {\left[ {{x_1},{x_2}, \cdots {x_{N - 1}}} \right]^{\rm{T}}} $

    那么

    $ {\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{BX}} = 12a\left( {x_1^2 + {{\left( {{x_1} - {x_2}} \right)}^2} + \cdots {{\left( {{X_{N - 2}} - {X_{N - 1}}} \right)}^2} + x_{N - 1}^2} \right) $

    因为a>0,所以XTBX>0,同理XTAX>0,由λXTAX=XTBX $ \Rightarrow $ λ>0且λ为实数.

    引理2   如果z>0,则有

    $ \left| {\frac{{12 - 6z + {z^2}}}{{12 + 6z + {z^2}}}} \right| < 1 $

      当z>0时,有

    $ \begin{array}{*{20}{c}} {288z + 24{z^3} > 0 \Rightarrow }\\ {144z + 12{z^3} > - 144z - 12{z^3} \Rightarrow }\\ {{z^4} + 60{z^2} + 144z + 12{z^3} > {z^4} + 60{z^2} - 144z - 12{z^3} \Rightarrow }\\ {{{\left( {12 - 6z + {z^2}} \right)}^2} < {{\left( {12 + 6z + {z^2}} \right)}^2}} \end{array} $

    由此,我们知

    $ \left| {\frac{{12 - 6z + {z^2}}}{{12 + 6z + {z^2}}}} \right| < 1 $

    定理1  格式(12)是无条件稳定的.

      假设λi(i=1,2,…)是矩阵C的任意特征值,则$(12-6\tau {{\lambda }_{i}}+{{\tau }^{2}}\lambda _{i}^{2}){{(12+6\tau {{\lambda }_{i}}+{{\tau }^{2}}\lambda _{i}^{2})}^{-1}} $$(12\mathit{\boldsymbol{I}}-6\tau \mathit{\boldsymbol{C}}+{{\tau }^{2}}{{\mathit{\boldsymbol{C}}}^{2}}){{(12\mathit{\boldsymbol{I}}+6\tau \mathit{\boldsymbol{C}}+{{\tau }^{2}}{{\mathit{\boldsymbol{C}}}^{2}})}^{-1}} $的特征值,根据引理1和引理2我们可以得矩阵C的特征值都为正的实数,且

    $ \left| {\left( {12 - 6\tau {\lambda _i} + {\tau ^2}\lambda _i^2} \right)\left( {12 + 6\tau {\lambda _i} + {\tau ^2}\lambda _i^2} \right)} \right| < 1 $

    因此,式(12)是无条件稳定的.

    3 数值算例

    考虑无内热源的对流传热控制方程模型模拟温度场变化的分布初边值问题

    $ \frac{{\partial T}}{{\partial t}} + 0.1\frac{{\partial T}}{{\partial x}} = 0.01\frac{{{\partial ^2}T}}{{\partial {x^2}}},\left( {x,t} \right) \in \left[ {0,1} \right] \times \left[ {0,1} \right] $

    该问题的精确解为

    $ T\left( {x,t} \right) = {{\rm{e}}^{5x - \left( {0.25 + 0.01{\pi ^2}} \right)t}}\sin \left( {\pi x} \right) $

    相应的初值条件为

    $ T\left( {x,0} \right) = {{\rm{e}}^{{\rm{5}}x}}\sin \left( {\pi x} \right),T\left( {0,t} \right) = T\left( {1,t} \right) = 0 $

    其中:边界条件x=0,x=1时温度场恒温;T(xt)表示温度场.通过精确解可以看出这是一个温度场随着空间位置和时间周期变化的对流扩散问题,对于上述算例,将区域均匀剖分,步长h=0.05,τ=0.001,在t=0.1时将本文格式和经典C-N,CD,FUD格式的计算结果进行比较(数值解与精确解的绝对误差比较见图 1表 1).

    图 1 t=0.1时数值格式模拟结果的比较
    表 1 t=0.1时绝对误差比较
    4 结论

    针对一维线性对流扩散方程提出了一种新的高精度紧致差分格式,其截断误差为O(τ4+h4),通过理论分析证明了该格式是无条件稳定的.通过数值算例将本文格式与FUD,CD,C-N格式的计算结果进行了比较,结果表明本文格式的计算结果精度较高.

    参考文献
    [1] 葛永斌, 田振夫, 詹咏, 等. 求解扩散方程的一种高精度紧致隐式差分方法[J]. 上海理工大学报, 2005, 27(2): 107-110.
    [2] 葛永斌, 田振夫, 吴文权. 含源项非定常对流扩散方程的高精度紧致隐式差分方法[J]. 水动力学研究与进展, 2006, 21(5): 619-625. DOI:10.3969/j.issn.1000-4874.2006.05.010
    [3] 田振夫. 一维对流扩散方程的四阶精度有限差分法[J]. 宁夏大学学报(自然科学版), 1995, 16(1): 30-35.
    [4] 杨志峰, 周雪漪, 许协庆. 对流-扩散方程的一种高精度优化差分格式[J]. 水动力学研究与进展(A辑), 1991, 6(1): 113-119.
    [5] 杨志峰, 陈国谦. 含源扩散方程的一种高精度差分方法[J]. 北京师范大学学报(自然科学版), 1992, 28(3): 315-316.
    [6] 王煊, 杨志峰. 基于非均匀网格求解非线性对流扩散问题的一种高精度差分格式[J]. 北京师范大学学报(自然科学版), 2003, 39(1): 131-137. DOI:10.3321/j.issn:0476-0301.2003.01.025
    [7] 田芳, 田振夫. 非均匀网格上求解对流扩散问题的高阶紧致差分方法[J]. 宁夏大学学报(自然科学版), 2009, 30(3): 209-212. DOI:10.3969/j.issn.0253-2328.2009.03.002
    [8] 林建国, 许维德, 陶尧森. 含源项非定常非线性对流扩散方程的三次样条四阶差分格式[J]. 水动力学研究与进展(A辑), 1994, 9(5): 599-602.
    [9] 张小峰, 张红武. Crank-Nicolson格式精度的改进[J]. 水科学进展, 2001, 12(1): 33-38. DOI:10.3321/j.issn:1001-6791.2001.01.006
    [10] DING Heng-fei, ZHANG Yun-xin. A New Difference Scheme with High Accuracy and Absolute Stability for Solving Convection-Diffusion Equations[J]. J Comput and Appl Math, 2009, 230(2): 600-606. DOI:10.1016/j.cam.2008.12.015
    [11] 马明书. 一维抛物型方程的一个新的高精度显示差分格式[J]. 数值计算与计算机应用, 2001, 22(2): 156-160. DOI:10.3969/j.issn.1000-3266.2001.02.010
    [12] 陆金甫, 关治. 偏微分方程数值解法[M]. 2版. 北京: 清华大学出版社, 2004.
    [13] 詹涌强, 谭志明. 求解抛物型方程的一个高精度隐格式[J]. 西南大学学报(自然科学版), 2013, 35(11): 81-85.
    A High-Order Compact Difference Scheme for Solving Convection-Diffusion Equations
    LUO Chuan-sheng1, LI Chun-guang2, DONG Jian-qiang1, JING He-fang3     
    1. School of Mathematics and Information Science, Beifang University of Nationalities, Yinchuan 750021, China;
    2. School of Civil Engineering, Beifang University of Nationalities, Yinchuan 750021, China;
    3. Institute of Numerical Calculation and Engineering Application, Beifang University of Nationalities, Yinchuan 750021, China
    Abstract: In this paper, based on exponential transform, the convection diffusion equation is transformed into a diffusion equation, thus eradicating the advection term, which is hard to treat in numerical solution. A high-order accurate implicit compact difference scheme is constructed for solving the one dimensional parabolic equation by the fourth-order pade' formula combined with time extension of the $\frac{1}{3} $-Simpson formulas. The truncation error of the scheme is O(τ4+h4). A theoretical analysis shows that the scheme is unconditionally stable. Numerical experiments verify the accuracy and reliability of the present scheme.
    Key words: high-order compact scheme    unconditionally stable    exponential transformation    numerical calculation    extension of the $\frac{1}{3} $-Simpson    
    X