求解一類交界面問題的模態(tài)基函數(shù)譜元法數(shù)值實驗

邵文婷

(上海第二工業(yè)大學(xué)理學(xué)院,上海201209)

求解一類交界面問題的模態(tài)基函數(shù)譜元法數(shù)值實驗

邵文婷

(上海第二工業(yè)大學(xué)理學(xué)院,上海201209)

采用模態(tài)基函數(shù)形式下的譜元法,對一類一維和二維交界面問題進行了數(shù)值計算研究。不僅考慮了具有間斷系數(shù)的方程,還涉及了右端帶有奇性源項的方程。對于二維問題的數(shù)值計算,研究了四邊形單元譜元法和三角單元譜元法的離散格式。通過一些具有精確解的數(shù)值算例,驗證了基于模態(tài)基函數(shù)譜元法求解交界面問題的可行性和有效性。與其他數(shù)值方法相比,譜元法對于一維問題可以實現(xiàn)指數(shù)階收斂精度,對于二維問題也得到了較高的計算精度。

譜元法;四邊形單元;三角單元;交界面問題;間斷系數(shù);奇性源項

0 引言

具有間斷系數(shù)或者右端帶有奇性源項的方程稱為交界面問題。這類問題起源于許多應(yīng)用領(lǐng)域,例如2種不同材料或者相同材料在不同狀態(tài)下物理機制的研究[1]。由于方程所含的奇性,使得真解的光滑性較差,甚至在求解區(qū)域內(nèi)部會出現(xiàn)函數(shù)值間斷面。盡管在有限元、有限差分和有限體方法上,目前已提出了不少有效的數(shù)值離散格式,但對于高精度數(shù)值方法,如譜方法和徑向基函數(shù)方法在這類問題上的計算研究工作并不多見。Jung等[2]提出了基于最小二乘的Legendre譜配點法和徑向基函數(shù)方法來求解一維橢圓交界面問題。

對于譜方法或者擬譜方法,早期的工作[3]主要集中在單個區(qū)域上的Chebyshev譜方法和Fourier譜方法。由于這些高精度方法本身在建立離散格式時采用的基函數(shù)為定義在整個計算區(qū)域上的正交多項式,因此對真解在求解區(qū)域上局部呈現(xiàn)奇性的問題并不適用,對光滑性較差的函數(shù)的展開逼近具有一定的局限性。另一方面,由于本文所考慮的問題真解在物理區(qū)域交界面處具有奇性,從而破壞了光滑性,使得采用基于全局逼近的譜方法在進行數(shù)值計算時得不到指數(shù)階收斂精度。近年來,不少計算數(shù)學(xué)研究者致力于研究適用于多個區(qū)域上進行求解的高精度方法。譜元法(SEM)[4]正是在這種迫切需求下產(chǎn)生的,它將有限元法的靈活網(wǎng)格剖分技術(shù)融合于譜方法中,使得譜方法不再局限于整體逼近,對問題的真解在局部區(qū)域上呈現(xiàn)不同性態(tài)時也可以進行計算,同時保留指數(shù)階收斂精度。譜元法不僅保留了譜方法的高精度優(yōu)勢,同時具備了有限元法強健的區(qū)域適應(yīng)性。其基本想法是將求解區(qū)域剖分成若干個互不重疊的子單元,在每個單元上采用譜方法進行離散。同有限元法類似,將單元剛度和質(zhì)量矩陣按一定的順序集成總體剛度和質(zhì)量矩陣,最后求解離散得到的線性代數(shù)方程組。根據(jù)基函數(shù)的選取,SEM分為2種形式:節(jié)點基函數(shù)[5]和模態(tài)基函數(shù)[6]。有別于選取Gauss-Lobatto-Legendre(GLL)積分節(jié)點或者Chebyshev積分節(jié)點上建立的Lagrange多項式型節(jié)點基函數(shù),模態(tài)基函數(shù)是采用Legendre正交多項式和自然定義在三角區(qū)域上的Dubiner正交多項式[7],分別建立四邊形單元譜元法(QSEM)和三角單元譜元法(TSEM)。模態(tài)基函數(shù)相對于節(jié)點基函數(shù)的優(yōu)勢在于,充分利用了基函數(shù)的正交性,使得采用較少的計算量就可以有效地計算單元剛度和質(zhì)量矩陣,并且保證了高精度。

