二維拋物型奇異攝動問題的移動網(wǎng)格方法

周 琴

(湖南涉外經(jīng)濟學(xué)院信息與機電工程學(xué)院，長沙 410205)

0 引言

二維拋物型奇異攝動問題在自然科學(xué)、工程等領(lǐng)域有著廣泛的應(yīng)用背景，它可以描述物質(zhì)擴散時的濃度變化規(guī)律，土壤力學(xué)中的滲流問題等等。該類問題的解通常在局部范圍內(nèi)會產(chǎn)生劇烈的變化，在均勻網(wǎng)格上獲取數(shù)值解不能很好地逼近實際問題。目前求解奇異攝動問題的非均勻網(wǎng)格方法主要有自適應(yīng)移動網(wǎng)格方法[?,?,?,?]和層適應(yīng)網(wǎng)格方法[5–6]等。其中自適應(yīng)移動網(wǎng)格方法通過網(wǎng)格的迭代變化使解變化較大的區(qū)域集中較多的網(wǎng)格點，能較好地反映解的性質(zhì)。

奇異攝動問題的自適應(yīng)移動網(wǎng)格方法研究中，關(guān)于二維問題的研究成果相對較少。文獻[?]針對一維奇異攝動對流擴散方程組問題，利用有限差分方法，提出了基于等弧長分布的自適應(yīng)移動網(wǎng)格方法，并給出了一階后驗誤差估計。文獻[2–3]分別對一維奇異攝動時變對流擴散問題和一階奇異攝動初值問題進行分析，建立了基于等分布弧長的移動網(wǎng)格方法。文獻[?]提出將適合問題的層適應(yīng)網(wǎng)格作為初始網(wǎng)格來進行網(wǎng)格移動，比較了兩種移動網(wǎng)格方法。文獻[7–8]針對一階的雙曲守恒律方程提出了基于等分布原理的移動網(wǎng)格方法。文獻[?,?,?,?,?]研究了自適應(yīng)移動網(wǎng)格方法在流體力學(xué)等問題中的應(yīng)用.

我們考慮一類含源項二維拋物型奇異攝動問題其中0＜ε ?1, Ω=[a1,b1]×[a2,b2]。函數(shù)f(x,y,t)為已知的源匯函數(shù)，設(shè)它滿足一定的正則性條件。本文將提出計算問題(??)的自適應(yīng)移動網(wǎng)格算法，并給出兩個算例的數(shù)值實驗結(jié)果。

1 數(shù)值格式

本節(jié)構(gòu)造移動網(wǎng)格上求解問題(??)的數(shù)值格式。設(shè)某時刻物理平面上網(wǎng)格分布為{(xj,k,yj,k)}, j= 0,1,··· ,N, k= 0,1,··· ,M，對應(yīng)數(shù)值解為{uj,k}。由于自適應(yīng)移動網(wǎng)格方法中的網(wǎng)格是迭代變化的，網(wǎng)格移動后產(chǎn)生的新網(wǎng)格為不規(guī)則的四邊形網(wǎng)格，不適用常規(guī)的數(shù)值方法作數(shù)值計算。因此，在某個時間層迭代的最終網(wǎng)格形成后，需要將不規(guī)則的四邊形網(wǎng)格轉(zhuǎn)化為計算平面上均勻分布的網(wǎng)格，以便于后續(xù)計算。記計算平面上均勻分布的網(wǎng)格為{(ξj,ηk)}, j= 0,1,··· ,N, k= 0,1,··· ,M。在計算過程中，計算平面上的網(wǎng)格總是保持不變,用來輔助求解問題。記

2 網(wǎng)格迭代方法和解的重新分布

3 移動網(wǎng)格算法和數(shù)值實驗

基于前文構(gòu)造的數(shù)值格式及網(wǎng)格迭代時解的重新分布方法，下面給出移動網(wǎng)格方法求解模型方程(??)的算法。

第4 步 如果tm+1≤T，令m=m+1，轉(zhuǎn)第2 步；否則，算法終止，得到T時刻的數(shù)值解。

下面根據(jù)上述算法對兩個算例進行數(shù)值實驗。

圖1 t=0.5, ε=0.001 時，例1 的自適應(yīng)移動網(wǎng)格分布

圖2 為例2 中t= 0.2, ε2= 0.0001 時，空間網(wǎng)格數(shù)分別為50×50 和100×100 節(jié)點的自適應(yīng)網(wǎng)格分布。該問題的解在邊界x=1, y=1 附近變化比較劇烈，由圖2 可見自適應(yīng)網(wǎng)格在該邊界層比較集中。這些結(jié)果充分說明了對于二維拋物型奇異攝動問題，通過移動網(wǎng)格方法求解能自適應(yīng)地進行局部加密，能較好地體現(xiàn)解在局部區(qū)域的特性。

圖2 t=0.2, ε2 =0.0001 時，例2 的自適應(yīng)移動網(wǎng)格分布

來度量。圖3(a)和圖3(b)分別為50×50 空間網(wǎng)格點的自適應(yīng)移動網(wǎng)格和均勻網(wǎng)格上求解例1的誤差圖像(取t=0.5, ε=0.001 時)。為了便于比較，圖3(a)選取的網(wǎng)格節(jié)點為計算平面上均勻分布的網(wǎng)格，該網(wǎng)格是由物理平面上的不規(guī)則網(wǎng)格經(jīng)過坐標變換映射產(chǎn)生的。比較兩個圖象可以看出，移動網(wǎng)格上求解的精度優(yōu)于均勻網(wǎng)格上求解的精度。

圖3 t=0.5, ε=0.001 時，例1 中自適應(yīng)移動網(wǎng)格與均勻網(wǎng)格上求解的逐點誤差比較

表1 至表4 分別列出了兩個算例中t= 0.1,0.2,0.3,0.4,0.5 時刻，網(wǎng)格剖分數(shù)分別取50×50、100×100、200×200 時，在自適應(yīng)移動網(wǎng)格上獲取的數(shù)值解與真解的L2誤差和最大模誤差。小參數(shù)ε的取值如表中所示。由表中數(shù)據(jù)可知該方法對于參數(shù)ε較小的情況是有效的，具有一定的收斂性。

表1 移動網(wǎng)格上求解例1 的L2 誤差(ε=0.001)

表2 移動網(wǎng)格上求解例1 的最大模誤差(ε=0.001)

表3 移動網(wǎng)格上求解例2 的L2 誤差(ε2 =0.0001)

表4 移動網(wǎng)格上求解例2 的最大模誤差(ε2 =0.0001)

綜上，對于含源項二維拋物型奇異攝動問題(??)，使用自適應(yīng)移動網(wǎng)格方法求解既能體現(xiàn)解在局部區(qū)域的特性，同時又能保證理想的求解精度。說明該方法對于二維拋物型奇異攝動問題的求解具有有效性和優(yōu)越性。