一种模拟输移扩散问题的高精度差分方法

摘要

本发明涉及环境、水利等工程领域的数值模拟方法，适用于多类流体流动和传热传质问题、特别是输移扩散问题系列典型模型方程的高精度差分方法，通用性较强。本发明针对当前水环境模拟的常规有限差分格式只具有较低的一阶或二阶精度，而现有的多步差分等高精度方法一般涉及网格点较多，构造方式过于复杂，且多基于十分规则的均匀矩形网格等不足，提出了一种非均匀网格下的高精度差分方法。该方法采用待定系数法确定差分方程中待求点和相邻点之间的影响系数，使差分格式达到该网格结构下的最高阶精度。由于采用非均匀网格系统构造差分格式，网格间距的疏密可根据待求物理量变化梯度的大小灵活地调整，有效地提高了模拟精度和效率。

说明书

技术领域

本发明涉及环境、水利等工程领域的数值模拟方法，适用于多类流体流动和传热传质问 题、特别是水环境中输移扩散问题系列典型模型方程在非均匀网格下的高精度数值求解。

背景技术

在环境、水利、化工、能源等许多工程领域，流体的动量、热量和质量等基本参量在流 动过程中的输移扩散问题备受关注。各类流体流动和传热传质问题的根本解决，都必须基于 对流体环境中输移扩散规律的科学把握，以便了解流体或其中含有物质随时空的变化情况。 数值模拟方法是目前研究这类问题的常用手段。

综观国内外的环境水质模型研究可以发现，实际模拟计算中的常规数值方法(如迎风格 式、中心差分格式等)普遍精度较低，只具有一阶或二阶精度，使模拟结果往往不能准确真 实地反映客观实际。高精度方法能使计算过程中的收敛速度和结果精度明显提高，有效地提 高计算流体力学解决实际工程问题的能力，这种优势在计算大区域或复杂区域环境水流流动 问题时显得尤为必要。因此，高精度数值方法的研究可丰富和完善计算流体力学理论和方法， 研究成果具有向实际工程推广的应用价值。

然而，经检索发现，目前应用较为广泛的高精度有限差分方法以高阶多步方法为代表(Rai M.M.，1987；唐晨等，2004)，一般涉及网格点较多，构造方法较为复杂。此外，陈国谦等 (1994)提出的高精度摄动差分法仅涉及到相邻网格点，大大简化了差分格式的构造和求解 方法，但和多步方法一样，目前仅应用于十分规则的均匀矩形网格，在水环境数值模拟计算 中具有很大的局限性。实际计算中根据需要调整网格间距来生成非均匀网格，既能灵活地适 应各种不规则边界，有效地解决待求量(如浓度)突变等复杂问题，又能显著提高计算效率， 使高精度差分方法在解决实际工程问题时更加实用。随着现代科学技术的不断发展，如何简 便高效地获取非均匀网格下输移扩散问题诸多模型方程的高精度数值解，以满足水污染防治 对水环境精细模拟的迫切需要，已成为计算流体力学的一个重要研究领域。

发明内容

本发明的主要目的是提高水环境中输移扩散问题的模拟精度，克服常规数值方法精度较 低和高精度差分方法不便于实际应用的缺陷，提出输移扩散问题的实用精细数值模拟方法。

本发明针对当前水环境中输移扩散问题的数值模拟常采用的有限差分格式只具有一阶或 二阶精度，不能满足水环境精细模拟的要求；而现有的多步差分等高精度方法一般涉及网格 点较多，或计算较为复杂，且多基于十分规则的均匀矩形网格的现状，提出了一种非均匀网 格下的高精度差分方法。通过在紧致的网格系统中，设定待求点和相邻点之间的相互影响关 系，在使差分格式达到尽可能高的精度的前提下，求解代数方程组以确定影响系数，从而确 定最终的差分格式。该方法包括非均匀网格系统中，典型模型方程(扩散方程、对流扩散方 程、Laplace方程)的高精度差分方法，以及源项的高精度离散化差分方法。

与现有方法和技术相比，本方法具有以下优点：

(1)模拟精度高：常规的有限差分方法只具有一阶或二阶精度，本方法由于在确定的网 格形式下，以能达到该网格形式下的最高阶精度为前提，采用反推结点之间影响系数的方法 确定差分格式，一般可达到四到六阶精度。在实际应用中，甚至能取得比陈国谦等(1994)提 出的四阶摄动差分法更精确的模拟结果。