如果求解區(qū)域采用的網(wǎng)格剖分與問題給定的交界面相吻合,使得真解具有奇性的區(qū)域落在網(wǎng)格單元之間的邊界處,由于譜元法本身的逼近方式在單元邊界處的連續(xù)性僅為C0,當(dāng)真解在交界面處其導(dǎo)數(shù)出現(xiàn)跳躍時,而離散格式恰好允許這樣的非光滑性,因此采用譜元法求解這類交界面問題可以得到有效的數(shù)值結(jié)果?；谏鲜鱿敕?本文的主要工作是通過一系列具有精確解的數(shù)值算例,采用模態(tài)基函數(shù)形式下的譜元法進行求解交界面問題的數(shù)值實驗。對于二維問題,考慮了四邊形單元和三角單元2種譜元法。選取的算例中既涉及具有間斷系數(shù)的方程又涉及右端帶有奇性源項的方程。通過數(shù)值實驗發(fā)現(xiàn),對于一維問題,可以得到具有指數(shù)階收斂精度的數(shù)值解,對于二維問題也得到了較為滿意的數(shù)值結(jié)果。

1 譜元法離散格式

1.1 一維問題

首先,考慮如下一維交界面邊值問題

式中:r為常數(shù);β為分片常數(shù)。當(dāng)x∈?+時,β=β+;當(dāng)x∈?-時,β=β-。

對求解區(qū)域?進行網(wǎng)格剖分T,記?e為網(wǎng)格

?e上進行如下積分:

將單元?e變換至標(biāo)準(zhǔn)單元?st=[-1,1],那么可以將上式中每個單元?e上涉及的積分轉(zhuǎn)化為在?st上進行計算。假設(shè)在?st上,真解具有如下的N階多項式展開逼近:

其中:φk(ξ)=Lk(ξ)-Lk+2(ξ);Lk(ξ) 為 Legendre多項式。

采用高斯積分公式進行數(shù)值積分,最后得到展開系數(shù)滿足的如下線性代數(shù)方程組

其中:M為由單元剛度和質(zhì)量矩陣裝配后得到總體剛度和質(zhì)量矩陣構(gòu)成的系數(shù)矩陣,u為展開系數(shù)向量,u=,k=0,1,···,N,e=1,2,···,Nel;F為由右端項fe(x)確定的已知向量。

1.2 四邊形單元(QSEM)

對于二維問題,首先考慮在矩形求解區(qū)域上采用四邊形單元進行網(wǎng)格剖分,即:? =算,與一維問題相同,引入標(biāo)準(zhǔn)單元?st=[-1,1]2,從而使得在離散過程中每個單元?e上所涉及的微分和積分計算均可以在標(biāo)準(zhǔn)單元?st上一致進行。

假設(shè)在?st上,真解具有N階多項式展開:數(shù)。采用四邊形單元網(wǎng)格剖分的一大優(yōu)點在于,二維模態(tài)基函數(shù)可以由一維基函數(shù)進行張量積得到,ψi,j(ξ,η)= ψi(ξ)ψj(η),即[3]:

充分利用Legendre多項式的正交性,采用上述基函數(shù)建立的單元Laplacian算子矩陣和質(zhì)量矩陣具有如圖1所示的稀疏結(jié)構(gòu)。

圖1 QSEM離散得到的單元Laplacian算子矩陣和質(zhì)量矩陣的稀疏結(jié)構(gòu),多項式階數(shù)N=9Fig.1 Sparse structure of the elemental Laplacian and mass matrices of QSEM with the modal bases of order N=9

1.3 三角單元(TSEM)

假設(shè)任意三角單元?e的3個頂點坐標(biāo)為:{(xA,yA),(xB,yB),(xC,yC)},單元?e上的積分運算都可以通過如下坐標(biāo)變換轉(zhuǎn)化至參考三角單元τ={(ξ1,ξ2)|-1 ≤ξ1,ξ2;ξ1+ξ2≤ 0}上進行:

為了進一步計算方便,引入如下Duffy變換[8]及其逆變換

首先,通過式(8)、(9)將任意三角單元?e變換至參考三角單元τ,再借助Duffy變換,將參考三角單元τ變化至標(biāo)準(zhǔn)單元?st上,進行數(shù)值積分等相關(guān)計算(見圖2)。

圖2 任意三角單元(a)變換至參考三角單元τ(b),再變換至標(biāo)準(zhǔn)單元?st(c)進行數(shù)值積分Fig.2 Each triangular element(a)is mapped into the reference triangular τ;(b)then mapped into the standard element?st;(c)f i nally all local operations are evaluated in ?st

根據(jù)圖2(b)標(biāo)記的參考三角單元τ的頂點順序,給出如下三角單元的模態(tài)基函數(shù)[3]:

需要指出的是,Duffy變換是將標(biāo)準(zhǔn)單元?st的1條邊壓縮至三角形的1個頂點。如圖2所示,邊CC′通過變換后壓縮至頂點C,頂點C′退化,使得數(shù)值積分時在η2方向上采用的GLL積分節(jié)點在頂點C處過于密集。如果在η1方向上采用GLL積分節(jié)點,而在η2方向上采用Gauss-Radau積分節(jié)點,那么頂點退化的問題不會對數(shù)值積分引入奇性。

