宋 洋,黃 偉

(上海大學(xué)理學(xué)院,上海200444)

在研究某些實(shí)際問題時(shí),往往需要考慮兩同心球所介區(qū)域內(nèi)偏微分方程模型的數(shù)值求解.例如氣象科學(xué)、海洋科學(xué)、地球物理和天體物理等領(lǐng)域,研究某些問題時(shí)常需要數(shù)值模擬流體在球形區(qū)域內(nèi)的運(yùn)動(dòng),而這類問題可轉(zhuǎn)化為球形區(qū)域內(nèi)偏微分方程的數(shù)值求解[1-2].隨著計(jì)算機(jī)技術(shù)的飛速發(fā)展,作為求解偏微分方程的數(shù)值方法之一的譜方法越來越受到人們的重視.譜方法的突出優(yōu)點(diǎn)是高精度.因此,本工作采用譜方法數(shù)值求解兩同心球所介區(qū)域內(nèi)的偏微分方程問題.

20世紀(jì)90年代,Cao等[3]和Guo等[4-5]發(fā)展了以球面調(diào)和函數(shù)為基函數(shù)的正交逼近方法,并應(yīng)用于求解球面上的渦度方程和低馬赫流動(dòng)方程等.這奠定了球面上譜方法的數(shù)學(xué)基礎(chǔ),也為進(jìn)一步研究?jī)蓚€(gè)同心球所介區(qū)域的問題提供了思路和手段.Guo等[6]和黃偉等[7]通過Jacobi譜方法和球面調(diào)和譜方法的結(jié)合解決了單位球內(nèi)Navier-Stokes方程的數(shù)值求解問題.夏文杰等[8]和鄧紅梅等[9]分別將混合Legendre-球面調(diào)和譜方法和混合Legendre-球面調(diào)和擬譜方法應(yīng)用于兩同心球所介區(qū)域上Fisher型方程的數(shù)值求解.

Navier-Stokes方程在研究不可壓縮流體運(yùn)動(dòng)中起著重要的作用.目前,已有許多算法用于Navier-Stokes方程的數(shù)值求解[6-7,10-18]本工作主要研究三維空間中兩同心球所介區(qū)域內(nèi)采用混合Legendre-球面調(diào)和譜方法數(shù)值求解Navier-Stokes方程的初邊值問題.

式中,U(r,λ,θ,t),P(r,λ,θ,t),f(r,λ,θ,t),U0(r,λ,θ)分別表示速度、壓力和密度之比、源項(xiàng)和初始速度,ν＞0為動(dòng)力黏性系數(shù),常數(shù)T＞0.

利用譜方法求解Navier-Stokes方程時(shí),找到滿足散度為0的正交基難度較大.為避免構(gòu)造散度自由的基函數(shù)的困難,文獻(xiàn)[6-7]在構(gòu)造數(shù)值求解單位球內(nèi)Navier-Stokes方程的混合譜格式時(shí),將不可壓縮條件?·U=0用人工壓縮方程

代替.盡管人工壓縮條件的引入會(huì)增加一個(gè)O(ε)階的逼近誤差,但能避免構(gòu)造散度自由的基函數(shù),也不需要通常在投影方法中附加的非物理邊界條件

與文獻(xiàn)[8]相比較,本工作所面對(duì)的問題更復(fù)雜,其中方程的形式從數(shù)量方程變?yōu)橄蛄糠匠?而且方程中主要的非線性項(xiàng)由U(U?1)變?yōu)?U·?)U.算法的實(shí)現(xiàn)中涉及非線性項(xiàng)的處理往往占用主要的計(jì)算資源,如果說前者的處理比較直接,那么后者就需要考慮在內(nèi)存和時(shí)間受到一定限制的條件下,如何平衡數(shù)組的大小與計(jì)算重復(fù)度的關(guān)系.本工作在具體的程序中較好地解決了這個(gè)問題.另外,由于空間自變量選擇了更適合區(qū)域形狀的球面坐標(biāo)形式,動(dòng)量方程中的速度應(yīng)按球面坐標(biāo)進(jìn)行分解,但是這種情況下對(duì)應(yīng)的3個(gè)分量方程形式上存在較大的差異,無論是數(shù)值分析還是實(shí)際計(jì)算都會(huì)面臨不小的挑戰(zhàn).為了解決上述問題,本工作采用速度的笛卡爾分解,使3個(gè)分量方程在形式上保持大部分的相似,而在相應(yīng)的程序設(shè)計(jì)上,則通過兩種坐標(biāo)系之間要素的轉(zhuǎn)換來匹配梯度和散度在運(yùn)算中出現(xiàn)的形式.