(2)求解效率高：本发明采用非均匀网格，网格间距可根据需要灵活调整，易于实现在 模拟的物理量变化梯度大的地方网格加密，反之网格变稀，故能有效地提高模拟效率。在实 际应用中，通过选取适当的网格间距，可取得比相同网格点数的均匀网格格式精度更高的模 拟结果。

具体实施方式

1.非均匀网格系统中，扩散方程的高精度差分方法

1.1无源项扩散方程的高精度差分方法

一维扩散问题的数学提法为：

$\frac{\partial T}{\partial t}=D\frac{{\partial }^{2}T}{{\partial x}^2}+f(x,t),(0\le x\le L,0<t<+\infty )---\left(1\right)$

T(0，t)＝g₁(t)，T(L，t)＝g₂(t)    (0＜t＜+∞)                    (2)

T(x，0)＝μ(x)    (0≤x≤L)                (3)

其中T(x，t)为温度、浓度等待求量；D为扩散系数(不失一般性，假设D为常数)；f(x，t)为 源函数；x，t为时空变量。当f(x，t)＝0和f(x，t)≠0时，分别称(1)式为无源扩散方程(或称 齐次方程)和含源扩散方程(或称非齐次方程)。本节先介绍f(x，t)＝0的情况，得出非均匀 网格上建立高精度差分格式的一般方法。

本发明方法首先在构造的局部网格内，近似设定相邻结点间的代数关系；根据各相邻结 点值在中心点展开的泰勒级数以及原微分方程，由最高阶相容性条件来推求该代数式中的待 定系数，以求达到最高精度；最后求解包含各结点值的代数方程组而得到数值解。

图1所示为一维方程非均匀网格布局结构示意图，假定待求量T在结点P(i，j)的值T_P与 其相邻结点值T_SW、T_SE、T_SC、T_W及T_E有关，考虑到网格结点待求量之间具有必然联系，构 造非均匀网格上无源扩散方程仅含两时层六结点的差分关系式如下。

T_P＝C_SCT_SC+C_SWT_SW+C_SET_SE+C_WT_W+C_ET_E                    (4)

其中C_SC、C_SW、C_SE、C_W、C_E为各结点上的待定差分系数。与在均匀网格上的情形不同， 这里C_SW≠C_SE，C_W≠C_E，即上、下游差分系数对中心点P(i，j)的贡献不同。

为使(4)式达到最高精度，利用原微分方程$\frac{\partial T}{\partial t}=D\frac{{\partial }^{2}T}{{\partial x}^2}$和泰勒级数反推各结点上的待 定差分系数，得：

$C_{\mathrm{SC}}=\frac{h_1^2+{3h}_1{h}_2+h_2^2-6\mathrm{D\tau }}{h_1^2+{3h}_1{h}_2+h_2^2+6\mathrm{D\tau }}---\left(5\right)$

$C_{\mathrm{SW}}=\frac{h_2(h_1^2+h_1{h}_2-h_2^2+6\mathrm{D\tau })}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+6\mathrm{D\tau })}---\left(6\right)$

$C_{\mathrm{SE}}=\frac{h_1(-h_1^2+h_1{h}_2+h_2^2+6\mathrm{D\tau })}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+6\mathrm{D\tau })}---\left(7\right)$

$C_W=\frac{-h_2(h_1^2+h_1{h}_2-h_2^2-6\mathrm{D\tau })}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+6\mathrm{D\tau })}---\left(8\right)$

$C_E=\frac{h_1(h_1^2-h_1{h}_2-h_2^2+6\mathrm{D\tau })}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+6\mathrm{D\tau })}---\left(9\right)$

其中h₁、h₂分别为点P(i，j)与上、下游相邻结点间的空间步长；τ为时间步长。

在此种实现非均匀网格的方式下，因h₁＝h₂，从而有C_SW＝C_SE，C_W＝C_E，故均匀网 格差分格式可视为非均匀网格差分格式的特例。

1.2源项的离散化方法

对于含源项的偏微分方程，当f(x，t)≠0时，对源项的处理方法直接影响着差分格式的 精度。本发明中，差分格式的构造可在相应无源方程格式的基础上增加一离散源项，利用原 微分方程，由相容性条件反演确定源项差分系数并保持原齐次方程格式的精度。非均匀网格 系统中，含源扩散方程的紧致差分格式可构造如下。

