蘭 斌, 王 濤

(1- 北方民族大學數(shù)學與信息科學學院，銀川 750021;2- 寧夏智能信息與大數(shù)據(jù)處理重點實驗室，銀川 750021)

1 引言

對流擴散方程是一類基本的運動方程，其廣泛應用于很多領(lǐng)域，如環(huán)境科學、能源開發(fā)、流體力學和電子科學等.在考慮輻射流體運動過程中，隨著溫度上升，光輻射逐漸在能量傳輸中占據(jù)主導地位，引起物質(zhì)狀態(tài)發(fā)生變化，變化的物態(tài)又反過來影響能量的傳輸.這是一類多介質(zhì)大變形非線性耦合特征鮮明的復雜問題，研究內(nèi)涵十分豐富，目前嚴格的理論研究難度大、進展較少，且實驗研究受實驗條件和成本的制約，因此，數(shù)值模擬是主要的研究手段.另外，在實際應用中，對具有極高密度比的物質(zhì)界面分辨精度的要求非?？量?，通常需要在拉氏(Lagrange)網(wǎng)格上進行計算；然而，在這個過程中，網(wǎng)格可能會隨著流體的運動而發(fā)生扭曲和變形，甚至有凹網(wǎng)格和翻轉(zhuǎn)網(wǎng)格的出現(xiàn)，這種大變形網(wǎng)格給具有間斷問題的非線性能量方程組的數(shù)值求解帶來極大困難.而現(xiàn)存的很多格式在求解該類問題時存在很多問題，諸如計算過程不收斂或收斂很慢，計算精度大幅度下降，出現(xiàn)非物理解或數(shù)值振蕩等等，最后導致數(shù)值結(jié)果失真而失去意義.另外，為適應一些實際問題的求解，諸如能量傳輸問題，要求溫度非負；質(zhì)量傳輸問題，要求濃度非負等等.因此，對于離散格式來說，保正性是一個關(guān)鍵點.在輸運與熱傳導中，一個不滿足保正性要求的格式很可能會違反熱力學第二定律的熵條件，引發(fā)熱流從低溫區(qū)流向高溫區(qū).Sharma 和Hammett[1]指出，在溫度變化較大的區(qū)域內(nèi)，容易導致溫度出負；但對于格式的健壯性來說，溫度的正性是必須的.

采用經(jīng)典的有限差分法或標準的有限元法往往得不到我們想要的結(jié)果.比如，傳統(tǒng)的迎風格式會因計算精度過低而易出現(xiàn)假擴散現(xiàn)象；而中心差分格式在一定條件下易引起非物理性的數(shù)值振蕩等等.為適應不規(guī)則區(qū)域問題的求解，采用有限體積法，即以方程的守恒形式出發(fā)，對計算域作一系列網(wǎng)格單元(控制體)剖分，然后在每個單元上進行積分，通過散度定理將其轉(zhuǎn)化為單元邊上的積分形式，再對邊上的通量作相應的離散.此方法能夠適應于各種守恒型問題的數(shù)值模擬，它能始終保持數(shù)值通量的局部守恒性，這一特性能夠清晰地顯示出它的優(yōu)越性.

近年來，科研工作者對扭曲網(wǎng)格上對流擴散方程滿足保正性要求的離散格式做了大量的工作.Le Potier[2]在非結(jié)構(gòu)三角形網(wǎng)格上提出了一個適用于強各向異性擴散系數(shù)的單元中心型非線性格式，但該格式只有當時間步長充分小的時候，系數(shù)矩陣才能在網(wǎng)格沒有限制性條件下滿足單調(diào)性.Lipnikov 等人[3]發(fā)展了Le Potier 的非線性格式，但卻依賴于擴散系數(shù)以及網(wǎng)格單元本身.后來，Yuan 和Sheng[4]提出任意星形網(wǎng)格上求解擴散方程的保正格式，且不需要對所選取的多邊形網(wǎng)格附加任何的正則性等限制性條件，同時，可很好地處理間斷、各向異性等問題，并顯式的給出了單元邊上法向通量的離散表達式，因此，從總體上大大地改進和推廣了上述的格式.Lipnikov 等人[5]基于文獻[4]中關(guān)于向量的分解辦法，給出一種不需要任何頂點輔助未知量的非線性兩點通量逼近方法.文獻[6]指出，已存在的部分單元中心型保正格式要求假設(shè)輔助未知量非負，但這并不一定總是成立，尤其是涉及到大變形網(wǎng)格情形.后來，文獻[7]構(gòu)造了一新的完全保正格式，即在無需假設(shè)輔助未知量非負的前提下可滿足保正性要求.關(guān)于對流項的離散，也出現(xiàn)了大量的工作[8-16]，他們在離散過程中，為避免數(shù)值振蕩，大都需要引入限制子技術(shù)，而部分問題的求解伴有降階情形.因此，文獻[17]提出了一種新的不含限制子技術(shù)的對流擴散方程的保正格式.

現(xiàn)考慮如下一維穩(wěn)態(tài)對流擴散方程

其中? 代表有界區(qū)域[l1,l2]，?? 是其邊界，κ(x), b(x)分別是擴散項和對流項的系數(shù)，f(x)為源項.