1 全離散混合譜格式

本工作以不可壓縮條件被人工壓縮條件(見式(2))替換后形成的方程形式為基礎(chǔ)給出相應(yīng)的全離散混合Legendre-球面調(diào)和譜格式.為了方便描述,首先引入一些必要的記號(hào).

令 χ(r)=r2,空間={v|v在 ? 上可測(cè),且＜ ∞}具有下列的內(nèi)積和范數(shù)：

考慮相關(guān)的正交投影.令M 和N為非負(fù)整數(shù),PM(I)表示所有次數(shù)不超過M 的代數(shù)多項(xiàng)式在I上的限制所組成的集合,特別地,WN(S)表示由中所有實(shí)值函數(shù)組成的子集.令

基于簡(jiǎn)化實(shí)際計(jì)算和數(shù)值分析的角度考慮,空間自變量用球面坐標(biāo)表示,速度向量按笛卡爾坐標(biāo)分解,即采用混合坐標(biāo)的處理方式.

若記ej為xj方向上的單位向量,v(j)表示向量函數(shù)在xj方向上的投影,那么相應(yīng)地,向量函數(shù)空間內(nèi)積正交投影

下面構(gòu)造全離散混合譜格式.為描述方便,令

令τ表示時(shí)間t方向上的步長(zhǎng),則在構(gòu)造格式時(shí),需要3個(gè)參數(shù)σ1,σ2,σ3,且滿足0≤σ1,σ2,σ3≤1.若記

則求解人工壓縮條件下Navier-Stokes方程的全離散混合Legendre-球面調(diào)和譜格式就轉(zhuǎn)換為求使對(duì)一切滿足如下等式

當(dāng)σ1=σ2=σ3=0時(shí),式(11)是顯格式;否則,在每個(gè)時(shí)間層上和p的值由線性方程組確定.特別地,取和w為所在子空間的基函數(shù),則該格式轉(zhuǎn)化為用基函數(shù)的組合形式表達(dá)和p時(shí)的組合系數(shù)所滿足的線性方程組.在實(shí)際計(jì)算中,選擇球面調(diào)和函數(shù)的實(shí)部和虛部作為WN(S)的基函數(shù),以此取代球面調(diào)和函數(shù).這樣不僅避免了復(fù)數(shù)運(yùn)算,節(jié)省了內(nèi)存和時(shí)間,而且在精度要求比較高的情況下為雙精度型變量的引入創(chuàng)造了條件.

2 數(shù)值結(jié)果

對(duì)參數(shù)σ1=0,σ2=σ3=1情況下的式(11)進(jìn)行數(shù)值實(shí)驗(yàn),并選取試驗(yàn)函數(shù)

(x1,x2,x3)T與(r,λ,θ)T滿足如下關(guān)系:

為了確保數(shù)值誤差、初始定解條件和右端項(xiàng)展開式系數(shù)的計(jì)算精確度,利用Gauss型求積公式,其中和并分別用和表示的零點(diǎn)和相應(yīng)的Christoffel數(shù).此外,令

圖1表示ν=10?4,ε=10?7和τ=τj時(shí)速度和壓力的相對(duì)誤差在t=1時(shí)隨M,N,τ的變化.由圖1可以看出:即使在M和N很小的情況下,數(shù)值結(jié)果也有較高的精度;當(dāng)τ減小或M=N增加時(shí),數(shù)值計(jì)算的誤差快速遞減.

通過變化小參數(shù)ε,比較t=1,M=N=7,ν=10?4及τ=2×10?4時(shí)的數(shù)值誤差EM,N,u(t)和EM,N,p(t),結(jié)果如圖2所示.由圖2可以發(fā)現(xiàn),對(duì)極大或極小的ε,數(shù)值精度可能會(huì)降低.這與文獻(xiàn)[7]中的結(jié)論一致,即當(dāng)ε與τ同階時(shí),可以得到較好的收斂速度.

圖1 速度和壓力的相對(duì)誤差隨M,N和τ的變化Fig.1 Relative errors of velocity and pressure varies with M,N and τ

圖2 速度和壓力的相對(duì)誤差隨ε的變化Fig.2 Relative errors of velocity and pressure varies with ε

3 結(jié)束語

本工作提出了一種求解兩同心球所介區(qū)域上Navier-Stokes方程的全離散Legendre-球面調(diào)和譜方法.數(shù)值實(shí)驗(yàn)的結(jié)果表明,本算法具有較高的精度.另外,本算法也可應(yīng)用于求解其他球形區(qū)域上的相關(guān)問題.