T_P＝C_SCT_SC+C_SWT_SW+C_SET_SE+C_WT_W+C_ET_E+C_ff_P                    (10)

式中C_f为源项的离散化系数，f_p为源函数f在点P(i，j)的值，其它符号意义同前。

为使(10)式达到最高精度，利用原微分方程$\frac{\partial T}{\partial t}=D\frac{{\partial }^{2}T}{{\partial x}^2}+f(x,t)$和泰勒级数反推各结点 上的待定差分系数，得C_SC、C_SW、C_SE、C_W、C_E表达式同(5)-(9)式，而

C_f＝C′_f+δ_f/f_P                    (11)

C′_f＝τ(C_SC+C_SW+C_SE)                (12)

${\delta }_f=-{\left(\frac{\partial f}{\partial x}\right)}_p·\tau ·(h_1{C}_{\mathrm{SW}}-h_2{C}_{\mathrm{SE}})+\frac{\tau }{2}(h_1^2{C}_{\mathrm{SW}}+h_2^2{C}_{\mathrm{SE}}){\left(\frac{{\partial }^{2}f}{{\partial x}^2}\right)}_p$

                 (13)

$-\frac{{\tau }^2}{2}(C_{\mathrm{SC}}+C_{\mathrm{SW}}+C_{\mathrm{SE}}){(\frac{\partial f}{\partial t}+D\frac{{\partial }^{2}f}{{\partial x}^2})}_p$

(10)式即为本发明提出的一维含源扩散方程(1)式的差分格式。本发明建立的差分格式(10) 考虑了源项f随时间的非均匀变化和在空间上的非线性分布，并通过(13)式中的δ_f对源项的 差分系数C_f进行修正，以确保差分格式的较高精度。当f(x，t)＝0时，(10)式中C_ff_p＝0，而 各结点上的差分系数仍由(5)-(9)式确定，这说明无论扩散方程(1)中是否包含源项，用本发明 方法构造的差分格式基本结构相同，差别只在于构造含源方程的差分格式时尚需对附加源项 进行处理，因此本发明方法对齐次或非齐次方程均具有普遍性，构造方法具有通用性，便于 推广应用。离散源项的这一高精度方法不仅适用于扩散方程，对其它形式的方程同样适用。

根据差分理论可以证明，当τ＝O(h²)时，(10)式具有三至四阶精度，可将其截断误差表 示为

$R_n=O(h^{3+\delta (h_1,h_2)},{\tau }^2)---\left(14\right)$

当h₁＝h₂时(即均匀网格)，δ(h₁，h₂)＝1；当h₁≠h₂时，0＜δ(h₁，h₂)＜1，格式的精度主 要取决于空间网格的分布情况。在实际流动计算中，可以根据数值经验来调整离散网格各结 点上下游的空间步长h₁和h₂，协调计算区域内最大误差和平均误差之间的关系，以获得比均 匀网格格式更精确的计算结果。

采用Fourier方法对(10)式的稳定性进行分析，可知格式具有无条件稳定性。

2.非均匀网格系统中，对流-扩散方程的高精度差分方法

包括在非均匀网格中，有效求解线性与非线性、定常与非定常对流-扩散方程，且精度 和稳定性俱佳的高精度差分方法。

一维对流-扩散问题的数学提法为：

$\frac{\partial T}{\partial t}+u(x,t)\frac{\partial T}{\partial x}=D\frac{{\partial }^{2}T}{{\partial x}^2}+f(x,t),(a<x<b,t>0)---\left(15\right)$

T(x，0)＝T₀(x)    (a≤x≤b)                        (16)

T(a，t)＝g₁(t)，T(b，t)＝g₂(t)    (t＞0)                (17)

其中T为待求量(如污染物浓度等)；u、f分别为对流流速和源汇函数；D为扩散系数，这 里视为常数；x，t为时空变量。

当T＝污染物浓度，u＝Const时，(15)式为线性对流-扩散方程；当u＝T，D＝运动 粘性系数ν时，(15)式为非线性对流-扩散(Burgers)方程，通常视为流体流动(Naviers-Stokes 方程)的模型方程。

2.1非均匀网格上线性方程的差分方法

考虑一维线性对流-扩散方程

$\frac{\partial T}{\partial t}+a\frac{\partial T}{\partial x}=D\frac{{\partial }^{2}T}{{\partial x}^2}+f(x,t)---\left(18\right)$