本文擬構(gòu)造任意非等距網(wǎng)格上的保正格式，旨在揭示文獻[15-17]中關(guān)于2D 對流擴散方程保正格式構(gòu)造的機理.其中，擴散通量的離散在非等距網(wǎng)格上是線性的，而在等距網(wǎng)格上，當擴散系數(shù)為標量時可退化為標準的二階中心差分格式；而對流通量的離散，為避免數(shù)值振蕩，對區(qū)間端點值的逼近，可在上游單元中心處作Taylor 級數(shù)展開，然后利用相關(guān)輔助未知量完成梯度值的重構(gòu)，將計算精度提高到二階，并對出負單元作正性校正，以滿足格式的保正性要求.

2 通量離散

將區(qū)間[l1,l2]任意剖分為互不重疊的子區(qū)間Ii:= [xi?1/2,xi+1/2], i = 1,2,··· ,N，即

Ii∩Ij當i ?=j 時最多只有一個交點，如圖1 所示，其中x1/2=l1, xN+1/2=l2.

圖1 網(wǎng)格模板圖

定義hi= xi+1/2?xi?1/2, i = 1,2,··· ,N，而h = max{hi, i = 1,2,··· ,N}，對方程(1)在區(qū)間Ii上進行積分，即

利用格林公式，得到

其中

2.1 擴散通量

考慮擴散通量的離散.可定義方程(4)中的前兩項為

略去上兩式的高階項后，擴散通量的離散形式可定義為如下具有二階精度的表達式，即

在兩端點處給出通量的守恒型條件，即

為保證區(qū)間端點處通量的守恒性，即有

于是，可分別解得區(qū)間端點值，即

然后，將式(7)代入式(5)中，經(jīng)整理，可得

最終，擴散通量的二階離散格式可整理為

其中

易證，線性方程(10)中的系數(shù)Ai,i?1, Ai,i, Ai,i+1均為正，且格式滿足離散極值原理.

接下來，考慮邊界單元的情形.當i=1 時，由上述推導，可得

其中

同理，當i=N 時，有

其中

特別指出，如果是等距網(wǎng)格剖分，且擴散系數(shù)是一個常數(shù)κ0，那么式(10)可寫成

其中即格式退化為標準的二階中心差分格式.

2.2 對流通量

由方程(4)出發(fā)，有

其中

同樣地

其中

于是，為匹配精度，對端點值的逼近可利用Taylor 級數(shù)展開法，在上游單元中心處進行展開至保留梯度項，即

類似地，可得

略去式(19)-(22)中的高階項后，連續(xù)對流通量(16),(17)求和后的離散格式可表示為

其中

并且

其中ui?1/2, ui+1/2等節(jié)點值的計算可由式(7)給出.

再考慮邊界單元.當i=1 時，由上述推導，可得

其中

同理，當i=N 時，有

其中

注1 為達到保正性要求，方程(23)中的相關(guān)系數(shù)必須滿足保正性要求，即

2.3 保正格式

2.3.1 有限體積格式

結(jié)合擴散通量的離散(10),(12)和(14)以及對流通量的離散(23),(25)和(26)，方程(1)和(2)的有限體積格式可以寫成如下形式

其中fi=f(xi), i=1,2,··· ,N.

2.3.2 非線性方程組

令U 為待求未知量所組成的向量，很顯然，三對角方程組(27)-(29)所形成的系數(shù)矩陣是非線性的，記為A(U)，那么，該方程組可表示為

而非對稱系數(shù)矩陣A(U)滿足如下性質(zhì)：

1) A(U)的所有主對角元素為正；

2) A(U)的所有非主對角元素為非正；

3) A(U)中所有列和均非負，并且至少存在一個列和的值為正.

因此，A(U)是按列弱對角占優(yōu)的.我們采用Picard 非線性迭代方法求解非線性系統(tǒng)(30)：任意選取較小的值Enon> 0 以及任意初始向量U0> 0，使用非線性迭代指標s=1,2,···，重復執(zhí)行以下過程：

1) 求解A(Us?1)Us=F；

2) 如果

則計算終止.

對于線性方程組A(Us?1)Us= F，可采用Gauss-Seidel 迭代法進行求解，直到殘差小于給定的值Elin，則計算終止.

2.3.3 算法

在這里，將給出具體的算法描述.

1) 任取初始U0>0, Enon和Elin.

2) 當s=0 時，

3) 當s=1,2,··· 時，

e) 計算殘差‖A(Us)Us?F‖2；

f) 當‖A(Us)Us?F‖2≤Enon‖A(U0)U0?F‖2時，則跳至4)，否則返回至3).

4) 終止.

2.3.4 保正性

在證明保正性前，首先引入如下引理[4,18].

引理1 對于滿足aii> 0(1 ≤i ≤n)和aij≤0(1 ≤i, j ≤n, i ?= j)的不可約矩陣A=(aij)n×n，如果A 是按列弱對角占優(yōu)的，即

且至少有一個不等式是嚴格大于0 的，則矩陣A 是一個M 矩陣[19,20].

