非等距網格上一維對流擴散方程的保正格式

蘭 斌, 王 濤

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

1 引言

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

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

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

現考慮如下一維穩態對流擴散方程

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

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

2 通量離散

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

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

圖1 網格模板圖

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

利用格林公式，得到

其中

2.1 擴散通量

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

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

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

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

于是，可分別解得區間端點值，即

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

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

其中

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

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

其中

同理，當i=N 時，有

其中

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

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

2.2 對流通量

由方程(4)出發，有

其中

同樣地

其中

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

類似地，可得

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

其中

并且

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

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

其中

同理，當i=N 時，有

其中

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

2.3 保正格式

2.3.1 有限體積格式

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

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

2.3.2 非線性方程組

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

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

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

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

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

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

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 是按列弱對角占優的，即

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

于是，有如下定理成立.

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

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

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

Uk≥0, ?k ≥0.

3 數值實驗

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

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

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

x:=x+σξxh,

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

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

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

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

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

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

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

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

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

圖2 非等距網格上計算結果示意圖

4 結論

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

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