其中a为常数。

对流-扩散方程数值解法的一个基本困难是不对称对流算子的出现。为避免处理对流算 子，拟通过综合变换法将(18)式变换为含源扩散方程，来建立非均匀网格上(18)式的高精度差 分格式，具体如下。

设

代入方程(18)，将其变换为扩散方程，最后再通过反变换，得非均匀网格上对流-扩散方程(18) 式的差分格式：

C_PT_P＝C_SCT_SC+C_SWT_SW+C_SET_SE+C_WT_W+C_ET_E+C_ff_P+δ_f                (20)

其中

C_P＝C_W+C_E+C_SC+C_SW+C_SE                            (21)

$C_W=C_W^{\prime \prime }\psi \left(h_1\right)=C_W^{\prime \prime }\exp \left(\frac{{\mathrm{ah}}_1}{2D}\right)---\left(25\right)$

$C_E=C_E^{\prime \prime }\psi (-h_2)=C_E^{\prime \prime }\exp (-\frac{{\mathrm{ah}}_2}{2D})---\left(26\right)$

C_f＝τ(C″_SC+C″_SW+C″_SE)                        (27)

${\delta }_f=(C_{\mathrm{SC}}^{\prime \prime }+C_{\mathrm{SW}}^{\prime \prime }+C_{\mathrm{SE}}^{\prime \prime })(-\frac{a^2{\tau }^2}{4D}f-\frac{1}{2}{\tau }^2\frac{\partial f}{\partial t}+\frac{1}{2}{\mathrm{a\tau }}^2-\frac{\partial f}{\partial x}$

$-\frac{1}{2}{\mathrm{D\tau }}^2\frac{{\partial }^{2}f}{{\partial x}^2})+(h_1{C}_{\mathrm{SW}}^{\prime \prime }-h_2{C}_{\mathrm{SE}}^{\prime \prime })(\frac{\mathrm{a\tau }}{2D}f-\tau \frac{\partial f}{\partial x})+---\left(28\right)$

$(h_1^2{C}_{\mathrm{SW}}^{\prime \prime }+h_2^2{C}_{\mathrm{SE}}^{\prime \prime })(\frac{a^2\tau }{{8D}^2}f-\frac{\mathrm{a\tau }}{2D}\frac{\partial f}{\partial x}+\frac{1}{2}\tau \frac{{\partial }^{2}f}{{\partial x}^2})$

式中C″_SC、C″_SW、C″_SE、C″_W、C″_E分别为扩散方程差分格式(10)中的相应结点上的差分系数C_SC、 C_SW、C_SE、C_W、C_E。这样对流-扩散方程(18)式的高精度差分格式就可通过其等价扩散方 程式来间接获取，使较难处理的对流项不再出现。

易知，差分格式(20)保持了扩散方程差分格式(10)的三至四阶精度。当a＝0时，差分格 式(20)与含源扩散方程格式(10)完全相同。因此，非均匀网格上扩散方程的高精度差分格式可 视为对流-扩散方程格式的特例。

由(20)式可见，本发明格式充分体现了非定常对流-扩散问题的“历经”性：时间步长愈 大，上时段贡献(C_SW、C_SE和C_SC)愈小且逐渐趋于零。通过分析可以发现，当h₁＝h₂时(即均 匀网格情况)，格式的差分系数可随流速a的方向、大小的改变自动调整，上游点系数(a＞0时 为C_SW、C_W，反之为C_SE、C_E)恒大于下游点系数(a＞0时为C_SE、C_E，反之为C_SW、C_W)， 且随|a|h值增大，下游影响减小，终至可忽略不计，说明格式很好地反映了对流的“迎风”效应。 然而，当h₁≠h₂时(即在非均匀网格情况下)，C_W与C_E之间的关系问题变为一个关于h₁、h₂的非线性问题。当a、D和τ一定时，C_W与C_E的比值随h₁、h₂的取值变化而变化，格式在一 定范围内具有迎风性。

对(20)式进行稳定性分析，可知它同样具有无条件稳定性。

2.2非均匀网格上非线性方程的差分方法

考虑一维非线性非定常对流-扩散问题

$2A(x,t)\frac{\partial T}{\partial x}+\mathrm{Re}\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{{\partial x}^2}+f(x,t)---\left(29\right)$

其中A、f分别为对流系数和源汇函数；Re为雷诺数。

