邵文婷

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

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

邵文婷

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

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

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

0 引言

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

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

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

1 譜元法離散格式

1.1 一維問題

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

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

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

?e上進(jìn)行如下積分:

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

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

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

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

1.2 四邊形單元(QSEM)

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

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

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

圖1 QSEM離散得到的單元Laplacian算子矩陣和質(zhì)量矩陣的稀疏結(jié)構(gòu),多項(xiàng)式階數(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個(gè)頂點(diǎn)坐標(biāo)為:{(xA,yA),(xB,yB),(xC,yC)},單元?e上的積分運(yùn)算都可以通過如下坐標(biāo)變換轉(zhuǎn)化至參考三角單元τ={(ξ1,ξ2)|-1 ≤ξ1,ξ2;ξ1+ξ2≤ 0}上進(jìn)行:

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

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

圖2 任意三角單元(a)變換至參考三角單元τ(b),再變換至標(biāo)準(zhǔn)單元?st(c)進(jìn)行數(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)記的參考三角單元τ的頂點(diǎn)順序,給出如下三角單元的模態(tài)基函數(shù)[3]:

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

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

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

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

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

精確解為

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

圖3 TSEM離散得到的單元Laplacian算子矩陣和質(zhì)量矩陣的稀疏結(jié)構(gòu),多項(xiàng)式階數(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誤差對(duì)比,μ=10,算例1Tab.1 Comparison of the L2norm errors for LCM and SEM withμ=10,Example 1

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

圖4 SEM計(jì)算得到的數(shù)值解和節(jié)點(diǎn)誤差,μ=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計(jì)算得到的數(shù)值解和節(jié)點(diǎn)誤差,μ=-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進(jìn)行數(shù)值測(cè)試,將求解區(qū)域?剖分為2 個(gè)子區(qū)域 [0,α]×[0,1],[α,L]×[0,1],采用 QSEM進(jìn)行數(shù)值離散。在表3中,將其求解得到的L∞和L2誤差與LCM[2]進(jìn)行了對(duì)比。由于LCM的微分矩陣條件數(shù)隨著配點(diǎn)個(gè)數(shù)的增加而增長(zhǎng)得非常迅速,使得精度嚴(yán)重受到舍入誤差的影響,而本文采用的QSEM可以達(dá)到計(jì)算機(jī)的容忍精度。

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

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

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

精確解為:

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

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

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

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

圖6 β+/β-取值非常小和非常大時(shí),TSEM計(jì)算得到的數(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時(shí),采用TSEM計(jì)算得到的L∞和L2誤差,圖7所示為對(duì)應(yīng)2組參數(shù)取值的數(shù)值解。

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

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

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

3 結(jié) 語

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

[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–),女,上海人,講師,博士,主要研究方向?yàn)槠⒎址匠虜?shù)值解。E-mail:wtshao@sspu.edu.cn。

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