曹富軍，袁冬芳，葛永斌

(1.中國科技大學數(shù)學科學學院，安徽合肥230000;2.內(nèi)蒙古科技大學數(shù)理與生物工程學院，內(nèi)蒙古包頭014010;3.北京應用物理與計算數(shù)學研究所計算物理實驗室，北京100088;4.寧夏大學數(shù)學計算機學院，寧夏銀川750021)

0 引言

考慮具有如下形式的2維對流擴散方程:

其中u(x，y)是未知量，c(x，y)和 d(x，y)分別為 x和y方向的對流系數(shù)，f(x，y)為計算域上足夠光滑的函數(shù).

求解對流擴散方程的方法較多，如傳統(tǒng)的中心差分格式、迎風格式等，然而，它們都不能夠滿足實際計算所需要的精度.緊致差分格式由于精度高并且具有小的離散子域、對單元的敏感性小和對邊界單元的處理無特殊困難等特點越來越受到人們的重視［1-7］.然而，大多數(shù)緊致格式是基于均勻網(wǎng)格上提出的，對于求解大梯度或者邊界層問題，均勻網(wǎng)格上的緊致格式很難得到理論上的高精度數(shù)值解.因此，J.C.Kalita等［1］直接在物理區(qū)域上建立了求解2維定常對流擴散方程非均勻網(wǎng)格的高精度緊致差分格式，由于該方法在物理區(qū)域上直接離散偏微分方程，沒有采用坐標變換，因此減少了坐標變換帶來的額外計算量，并且數(shù)值結果表明該方法具有4階精度.

另一方面，多重網(wǎng)格方法已經(jīng)從理論上被證明至少對于求解線性橢圓型問題是一種最優(yōu)化的數(shù)值計算方法［2，8-11］.許多學者致力于高精度緊致差分格式與多重網(wǎng)格方法相結合求解偏微分方程并做出了不少有意義的工作.文獻［2］提出了非等距網(wǎng)格上2維泊松方程的網(wǎng)格部分半粗化策略(網(wǎng)格粗化過程僅在主控方向上進行，而在非主控方向上的網(wǎng)格步長則保持不變)，數(shù)值實驗表明:對于各向異性問題，由于減少了非主控方向的網(wǎng)格數(shù)，該策略下的多重網(wǎng)格算法比網(wǎng)格完全粗化策略下的多重網(wǎng)格算法，在不損失精度的前提下，大大節(jié)約了計算時間，提高了計算效率.文獻［9］將部分半粗化策略推廣到非等距網(wǎng)格上3維泊松方程的求解;文獻［10］基于J.C.Kalita等提出的求解2維對流擴散方程的高精度緊致格式，提出相應的完全粗化的多重網(wǎng)格算法，并通過對邊界層問題的求解，展現(xiàn)了非均勻網(wǎng)格上多重網(wǎng)格算法的效率.文獻［11］將該方法推廣到3維，建立了3維非均勻網(wǎng)格上完全粗化多重網(wǎng)格算法的插值算子和限制算子，并求解了3維泊松方程.該方法的缺點是在所有坐標軸方向上采用相同的網(wǎng)格數(shù)，這對于邊界層或大梯度僅在1個方向的問題，無疑會浪費計算量，因為在該方向并不需要網(wǎng)格加密.基于此，本文將設計一種非均勻網(wǎng)格上的網(wǎng)格部分半粗化策略，并建立其相應的多重網(wǎng)格算法，然后通過數(shù)值算例驗證方法的精確性和可靠性.

1 HOC格式

考慮矩形區(qū)域Ω=［a1，a2］×［b1，b2］，將求解區(qū)域［a1，a2］和［b1，b2］分別剖分a1=x0＜ x1＜x2＜…＜xm-1＜xm=a2，b1=y0＜y1＜y2＜…＜yn-1＜ yn=b2.并且定義:xf=xi+1-xi，xb=xi-xi-1，yf=yi+1-yi，yb=yi-yi-1，1≤i≤m-1，1≤j≤n-1.記 xi-xi-1= θlxhx，xi+1-xi= θrxhx，yj-yj-1=θlyhy，yj+1-yj= θryhy，其中 hx=(a2-a1)/n，hy=(b2-b1)/m.進一步地，記 β∧= θr∧+ θl∧，γ∧=θr∧-θl∧，α∧= θl∧θr∧，∧代表 x 或 y. 當且僅當θl∧=θr∧=1時網(wǎng)格剖分為均勻剖分.

通過以上的標記，x和y方向上的1階和2階差分算子可以寫為

2 部分半粗化多重網(wǎng)格方法