对于非定常问题，在某一确定时层，可视t为常数，而各待求量和函数值仅与x有关， 并将时间导数项视为定常问题的源项来近似地处理。

在(29)式中，令

$S(x,t)=f(x,t)-\mathrm{Re}\frac{\partial T}{\partial t}---\left(30\right)$

并进一步将原方程简化为定常问题(31)式来进行分析。

$2A\left(x\right)\frac{\partial T}{\partial x}=\frac{{\partial }^{2}T}{{\partial x}^2}+S\left(x\right)---\left(31\right)$

下面我们首先介绍非均匀网格上(31)式的高精度有限差分方法的建立。

2.2.1定常流动问题

令

T(x)＝θ(x)ψ(x)                    (32)

代入方程(31)，消除关于θ的对流因子，将(31)式变换为等价扩散方程，最后通过反变换得 非均匀网格上求解定常非线性对流-扩散方程(31)式的相应差分格式：

C_iT_i＝C_i-1T_i-1+C_i+1T_i+1-C_SS_i-δ_S        (33)

其中

C_i＝C_i-1+C_i+1                (34)

$C_{i-1}=C_{i-1}^{\prime }·\exp [A_i{h}_1-\frac{1}{2}\frac{\partial A}{\partial x}{h}_1^2+\frac{1}{3!}\frac{{\partial }^{2}A}{{\partial x}^2}{h}_1^3]---\left(35\right)$

$C_{i+1}=C_{i+1}^{\prime }·\exp [-A_i{h}_2-\frac{1}{2}\frac{\partial A}{\partial x}{h}_2^2-\frac{1}{3!}\frac{{\partial }^{2}A}{{\partial x}^2}{h}_2^3]---\left(36\right)$

$C_S=-\frac{C_{i-1}^{\prime }{h}_1^2}{2}-\frac{C_{i+1}^{\prime }{h}_2^2}{2}-\frac{1}{24}{C}_{i-1}^{\prime }{h}_1^{4}P-\frac{1}{24}{C}_{i+1}^{\prime }{h}_2^{4}P---\left(37\right)$

${\delta }_S=(\frac{\partial S}{\partial x}-\mathrm{SA})(\frac{C_{i-1}^{\prime }{h}_1^3}{6}-\frac{C_{i+1}^{\prime }{h}_2^3}{6})+{(A}^{2}S-\frac{\partial A}{\partial x}S-2A\frac{\partial S}{\partial x}$

                    (38)

$+\frac{{\partial }^{2}S}{{\partial x}^2})(-\frac{C_{i-1}^{\prime }{h}_1^4}{24}-\frac{C_{i+1}^{\prime }{h}_2^4}{24})$

其中

$C_{i-1}^{\prime }=\frac{-(h_2+\frac{1}{6}{h}_2^{3}P+\frac{1}{12}{h}_2^4\frac{\partial P}{\partial x})}{D_c}---\left(39\right)$

$C_{i+1}^{\prime }=\frac{-(h_1+\frac{1}{6}{h}_1^{3}P-\frac{1}{12}{h}_1^4\frac{\partial P}{\partial x})}{D_c}---\left(40\right)$

$D_c=-(h_2+\frac{1}{6}{h}_2^{3}P+\frac{1}{12}{h}_2^4\frac{\partial P}{\partial x})(1+\frac{1}{2}{h}_1^{2}P+\frac{1}{24}{h}_1^4{P}^2-\frac{1}{6}{h}_1^3\frac{\partial P}{\partial x})$

           (41)

$-(h_1+\frac{1}{6}{h}_1^{3}P-\frac{1}{12}{h}_1^4\frac{\partial P}{\partial x})(1+\frac{1}{2}{h}_2^{2}P+\frac{1}{24}{h}_2^4{P}^2+\frac{1}{6}{h}_2^3\frac{\partial P}{\partial x})$

在(33)式中，中心点系数等于相邻点系数之和，故相应系数矩阵对角占优，不会出现伪 物理振荡解。该差分格式具有迎风性，即随Re数的增大，下游影响减弱，终至可忽略不计。 另外，当A为常数时，上述推导中所有关于A的导数项为零，而原表达式仍然适用，说明本 发明方法对线性和非线性方程均具有普遍适应性。同样，本发明格式(33)的推导虽基于非均 匀网格，但均匀网格格式实质上是非均匀网格格式的特例，而后者在实际流场计算中更为实 用。理论分析可知，差分格式(33)式具有三至四阶精度，在均匀网格情况下取得四阶精度。