最后,圖3給出了單元Laplacian算子矩陣和質(zhì)量矩陣的稀疏結(jié)構(gòu)。

2 數(shù)值實驗

選取一系列帶有間斷系數(shù)或奇性源項的交界面問題,采用模態(tài)基形式下的譜元法來進行數(shù)值離散求解,通過計算結(jié)果來驗證譜元法求解這類問題的可行性和高效性。所有運算均借助Matlab 7數(shù)學(xué)軟件得以實現(xiàn),計算過程中所涉及的三角單元網(wǎng)格由Matlab自帶的PDE工具箱生成。

算例1 首先,考慮如下一維交界面問題,方程具有間斷系數(shù)且右端項含有奇性源Dirac函數(shù)[2]

精確解為

取μ=10,-10進行數(shù)值實驗。表1和表2列出了L2誤差隨著選取不同的多項式展開階數(shù)N和子區(qū)域個數(shù)Nel時的變化趨勢,并且與文獻[2]中采用的最小二乘Legendre譜配點法(LCM)得到的數(shù)值結(jié)果進行對比。從表中數(shù)據(jù)可以發(fā)現(xiàn),對該一維交界面問題,當(dāng)采用較少的未知量個數(shù)時,SEM得到了更高的精度,同時實現(xiàn)了指數(shù)階收斂速度。圖4和圖5分別呈現(xiàn)了μ=10,-10時,得到的數(shù)值解和節(jié)點上的誤差,這里取Nel=2,N=10。

圖3 TSEM離散得到的單元Laplacian算子矩陣和質(zhì)量矩陣的稀疏結(jié)構(gòu),多項式階數(shù)N=14Fig.3 Sparse structure of the elemental Laplacian and mass matrices of TSEM with the modal bases of order N=14

表1 LCM和SEM離散得到的L2誤差對比,μ=10,算例1Tab.1 Comparison of the L2norm errors for LCM and SEM withμ=10,Example 1

表2 LCM和SEM離散得到的L2誤差對比,μ=-10,算例1Tab.2 Comparison of the L2norm errors for LCM and SEM withμ=-10,Example 1

圖4 SEM計算得到的數(shù)值解和節(jié)點誤差,μ=10,Nel=2,N=10,算例1Fig.4 The numerical solution and pointwise errors of SEM withμ=10,Nel=2,N=10,Example 1

圖5 SEM計算得到的數(shù)值解和節(jié)點誤差,μ=-10,Nel=2,N=10算例1Fig.5 The numerical solution and pointwise errors of SEM withμ=-10,Nel=2,N=10,Example 1

算例2 考慮如下定義在?=(0,L)×(0,1)上的二維交界面問題[2]

精確解為:

其中:β-、β+、γ-、γ+、α、L、Ci的值與算例1中相同。

取μ=10進行數(shù)值測試,將求解區(qū)域?剖分為2 個子區(qū)域 [0,α]×[0,1],[α,L]×[0,1],采用 QSEM進行數(shù)值離散。在表3中,將其求解得到的L∞和L2誤差與LCM[2]進行了對比。由于LCM的微分矩陣條件數(shù)隨著配點個數(shù)的增加而增長得非常迅速,使得精度嚴重受到舍入誤差的影響,而本文采用的QSEM可以達到計算機的容忍精度。

表3 LCM[2]和SEM離散得到的L∞和L2誤差對比,μ=10,算例2Tab.3 Comparison of the L∞and L2norm errors for LCM and SEM withμ=10,Example 2

算例3 考慮如下二維交界面問題[9],采用TSEM進行數(shù)值求解

其中:交界面為曲線Γ=?-∩?+,?-=此外,[?u·n]=?u-·n-+?u+·n+,n±為?±的單位外法向量。

精確解為:

其中:α=3,r0=1/2,R=1/2。f,μ和g可以由精確解代入方程推導(dǎo)得到。

間斷系數(shù)β選取如下4組值進行數(shù)值測試:(1)β+=103,β-=1;(2)β+=106,β-=1;(3)β+=1,β-=103;(4)β+=1,β-=106。該算例的求解困難之處在于β+/β-的取值非常小或非常大。表4列出了當(dāng)多項式展開階數(shù)取為N=4,對應(yīng)不同大小的三角單元網(wǎng)格剖分,TSEM數(shù)值計算得到的L∞和L2誤差,其中h記為所有三角形單元的直徑最大值。從表中數(shù)據(jù)可以發(fā)現(xiàn),對于上述給出的4種間斷系數(shù)選取,無論β+/β-的值非常大還是非常小,TSEM均得到了具有相同數(shù)量級的精度,對參數(shù)的選取并不敏感。圖6展示了β+/β-=10-6,106時的數(shù)值解。

算例4 最后,考慮如下定義在?=[-1,1]2上,同時具有間斷變系數(shù)以及右端奇性源項的問題[10]