文獻［2］采用非等距網(wǎng)格離散2維泊松方程并提出2種多重網(wǎng)格方法:(i)在主方向上使用線松弛方法代替點松弛迭代，因為主控方向上節(jié)點多計算量大，線松弛可以快速消除高頻誤差;(ii)采用部分半粗化的多重網(wǎng)格方法.由于HOC格式在不同方向上網(wǎng)格數(shù)也不一樣，所以粗化過程僅在主控方向上進行實施，而非主控方向上網(wǎng)格數(shù)保持不變.數(shù)值實驗表明該方法對于求解各向異性問題表現(xiàn)出極高的效率.文獻［11］將該方法推廣到3維Poisson方程各向異性問題的求解.本文利用文獻［2］提出的部分半粗化多重網(wǎng)格方法的思想，結合非均勻網(wǎng)格上的HOC格式來求解帶邊界層的2維對流擴散問題.對于傳統(tǒng)的多重網(wǎng)格方法，限制算子和插值算子都是基于完全粗化過程，在均勻網(wǎng)格上提出的.然而，當采用非均勻網(wǎng)格時，傳統(tǒng)基于均勻網(wǎng)格上的限制算子和插值算子不再使用，因此，文獻［11］構造了非均勻網(wǎng)格上的限制算子和插值算子，并建立了相應的多重網(wǎng)格算法.但是該方法只適用于所有空間方向網(wǎng)格數(shù)相同的情形，如果不同空間方向網(wǎng)格數(shù)不同，則需要重新構造限制算子和插值算子，并建立部分半粗化的多重網(wǎng)格算法.下面先建立部分半粗化網(wǎng)格策略下多重網(wǎng)格方法的限制算子和插值算子.

首先，假設Nx和Ny分別為x和y方向的網(wǎng)格數(shù)且Nx＞Ny.對于部分半粗化過程而言，需要在細網(wǎng)格上在x方向從第1條邊開始每隔1個刪除1條線進行粗化.由于在進行x方向的粗化過程中y方向的網(wǎng)格數(shù)不變.因此，當Nx=Ny時，采用完全粗化過程.整個過程為部分半粗化.標準的完全粗化和部分半粗化過程分別如圖1和圖2所示.

為了便于表述，記uij為方程的近似解，Ωh和ΩH分別為細網(wǎng)格和粗網(wǎng)格.為細網(wǎng)格到粗網(wǎng)格上的限制算子為從粗網(wǎng)格到細網(wǎng)格上的插值算子和rij分別為粗網(wǎng)格和細網(wǎng)格上對應的殘差.由于x和y方向的網(wǎng)格數(shù)不同，需要針對不同情況構造限制算子.

當Nx≠Ny(Nx＞Ny)時，只需要在x方向根據(jù)相鄰點進行加權平均就可構造限制算子(見圖3):

圖1 2個方向完全粗化

圖2 x方向部分半粗化

當Nx≠Ny時，只需要在x方向上構造插值算子(見圖4):

其中較大的黑點表示粗網(wǎng)格上的網(wǎng)格點，較小的黑點表示細網(wǎng)格上的網(wǎng)格點，且表示 2i表示2j.

圖3 單個方向的限制算子

圖4 單個方向的插值算子

當Nx=Ny時，需要同時考慮x和y方向的8個網(wǎng)格點進行加權平均構造完全粗化限制算子(見圖5):

圖5 非均勻網(wǎng)格限制算子面積率系數(shù)示例

當Nx=Ny時，考慮x和y方向的點構造完全粗化插值算子(見圖6):

圖6 非均勻網(wǎng)格插值算子面積率系數(shù)示例

3 數(shù)值算例

例1 考慮2維對流擴散方程［11］:

其中(x，y)∈［0，1］× ［0，1］，c(x，y)=Rex(x-1)(1-2y)，d(x，y)=Rey(y-1)(1-2x).該方程的精確解為 u(x，y)=e-σ(x-0.5)2-y2.

該問題在x=0.5處有陡峭的邊界層，因此，考慮x方向為主方向，采用如下的伸縮函數(shù)生成非均勻網(wǎng)格:

其中λx(-1＜λx＜1)為伸縮系數(shù)，控制x方向網(wǎng)格的分布.比如，當 -1＜λx＜0時，網(wǎng)格在x=0和x=1附近分布密集，λx越小，x=0和x=1附近的網(wǎng)格點越密;當0＜λx＜1時，網(wǎng)格在x=0.5附近分布密集，且λx越大，x=0.5附近的網(wǎng)格點越密.當λx=0時，得到物理區(qū)域上的均勻網(wǎng)格.

表1給出了當σ=1×103時均勻網(wǎng)格和非均勻網(wǎng)格上的最大絕對誤差、迭代次數(shù)和CPU時間比較.從表1中可以看出，非均勻網(wǎng)格上在從方向上網(wǎng)格數(shù)較少的情況下得到的結果比均勻網(wǎng)格從方向上網(wǎng)格大時還精確.比如，網(wǎng)格數(shù)為256×32的非均勻網(wǎng)格上計算結果比網(wǎng)格數(shù)為256×256的均勻網(wǎng)格上計算結果更精確.同時，在非均勻網(wǎng)格上計算迭代步數(shù)和CPU時間都基本與均勻網(wǎng)格上的保持一致.因此，在非均勻網(wǎng)格上采用部分半粗化方法能夠在不增加計算消耗的情況下提高問題求解的精度.