2.2.2非定常流动问题

数学表达式为(29)。利用2.2.1中的成果，并通过隐式格式来离散(30)式中的即

$\frac{\partial T}{\partial t}=\frac{T_i^{j+1}-T_i^j}{\tau }---\left(42\right)$

式中T_i^j和T_i^j+1的上标分别表示j和j+1时层。则(29)式的差分格式为

$(C_{i+1}+C_{i+1}-\frac{C_S\mathrm{Re}}{\tau })T_i^{j+1}=C_{i-1}{T}_{i-1}^{j+1}+C_{i+1}{T}_{i+1}^{j+1}-\frac{C_S\mathrm{Re}}{\tau }{T}_i^j-C_S{f}_i-{\delta }_S---\left(43\right)$

其中C_i-1、C_i+1分别由(35)、(36)式确定，C_S、δ_S分别由(37)和(38)式确定。

3.非均匀网格系统中，泊松方程的高精度差分方法

二维泊松方程的一般形式为

$\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2}=f(x,y),(x,y)\in D.---\left(44\right)$

其中T为待求量，x，y为空间坐标，f(s，y)为源分布函数。D是x-y平面上的有界区域，设 它的边界是分段光滑的曲线。当f(x，y)＝0时，(44)式称为拉普拉斯方程。

用两族平行直线对计算区域作矩形剖分，如图2所示为二维方程非均匀网格布局结构示 意图。构造非均匀网格上二维泊松方程包含九结点的差分关系式如下。

T_P＝C_NWT_NW+C_NCT_NC+C_NET_NE+C_WCT_WC+C_ECT_EC

                                                    (45)

+C_SWT_SW+C_SCT_SC+C_SET_SE+C_ff_P

其中C_NW、C_NC、C_NE、C_WC、C_EC、C_SW、C_SC、C_SE为各相邻结点上的差分系数，C_f为源项f(x，y) 的离散化系数，f_P为f(x，y)在中心点P(i，j)的值。

将(45)式中各结点值表示为点P(i，j)处的泰勒级数形式，注意在非均匀网格下中心点 P(i，j)在x和y方向距上下游结点的距离一般不相等，即h₁≠h₂，k₁≠k₂。反复利用原微分方程 (T_xx)_P＝(-T_yy)_P+f_P及其(T_xxxx)_P＝(T_yyyy)_P+(f_xx-f_yy)_P，(T_xxyy)_P＝(-T_yyyy)_P+(f_yy)_P等，得：

