蘭 斌, 王 濤

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

1 引言

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

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

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

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

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

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

2 通量離散

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

Ii∩Ij當(dāng)i ?=j 時(shí)最多只有一個(gè)交點(diǎn)，如圖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}，對(duì)方程(1)在區(qū)間Ii上進(jìn)行積分，即

利用格林公式，得到

其中

2.1 擴(kuò)散通量

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

略去上兩式的高階項(xiàng)后，擴(kuò)散通量的離散形式可定義為如下具有二階精度的表達(dá)式，即

在兩端點(diǎn)處給出通量的守恒型條件，即

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

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

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

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

其中

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

接下來(lái)，考慮邊界單元的情形.當(dāng)i=1 時(shí)，由上述推導(dǎo)，可得

其中

同理，當(dāng)i=N 時(shí)，有

其中

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

其中即格式退化為標(biāo)準(zhǔn)的二階中心差分格式.

2.2 對(duì)流通量

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

其中

同樣地

其中

于是，為匹配精度，對(duì)端點(diǎn)值的逼近可利用Taylor 級(jí)數(shù)展開(kāi)法，在上游單元中心處進(jìn)行展開(kāi)至保留梯度項(xiàng)，即

類(lèi)似地，可得

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

其中

并且

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

再考慮邊界單元.當(dāng)i=1 時(shí)，由上述推導(dǎo)，可得

其中

同理，當(dāng)i=N 時(shí)，有

其中

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

2.3 保正格式

2.3.1 有限體積格式

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

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

2.3.2 非線性方程組

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

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

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

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

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

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

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

2) 如果

則計(jì)算終止.

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

2.3.3 算法

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

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

2) 當(dāng)s=0 時(shí)，

3) 當(dāng)s=1,2,··· 時(shí)，

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

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

4) 終止.

2.3.4 保正性

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

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

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

于是，有如下定理成立.

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

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

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

Uk≥0, ?k ≥0.

3 數(shù)值實(shí)驗(yàn)

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

來(lái)定義計(jì)算解的截?cái)嗾`差，并取εnon= 1.0E ?6 及εlin= 1.0E ?10.而收斂階可由如下公式給出，即

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

x:=x+σξxh,

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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