于是，有如下定理成立.

定理1 令F ≥0, U0≥0，假設(shè)Picard 迭代中的線性方程組可以精確求解，則所有迭代向量Uk均為非負向量，即Uk≥0.

證明 在第2.3.2 小節(jié)中，我們給出了系數(shù)矩陣A(U)的一些性質(zhì).因此，它的轉(zhuǎn)置滿足引理1 中的條件，即AT(U)是一個M 矩陣，也即(AT(U))?1的所有元素均非負.因此，由(AT(U))?1= (A?1(U))T可知，A?1(U)也是非負的.如此，對任意的U ≥0, A(U)是一單調(diào)矩陣.

由于任意初始非負向量U0≥0，現(xiàn)假設(shè)正整數(shù)k0> 0，成立Uk0?1≥0.于是，系數(shù)矩陣A(Uk0?1)是單調(diào)矩陣，也即A?1(Uk0?1) ≥0 恒成立.當F ≥0 時，我們得到方程A(Uk0?1)Uk0= F 的解向量Uk0是完全非負的，即有Uk0≥0.因此，更具一般性地成立如下結(jié)論

Uk≥0, ?k ≥0.

3 數(shù)值實驗

本節(jié)將用兩組不同算例的數(shù)值結(jié)果來驗證本文格式的有效性.首先，引入L2范數(shù)

來定義計算解的截斷誤差，并取εnon= 1.0E ?6 及εlin= 1.0E ?10.而收斂階可由如下公式給出，即

其中N1和N2分別代表不同的網(wǎng)格單元數(shù)，而Error(N1)和Error(N2)分別代表其對應的L2誤差.對計算域作任意網(wǎng)格剖分，對任意內(nèi)節(jié)點坐標x，可取

x:=x+σξxh,

其中ξx為介于?0.5 到0.5 之間的隨機數(shù)，σ ∈[0,1]為擾動參數(shù)，h 為等距網(wǎng)格情形時的步長.

取?=[0,1]及b=1.0，問題的精確解以及擴散系數(shù)為

其中c=1,100; k0=1.0,10?6.

為了驗證格式的精確性，我們將分別在等距網(wǎng)格和非等距網(wǎng)格上進行問題的求解.表1 至表4 分別給出了當擴散系數(shù)連續(xù)和間斷時，在不同網(wǎng)格上數(shù)值解的L2誤差及其對應精度的比較.

從四個表中數(shù)據(jù)可以看出，無論是在等距網(wǎng)格上還是非等距網(wǎng)格上，當擴散系數(shù)連續(xù)(c = 1)和間斷(c = 100)時，均隨著擴散系數(shù)的降低，即問題的求解由擴散占優(yōu)逐漸向?qū)α髡純?yōu)情形轉(zhuǎn)變的過程中，數(shù)值解的精度均可達到二階.表明本文格式適用于任意網(wǎng)格上擴散占優(yōu)和對流占優(yōu)問題的求解.現(xiàn)取網(wǎng)格數(shù)為48 時，在四種情形下，數(shù)值解的最小值均為4.085E ?2，即格式滿足保正性要求.圖2 給出了非等距網(wǎng)格上，當網(wǎng)格點數(shù)取24，k0=10?6時，擴散系數(shù)連續(xù)(圖左)和擴散系數(shù)間斷(圖右)時的計算結(jié)果示意圖，可以看出，計算解與精確解吻合得很好.

表1 擴散系數(shù)連續(xù)時的計算結(jié)果(c=1, k0 =1.0)

表2 擴散系數(shù)連續(xù)時的計算結(jié)果(c=1, k0 =10?6)

表3 擴散系數(shù)間斷時的計算結(jié)果(c=100, k0 =1.0)

表4 擴散系數(shù)間斷時的計算結(jié)果(c=100, k0 =10?6)

圖2 非等距網(wǎng)格上計算結(jié)果示意圖

4 結(jié)論

本文構(gòu)造了任意非等距網(wǎng)格上一維穩(wěn)態(tài)對流擴散方程的保正格式. 其中，關(guān)于擴散項的離散，在等距網(wǎng)格上，當擴散系數(shù)為標量時，任意內(nèi)區(qū)間兩端點處的擴散通量可退化為標準的二階中心差分. 而對流項的離散，在處理區(qū)間端點值的逼近時，為避免數(shù)值振蕩而使其保持迎風特性，并提出一種新的方法使格式精度整體提高到二階. 該方法在上游單元中心處作泰勒級數(shù)展開，通過相關(guān)輔助未知量以完成梯度的重構(gòu)，并對出負情形作正性校正，使得格式滿足保正性要求.

該格式只含有區(qū)間單元中心未知量，并滿足區(qū)間端點處通量的局部守恒性. 數(shù)值實驗表明，本文格式對于處理擴散占優(yōu)、對流占優(yōu)問題，擴散系數(shù)連續(xù)和間斷情形均具有良好的適應性，并保持二階精度；并且，當網(wǎng)格退化為等距網(wǎng)格時，計算精度并未受到影響，這表明本文格式不受網(wǎng)格限制，可進行任意非等距網(wǎng)格剖分.