$C_{\mathrm{NW}}=\frac{h_2{k}_2(h_1^2+h_1{h}_2-h_2^2+k_1^2+k_1{k}_2-k_2^2)}{(h_1+h_2)(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(46\right)$

$C_{\mathrm{NC}}=\frac{k_2(h_1^2+3h_1{h}_2+h_2^2-k_1^2-k_1{k}_2+k_2^2)}{(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(47\right)$

$C_{\mathrm{NE}}=\frac{h_1{k}_2(-h_1^2+h_1{h}_2+h_2^2+k_1^2+k_1{k}_2-k_2^2)}{(h_1+h_2)(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(48\right)$

$C_{\mathrm{WC}}=\frac{h_2(-h_1^2-h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+{3k}_1{k}_2+k_2^2)}---\left(49\right)$

$C_{\mathrm{EC}}=\frac{h_1(h_1^2-h_1{h}_2-h_2^2+k_1^2+3k_1{k}_2+k_2^2)}{(h_1+h_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(50\right)$

$C_{\mathrm{SW}}=\frac{h_2{k}_1(h_1^2+h_1{h}_2-h_2^2-k_1^2+k_1{k}_2+k_2^2)}{\left(h_1+h_2\right)(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(51\right)$

$C_{\mathrm{SC}}=\frac{k_1(h_1^2+3h_1{h}_2+h_2^2+k_1^2-k_1{k}_2-k_2^2)}{(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(52\right)$

$C_{\mathrm{SE}}=\frac{h_1{k}_1(-h_1^2+h_1{h}_2+h_2^2-k_1^2+k_1{k}_2+k_2^2)}{(h_1+h_2)(k_1+k_2)(h_1^2+3h_1{h}_2+h_2^2+k_1^2+3k_1{k}_2+k_2^2)}---\left(53\right)$

$C_f=-\frac{F_P}{f_P}=-\frac{1}{2}[h_1^2(C_{\mathrm{NW}}+C_{\mathrm{SW}}+C_{\mathrm{WC}})+h_2^2(C_{\mathrm{EC}}+C_{\mathrm{NE}}+C_{\mathrm{SE}})]-\frac{{\delta }_f}{f_P}---\left(54\right)$

${\delta }_f=\frac{1}{6}[-h_1^3(C_{\mathrm{NW}}+C_{\mathrm{SW}}+C_{\mathrm{WC}})+h_2^3(C_{\mathrm{EC}}+C_{\mathrm{NE}}+C_{\mathrm{SE}})]\frac{\partial f}{\partial x}+\frac{1}{2}[h_1^2·$

$(k_1{C}_{\mathrm{NW}}-k_2{C}_{\mathrm{SW}})+h_2^2(k_1{C}_{\mathrm{NE}}-k_2{C}_{\mathrm{SE}})]\frac{\partial f}{\partial y}+\frac{1}{6}[h_1^3·(k_2{C}_{\mathrm{SW}}-$

$k_1{C}_{\mathrm{NW}})+h_2^3(k_1{C}_{\mathrm{NE}}-k_2{C}_{\mathrm{SE}})]\frac{{\partial }^{2}f}{\partial x\partial y}+\frac{1}{24}[h_1^4(C_{\mathrm{NW}}+C_{\mathrm{SW}}+C_{\mathrm{WC}})---\left(55\right)$

$+h_2^4(C_{\mathrm{EC}}+C_{\mathrm{NE}}+C_{\mathrm{SE}})](\frac{{\partial }^{2}f}{{\partial x}^2}-\frac{{\partial }^{2}f}{{\partial y}^2})+\frac{1}{4}[h_1^2(k_1^2{C}_{\mathrm{NW}}+k_2^2{C}_{\mathrm{SW}})$

$+h_2^2(k_1^2{C}_{\mathrm{NE}}+k_2^2{C}_{\mathrm{SE}})]\frac{{\partial }^{2}f}{{\partial y}^2}$

其中h₁、h₂和k₁、k₂分别为x、y方向上点P(i，j)距上、下游相邻结点的距离，如图2所示。

将(46)-(55)式代入(45)式，即可得到非均匀网格上数值求解二维泊松方程的差分格式。对 于一般非均匀网格而言，各结点上下游的网格间距h₁≠h₂，k₁≠k₂。由差分格式(45)式的 等价微分方程，可知其截断误差阶次不小于三阶；当h₁＝h₂＝h，k₁＝k₂＝k(即差分网格为 均匀网格)时，差分方程(45)的截断误差为R_n＝O(h⁴)或O(k⁴)，因此格式具有四阶精度；当 h₁＝h₂＝k₁＝k₂＝h，即采用均匀正方形差分网格时，格式可达六阶精度。

在实际流动计算中，需要根据所研究问题的特点和数值经验来调整各离散网格点上h₁、 h₂、k₁、k₂之间的关系，以获得最佳计算结果。

附图说明

1、附图1给出了本发明中一维方程高精度差分格式的非均匀网格布局结构示意图，包含两时 层六结点。其中P(i，j)为中心结点，SW、SE、SC、W和E为相邻结点。h₁、h₂为非均匀 网格的空间步长，τ为时间步长。

2、附图2给出了本发明中二维方程高精度差分格式的非均匀网格布局结构示意图，包含三时 层九结点。其中P(i，j)为中心结点，NW、NC、NE、WC、EC、SW、SC和SE为相邻结 点。h₁、h₂为非均匀网格下x方向的空间步长，k₁、k₂为y方向的空间步长。

实施例1

考虑如下无源项扩散方程的初边值问题

$\left(\begin{array}{cc}\frac{\partial T}{\partial t}=\frac{{\partial }^{2}T}{{\partial x}^2},& 0<x<1,t>0\\ T(x,0)=\sin \mathrm{\pi x},& 0\le x\le 1\\ T(0,t)=T(1,t)=0,& t>0\\ & \end{array}\right)$

该问题的解析解为$T(x,t)=e^{-{\pi }^{2}t}\sin \mathrm{\pi x}.$

为比较有据，选用一阶精度的显式格式、隐式格式作为参考格式，这些格式采用时间前差， 空间中差的格式，含两时层四结点，j时层已知，而j+1时层未知。由于显式格式为条件稳定， 稳定条件为$r=\frac{\mathrm{D\tau }}{\min {(h_1,h_2)}^2}\le \frac{1}{2},$故为保证收敛性须选用足够小的时间步长(计算中取τ＝ 10^-4)；隐式格式和本发明格式为无条件稳定，时间步长和空间步长的选取不受稳定条件的限 制。计算结果如表1所示。从表中可见，即使本发明格式在较粗的空间网格和时间网格 (τ＝0.04)下亦能取得比隐式格式在较细网格下(τ＝0.01)更为精确的结果；在相同的空间网格 下，即使本发明格式计算采用的时间步长为显式格式的100倍(二者分别为τ＝0.01和 τ＝10^-4)，仍然能取得精度高于后者的数值结果。因此，本发明方法明显优于传统差分方法， 在较粗的网格下亦能获得较高精度的数值解。

表1最大绝对误差变化过程比较

实施例2

考虑函数值可发生急剧变化的定常非线性问题(Burgers方程)

$u\frac{\partial u}{\partial x}=\frac{1}{\mathrm{Re}}\frac{{\partial }^{2}u}{{\partial x}^2}$

其解析解形式为

$u\left(x\right)=C_2\tanh (C_1-\frac{\mathrm{Re}}{2}{C}_{2}x)$

其中的系数C₁、C₂可由边界条件唯一确定。这里选定其解为

u(x)＝tanh[Re(1-2x)/4]

的情形，就区段(0≤x≤1)进行计算。为比较有据，选用中心差分格式或迎风差分格式(当 Reuh＜2时取中心差分格式，否则取迎风差分格式)为参考。其中非线性代数方程组以欠松 弛法迭代求解，迭代初值均从零值开始，而边界值给定。迭代收敛准则为$\max |u_i^{k+1}-u_i^k|\le ϵ,$其中k为迭代次数，ε取10^-6。

表2给出了Re数为20和200两种情况下的详细数值结果比较，其中非均匀网格点数N ＝10，结点分布详见表中。在同一非均匀网格系统中，本发明格式的计算结果均明显优于相 应参考格式，即使是在中间(x＝0.5左右)大梯度区域附近，本发明格式误差依然比参考格式 小几个数量级。说明本发明格式能较精确地模拟数值突变。

表2各网格点的误差比较

表3给出了在相同网格点数的均匀网格中(N＝20)，本发明格式与陈国谦等(1994)提出的 四阶精度摄动格式的比较。由表中可以看出，在相同的网格条件下，本发明格式能取得比四 阶摄动差分格式精度更高的计算结果，并同样在大Re数下有较高的分辨率，能较准确地反映 函数值的突然变化。

表3均匀网格下本发明格式与摄动格式数值解的比较

实施例3

考虑泊松方程

$\frac{{\partial }^{2}T}{{\partial x}^2}+\frac{{\partial }^{2}T}{{\partial y}^2}=-\frac{{\pi }^2}{4}(\frac{1}{L^2}+\frac{1}{H^2})\sin \left(\frac{\mathrm{\pi x}}{2L}\right)\sin \left(\frac{\mathrm{\pi y}}{2H}\right),$

0＜x＜2L，0＜y＜2H

T(x，0)＝T(x，2H)＝0，    0≤x≤2L

T(0，y)＝T(2L，y)＝0，    0≤y≤2H

其精确解为

$T=\sin \left(\frac{\mathrm{\pi x}}{2L}\right)\sin \left(\frac{\mathrm{\pi y}}{2H}\right)$

在同一非均匀网格下，运用本发明方法和五点差分方法进行计算，五点差分方法的网格 结构包含P、WC、EC、NC、SC五个点，计算结果如表4所示。从表中可以看出，对于含源 项情况，在相同的网格(9×9)下，本发明方法的计算精度仍明显高于五点差分方法。对于绝大 部分结点，本发明方法即使在粗网格(5×5)下亦能取得比细网格(9×9)下的五点方法更精确的计 算结果。随着网格点数的增加，本发明方法求解代数方程组的迭代次数逐渐变得比五点方法 小，说明本发明方法不仅能更真实地反映物理规律，而且能更快地逼近问题的真解。由此也 可以想见高精度差分方法在节省计算量，提高计算精度和效率方面的优越性和重要性，尤其 是对大型科学计算而言。

表4L/H＝0.5时各网格点的误差比较