表4 TSEM離散得到的L∞和L2誤差,算例3Tab.4 The L∞and L2errors of TSEM,Example 3

圖6 β+/β-取值非常小和非常大時,TSEM計算得到的數(shù)值解,算例3Fig.6 The numerical solutions of TSEM with the very small and big ratio β+/β-,Example 3

式中,f(x,y)=8r2+4并且

精確解為:及其滿足的Dirichlet邊界條件。

表5列出了b=10,C=0.1和b=-3,C=0.1時,采用TSEM計算得到的L∞和L2誤差,圖7所示為對應(yīng)2組參數(shù)取值的數(shù)值解。

數(shù)值結(jié)果表明,采用TSEM求解同時具有間斷變系數(shù)和奇性源項的交界面問題是可行有效的,得到了不錯的數(shù)值結(jié)果。需要指出的是對于第二組參數(shù)值,由于b=-3＜0,問題本身并不是典型的橢圓方程,給求解帶來一定的難度,然而還是得到了有效的計算精度。

表5 TSEM離散得到的L∞和L2誤差,算例4Tab.5 The L∞and L2errors of TSEM,Example 4

圖7 TSEM計算得到的數(shù)值解,算例4Fig.7 The numerical solutions of TSEM,Example 4

3 結(jié) 語

本文對采用模態(tài)基函數(shù)形式下的譜元法求解交界面問題進行了數(shù)值實驗。對于二維問題的網(wǎng)格剖分,不僅考慮了四邊形單元還考慮了三角單元。選取一系列具有間斷系數(shù)和右端帶有奇性源項的交界面問題來驗證譜元法的計算效果。數(shù)值結(jié)果表明,譜元法對于交界面問題的求解是可行并且有效的,對于一維問題實現(xiàn)指數(shù)階收斂精度,對于二維問題也得到了較為理想的計算精度。

[1] WU Y,TRUSCOTT S,OKADA M.A f i nite difference scheme for an interface problem[J].Japan J Indust Appl Math,2010,27(2):239-262.

[2] SHIN B C,JUNG J H.Spectral collocation and radial basis function methods for one-dimensional interface problems[J].Appl Numer Math,2011,61(8):911-928.

[3] SHEN J,TANG T.Spectral and high-order methods with applications[M].2nd ed.Beijing:Science Press,2006.

[4] KARNIADAKIS G,SHERWIN S.Spectral/hp element methods for CFD[M].2nd ed.New York:Oxford University Press,2013.

[5] DEHGHAN M,SABOURI M.A spectral element method for solving the Pennesbioheat transfer equation by using triangular and quadrilateral elements[J].Appl Math Model,2012,36(12):6031-6049.

[6] FAKHAR–IZADI F,DEHGHAN M.A spectral element method using the modal basis and its application in solving second-order nonlinear partial differential equations[J].Math Method Appl Sci,2015,38(3):478-504.

[7] DUBINER M.Spectral methods on triangles and other domains[J].J Sci Comput,1991,6(4):345-390.

[8] DUFFY M G.Quadrature over a pyramid or cube of integrands with a singularity at a vertex[J].Siam Journal on Numerical Analysis,1982,19(6):1260-1262.

[9] HOU S,LIU X D.A numerical method for solving variable coeff i cient elliptic equation with interfaces[J].J Comput Phys,2005,202(2):411-445.

[10]LEVEQUE R J,LI Z.The immersed interface method for elliptic equations with discontinuous coeff i cients and singular sources[J].Siam J Numer Anal,1994,31(4),1019-1044.

A Numerical Study of The Modal Bases Spectral Element Method for Solving Interface Problems

SHAO Wenting(School of Science,Shanghai Polytechnic University,Shanghai 201209,China)

The numerical solution of one and two dimensional interface problems using the modal bases spectral element method(SEM)was studied.Equations with the discontinuous coeff i cient or singular source were concerned.For two dimensional cases,the implementations of the quadrilateral element(QSEM)and the triangular element(TSEM)were both investigated.Several numerical examples with exact solutions were provided to illustrate the feasibility and effectiveness of the modal bases SEM.Compared with other numerical results,exponential accuracy was regained for the one dimensional case,and better accuracy was also obtained for the two dimensional case.

spectral element method;quadrilateral elements;triangular elements;interface problems;discontinuous coeff i cient;singular source

O 241.82

A

1001-4543(2017)04-0283-08

10.19570/j.cnki.jsspu.2017.04.006

2017-04-27

邵文婷(1987–),女,上海人,講師,博士,主要研究方向為偏微分方程數(shù)值解。E-mail:wtshao@sspu.edu.cn。

上海市自然科學(xué)基金項目(16ZR1412700),上海第二工業(yè)大學(xué)校重點學(xué)科建設(shè)項目(XXKZD1304)資助