表1 均勻和非均勻網(wǎng)格、迭代次數(shù)和CPU時間比較(Re=1000，σ=1×103)

表2給出了當σ=1×104時非均勻網(wǎng)格和均勻網(wǎng)格上計算誤差、迭代步和CPU時間比較.從2表中可以看出，非均勻網(wǎng)格上的計算結果精度遠遠高于均勻網(wǎng)格.對于x方向的網(wǎng)格數(shù)64和32來說，由于邊界層的存在，結果在邊界層附近變得非常陡峭，均勻網(wǎng)格上的計算結果精度較差，然而非均勻網(wǎng)格上仍然能夠得到較精確的結果.從結果來看，非均勻網(wǎng)格得到比均勻網(wǎng)格上更精確的結果且所需CPU時間僅為均勻網(wǎng)格上的1/8.

為了展示部分半粗化方法在整個區(qū)域上的求解精度和效率，圖7給出當網(wǎng)格數(shù)為64×16，σ=1×105時，方程的精確解(見圖7(b))，均勻網(wǎng)格上的數(shù)值解(見圖7(c))及非均勻網(wǎng)格上數(shù)值解(見圖7(d))的比較.從圖7(a)中可以看出非均勻網(wǎng)格在邊界層附近分布了足夠的網(wǎng)格，因此，其計算結果與精確解吻合的很好.然而，均勻網(wǎng)格上由于在邊界層附近的網(wǎng)格數(shù)不足，數(shù)值結果與精確解比較，表現(xiàn)出了很大的計算誤差.

表2 均勻和非均勻網(wǎng)格最大誤差、迭代次數(shù)和CPU時間比較(Re=1000，σ=1×104)

圖7 當σ=1×105，Re=1000時部分半粗化方法的精確解和數(shù)值解

4 結論

本文基于面積率公式構造了部分半粗化多重網(wǎng)格實施過程的插值算子和限制算子，并結合非均勻網(wǎng)格上的HOC格式與部分半粗化的多重網(wǎng)格方法對邊界層2維對流擴散方程進行求解.部分半粗化的多重網(wǎng)格方法根據(jù)問題邊界層所在的方向，允許不同方向上分布不同數(shù)目的網(wǎng)格，由于減少了無邊界層方向的網(wǎng)格數(shù)，因此節(jié)約了計算量和計算時間.數(shù)值實驗表明非均勻網(wǎng)格上的HOC格式和部分半粗化多重網(wǎng)格方法的結合對于帶邊界層的對流擴散方程的求解在精度和效率上表現(xiàn)出更大的優(yōu)勢.

［1］ Kalita J C，DassA K，Dalal D C.A transformation-free HOC scheme for steady convection-diffusion on non-uniform grids［J］.Int J Numer Meth Fluids，2004，44(1):33-53.

［2］ Zhang Jun.Multigridmethod and fourth-order compact scheme for 2D poisson equation with unequal grids size discretization ［J］.J Comput Phys，2002，179(1):170-179.

［3］Tian Zhengfu，Dai Shiqiang.High-order compact exponential finite differencemethods for convection-diffusion type problems［J］.J Comput Phys，2007，220(2):952-974.

［4］曹瑩，孔令華，王蘭，等.3維Helmholtz方程的4階緊致有限差分格式［J］.江西師范大學學報:自然科學版，2010，34(6):597-599.

［5］黃雪芳，郭銳，葛永斌.一維非定常對流擴散方程非均勻網(wǎng)格上的高精度緊靠致差分格式［J］.工程數(shù)學學報，2014(3):371-380.

［6］袁冬芳，曹富軍，石琳.對流擴散方程非均勻網(wǎng)格上四階緊致格式的預條件迭代法［J］.山西大學學報:自然科學版，2013(3):379-388.

［7］Wang Jie，Zhong Weijun，Zhang Jun.A generalmeshsize fourth-order compact difference discretization scheme for 3D poisson equation ［J］.Appl Math Comput，2006，183(2):804-812.

［8］Ge Lixin，Zhang Jun.Accuracy，robustness and efficiency comparison in iterative computation of convection diffusion equation with boundary layers［J］.Numerical Methods of Partial Differential Eqs，2000，16(3):379-394.

［9］Ge Yongbin，Cao Fujun，Zhang Jun.A transformation-free HOC scheme andmultigridmethod for solving the 3D poisson equation on nonuniform grids［J］.J Comput Phys，2013，234(1):199-216.

［10］ Ge Yongbin，Cao Fujun.Multigridmethod based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems ［J］.J Comput Phys，2011，230(10):4051-4070.

［11］Ge Yongbin.Multigridmethod and fourth-order compact difference discretization scheme with unequalmeshsizes for 3D poisson equation ［J］.J Comput Phys，2010，229(18):6381-6391.