胡明皓 王莉華

(同濟(jì)大學(xué)航空航天與力學(xué)學(xué)院,上海 200092)

無(wú)網(wǎng)格法[1-6]由于不需要?jiǎng)澐志W(wǎng)格,不存在網(wǎng)格類方法在求解大變形問題時(shí)容易出現(xiàn)的網(wǎng)格畸變問題,而且具有精度高、收斂率高等優(yōu)點(diǎn),近年來(lái)受到越來(lái)越多的關(guān)注,廣泛應(yīng)用于高速?zèng)_擊、爆炸等復(fù)雜問題.常用的無(wú)網(wǎng)格法主要分為兩類: 伽遼金型和配點(diǎn)型.由于進(jìn)行區(qū)域積分,伽遼金型無(wú)網(wǎng)格法具有較好的精度和穩(wěn)定性,然而其缺點(diǎn)是計(jì)算效率比較低.配點(diǎn)型無(wú)網(wǎng)格法通常更加簡(jiǎn)單高效[7-11],而且在一些簡(jiǎn)單問題中可以獲得較好精度[12],但是對(duì)于一些復(fù)雜問題,其精度和穩(wěn)定性明顯降低.Zhang 等[13-14]提出最小二乘配點(diǎn)法(least squares collocation method,LSCM),通過(guò)采用比源點(diǎn)個(gè)數(shù)更多的配點(diǎn)來(lái)構(gòu)建超定離散方程進(jìn)行計(jì)算,顯著提升計(jì)算精度和穩(wěn)定性,然而超定方程求解也大幅降低計(jì)算效率.分區(qū)配點(diǎn)法[15-17]通過(guò)將區(qū)域分成若干個(gè)子域,提高了計(jì)算效率,而且降低了離散矩陣的條件數(shù),提高了結(jié)果的穩(wěn)定性,但是這種方法也是基于最小二乘求解,無(wú)法避免超定方程計(jì)算的缺點(diǎn).最近,Wang 等[18-20]提出一種基于重構(gòu)核近似(reproducing kernel,RK)的穩(wěn)定配點(diǎn)法(stabilized collocation method,SCM),這種方法將問題域劃分為若干個(gè)規(guī)則子域,對(duì)強(qiáng)形式方程在子域內(nèi)進(jìn)行積分,由于滿足積分約束,可以實(shí)現(xiàn)精確積分.該方法提高了無(wú)網(wǎng)格法計(jì)算結(jié)果的精度和穩(wěn)定性.由于積分效率高,該方法仍然保持配點(diǎn)型無(wú)網(wǎng)格法的高效特性.

常用的無(wú)網(wǎng)格法采用的形函數(shù)通常都是有理式,不具有插值特性,難以像有限元法(finite element method,FEM)一樣方便準(zhǔn)確地施加本質(zhì)邊界條件.這個(gè)問題也是目前無(wú)網(wǎng)格法的研究熱點(diǎn)之一.在無(wú)網(wǎng)格方法中常見的處理本質(zhì)邊界條件的方法有拉格朗日乘子法[21]、罰函數(shù)法[22]和一些修正的方法[23-24].然而拉格朗日乘子法會(huì)增加未知數(shù)的數(shù)量,導(dǎo)致剛度矩陣不對(duì)稱,從而增加了計(jì)算難度.罰函數(shù)法雖可得到對(duì)稱的剛度矩陣,但其計(jì)算精度往往取決于罰參數(shù)的選取.另一些學(xué)者努力嘗試將有限元法與無(wú)網(wǎng)格法相結(jié)合來(lái)施加邊界條件.基于強(qiáng)形式配點(diǎn)法和有限元的單元,Gao 等[25]提出單元微分法(element differential method,EDM),之后改進(jìn)等參單元的構(gòu)建形式,可以由配點(diǎn)和相鄰節(jié)點(diǎn)構(gòu)建,提出自由單元配點(diǎn)法(free element collocation method,FECM)[26-27].在此基礎(chǔ)上進(jìn)一步弱化單元,提出有限線法(finite line method,FLM)[28],這種方法在求解低維和高維問題時(shí)均只需要在一個(gè)方向上構(gòu)建單元,降低了單元畸變的可能性.這幾種方法都是基于直接配點(diǎn)法,求解比較簡(jiǎn)便,雖然可以利用有限元單元的優(yōu)勢(shì)方便地施加本質(zhì)邊界條件,但是沒能避免直接配點(diǎn)法的缺點(diǎn)和完全消除有限元單元畸變的可能性.

為了結(jié)合無(wú)網(wǎng)格法無(wú)單元畸變的優(yōu)勢(shì)和有限元法能方便施加本質(zhì)邊界條件的特點(diǎn),Wang 等[29]在穩(wěn)定配點(diǎn)法的基礎(chǔ)上提出一種基于拉格朗日插值的穩(wěn)定配點(diǎn)法,稱為拉格朗日插值穩(wěn)定配點(diǎn)法(stabilized Lagrange interpolation collocation method,SLICM).這種方法繼承了穩(wěn)定配點(diǎn)法精度高和穩(wěn)定性好的優(yōu)勢(shì),能夠?qū)崿F(xiàn)精確積分,同時(shí)由于拉格朗日插值形函數(shù)具有Kronecker delta 性質(zhì),可以像有限元法一樣簡(jiǎn)便準(zhǔn)確地直接施加本質(zhì)邊界條件.Wang等[29]考慮了均勻離散點(diǎn)和對(duì)應(yīng)結(jié)構(gòu)化網(wǎng)格的非均勻離散點(diǎn)的離散布置方案,并沒有詳細(xì)討論離散點(diǎn)任意布置的情形.本文通過(guò)引入曲線拉格朗日插值形函數(shù),實(shí)現(xiàn)了基于拉格朗日插值的直接配點(diǎn)法和穩(wěn)定配點(diǎn)法的任意離散,進(jìn)一步提升了這兩種方法的適用范圍.

1 拉格朗日插值近似

1.1 拉格朗日插值形函數(shù)

將一個(gè)一維區(qū)域離散為若干個(gè)離散點(diǎn),基于其中部分離散點(diǎn)x1,x2,···,xi,···,xm(這一組離散點(diǎn)也可稱為節(jié)點(diǎn)),一維拉格朗日插值多項(xiàng)式可表示為

其中m表示離散點(diǎn)個(gè)數(shù).對(duì)于所有的i≠I,NI(xI)在x=xi處結(jié)果為0,當(dāng)x=xI時(shí),NI(xI)=1.因此,拉格朗日插值具有Kronecker delta 性質(zhì),可以表示為以下形式

圖1 和圖2 展示了一維2 節(jié)點(diǎn)和3 節(jié)點(diǎn)拉格朗日插值形函數(shù)中節(jié)點(diǎn)分布情況.二維和三維形函數(shù)可以通過(guò)一維形函數(shù)在不同方向上的張量積表示為如下形式

圖1 一維2 節(jié)點(diǎn)和3 節(jié)點(diǎn)拉格朗日插值形函數(shù)的節(jié)點(diǎn)分布圖Fig.1 Node distributions of the 1D 2-node and 3-node Lagrange interpolation shape functions

圖2 二維4 節(jié)點(diǎn)、6 節(jié)點(diǎn)和9 節(jié)點(diǎn)拉格朗日插值形函數(shù)的節(jié)點(diǎn)分布圖Fig.2 Node distributions of the 2D 4-node,6-node and 9-node Lagrange interpolation shape functions

對(duì)于二維問題x=(x,y),三維問題x=(x,y,z),下標(biāo)L由下標(biāo)I ,J和K的順序排列確定.二維情況下4 節(jié)點(diǎn)、6 節(jié)點(diǎn)和9 節(jié)點(diǎn)拉格朗日形函數(shù)中節(jié)點(diǎn)分布如圖2 所示.

圖3～圖5 展示了不同階數(shù)下的拉格朗日插值形函數(shù)及其導(dǎo)數(shù),其中p表示拉格朗日插值形函數(shù)的階數(shù).對(duì)于一維問題p=m-1,二維問題p=min{(m1-1),(m2-1)},三維問題p=min{(m1-1),(m2-1),(m3-1)}.

圖3 p=2 時(shí)一維拉格朗日插值形函數(shù)Fig.3 1D Lagrange interpolation shape function when p=2

圖4 p=3 時(shí)一維拉格朗日插值形函數(shù)Fig.4 1D Lagrange interpolation shape function when p=3

1.2 拉格朗日插值形函數(shù)的導(dǎo)數(shù)

一維拉格朗日插值形函數(shù)的一階導(dǎo)數(shù)和二階導(dǎo)數(shù)可以表示為以下形式

二維形函數(shù)的一階和二階導(dǎo)數(shù)也可以通過(guò)一維形狀函數(shù)的張量積得到,可以表示為

同理,三維形函數(shù)的導(dǎo)數(shù)可表示為

1.3 拉格朗日插值近似

式(10)也可寫成如下形式

2 直接配點(diǎn)法(direct collocation method,DCM)

邊值問題可以表示為以下的一般形式

其中 Ω 表示待求問題的開區(qū)域,Γ 表示諾伊曼邊界,Π 表示狄利克雷邊界,整體求解區(qū)域=Ω∪Γ ∪Π.A,Bh和 Bg分別表示區(qū)域 Ω,Γ 和 Π 上的微分算子.u是未知量,f,h和g為對(duì)應(yīng)的源項(xiàng).

未知變量u可采用拉格朗日插值近似表示為

于是可以得到

將近似函數(shù)式(17)代入方程(14)～(16),并使其在域內(nèi)和邊界上滿足配點(diǎn)方程,可以得到如下離散形式

在直接配點(diǎn)法中,配點(diǎn)的位置與源點(diǎn)位置相同且數(shù)量一致.當(dāng)采用拉格朗日插值作為形函數(shù)時(shí),該方法被稱為拉格朗日插值配點(diǎn)法(Lagrange interpolation collocation method,LICM).系數(shù)矩陣a可以通過(guò)求解離散方程式(21)～式(23)來(lái)確定.離散方程式(21)～式(23)可改寫為以下的矩陣形式

式中子矩陣為

3 穩(wěn)定配點(diǎn)法(stabilized collocation method,SCM)

如圖6 所示,在穩(wěn)定配點(diǎn)法中,將域內(nèi)和邊界上布置若干個(gè)離散點(diǎn),以該離散點(diǎn)為中心構(gòu)建子域Ωl/Πl(fā)/Γl.子域的大小對(duì)于每個(gè)點(diǎn)可以是不同的,但為了簡(jiǎn)單起見一般使用相同大小的子域.子域的形狀對(duì)于二維問題一般選擇正方形,三維問題選擇立方體即可.

圖6 穩(wěn)定配點(diǎn)法中點(diǎn)的布置圖Fig.6 Points allocation in SCM

在穩(wěn)定配點(diǎn)法中,對(duì)方程式(14)～式(16)的強(qiáng)形式在相應(yīng)子域內(nèi)進(jìn)行積分,得到如下積分形式

將近似函數(shù)式(17)代入方程式(27)～式(29)可以得到

其中,p,q和r分別表示域內(nèi)、諾伊曼邊界和狄利克雷邊界上的離散點(diǎn).采用拉格朗日插值作為形函數(shù)的穩(wěn)定配點(diǎn)法稱為拉格朗日插值穩(wěn)定配點(diǎn)法(stabilized Lagrange interpolation collocation method,SLICM).同時(shí)對(duì)方程式(30)～式(32)等式兩邊做數(shù)值積分可得

其中,ω表示積分點(diǎn)個(gè)數(shù),wi表示積分權(quán)重,()l,()l和 ()l分別表示離散點(diǎn)pl ,ql和rl對(duì)應(yīng)子區(qū)域內(nèi)的積分點(diǎn).通?？梢允褂酶咚狗e分進(jìn)行數(shù)值積分運(yùn)算.式(33)～式(35)可以用矩陣形式表示如下

4 數(shù)值算例

本節(jié)采用一維、二維和三維3 個(gè)算例來(lái)評(píng)估基于拉格朗日插值的無(wú)網(wǎng)格直接配點(diǎn)法和穩(wěn)定配點(diǎn)法的精度和穩(wěn)定性.拉格朗日插值形函數(shù)的階數(shù)與影響域內(nèi)源點(diǎn)個(gè)數(shù)有對(duì)應(yīng)關(guān)系,3 個(gè)算例均采用二階形函數(shù),其影響域均為二階形函數(shù)所對(duì)應(yīng)的影響域范圍,即p階形函數(shù)在單個(gè)方向上對(duì)應(yīng)p+1個(gè)源點(diǎn).根據(jù)文獻(xiàn)[18]的研究結(jié)論,積分域大小和特征節(jié)點(diǎn)間距相距較近時(shí),可以得到比較高的精度.因此,本文算例中的積分域大小均參照特征節(jié)點(diǎn)間距.

4.1 一維直桿問題

考慮如圖7 所示左端固支的一維直桿的二階微分方程問題,其控制方程為

圖7 一維直桿示意圖Fig.7 Description of the 1D rod problem

其中EA=1,L=1,f(x)=3x2,解析解為.

圖8 為隨機(jī)布點(diǎn)的一維直桿離散圖,源點(diǎn)和配點(diǎn)布置在相同位置.圖9 展示了拉格朗日插值配點(diǎn)法、拉格朗日插值穩(wěn)定配點(diǎn)法和傳統(tǒng)重構(gòu)核近似配點(diǎn)法(reproducing kernel collocation method,RKCM)[30]在一維直桿問題中的數(shù)值計(jì)算結(jié)果和相應(yīng)的誤差,結(jié)果表明基于拉格朗日插值方法(LICM和SLICM)的求解精度明顯優(yōu)于基于重構(gòu)核近似方法(RKCM)的求解精度,拉格朗日插值穩(wěn)定配點(diǎn)法相比拉格朗日插值配點(diǎn)法有更好的計(jì)算精度,圖10的收斂性分析更清晰地展示了這一結(jié)論(圖標(biāo)中的數(shù)字為收斂率).由于源點(diǎn)和配點(diǎn)均隨機(jī)布置,收斂率和理論值略有差距,理論收斂率可參考文獻(xiàn)[29].圖11 比較了3 種方法的剛度矩陣條件數(shù),可以看出穩(wěn)定配點(diǎn)法的矩陣條件數(shù)更低,可以獲得更好的穩(wěn)定性.

圖8 一維直桿問題源點(diǎn)離散圖Fig.8 Discrete points of source points for the 1D rod problem

圖9 一維直桿問題計(jì)算結(jié)果及誤差Fig.9 Numerical solutions for the 1D rod problem

圖10 一維直桿問題收斂性分析Fig.10 Convergence comparisons for the 1D rod problem

圖11 一維直桿問題剛度矩陣條件數(shù)Fig.11 Condition number comparisons for the 1D rod problem

4.2 二維泊松問題

考慮二維空間中的泊松問題,控制方程和邊界條件如下

其中 Ω=(-1,1)×(-1,1),解析解為u(x,y)=exy.文獻(xiàn)[29]中研究了均勻離散和對(duì)應(yīng)結(jié)構(gòu)化網(wǎng)格布點(diǎn)的兩種離散方式,本文討論如圖12 所示任意非均勻離散模式.源點(diǎn)和配點(diǎn)布置在相同位置.

圖12 二維泊松問題源點(diǎn)離散圖Fig.12 Discrete points of source points for the 2D Poisson problem

圖13 和圖14 分別展示了采用拉格朗日插值配點(diǎn)法和拉格朗日插值穩(wěn)定配點(diǎn)法在域內(nèi)和邊界上的計(jì)算精度,邊界計(jì)算結(jié)果達(dá)到了計(jì)算機(jī)精度,表明這兩種方法都能夠準(zhǔn)確地施加本質(zhì)邊界條件.結(jié)合圖15的收斂性分析可以看出穩(wěn)定配點(diǎn)法對(duì)于離散點(diǎn)非均勻分布的二維泊松問題也能獲得更高的精度.圖16比較了兩種方法的剛度矩陣條件數(shù),穩(wěn)定配點(diǎn)法的條件數(shù)更小,表示該方法具有更好的穩(wěn)定性.

圖13 二維泊松問題域內(nèi)位移解Fig.13 Numerical solutions for the 2D Poisson problem

圖14 二維泊松問題邊界位移解Fig.14 Numerical solutions of boundary for the 2D Poisson problem

圖15 二維泊松問題收斂性分析Fig.15 Convergence comparisons for the 2D Poisson problem

圖16 二維泊松問題剛度矩陣條件數(shù)Fig.16 Condition number comparisons for the 2D Poisson problem

4.3 三維亥姆霍茲問題

考慮三維空間下亥姆霍茲方程,表示如下

其中 Ω=[-1,1]×[-1,1]×[-1,1],k=1,解析解為u(x,y,z)=sinx+siny+sinz.

將兩種方法進(jìn)一步擴(kuò)展至三維空間計(jì)算上述亥姆霍茲方程,可以得到類似的計(jì)算結(jié)果.圖17 展示了三維空間下的離散點(diǎn)布置圖,圖18 給出了域內(nèi)x=0上的計(jì)算結(jié)果.圖19 展示三維空間下的邊界計(jì)算結(jié)果,可以看出兩種方法對(duì)于邊界的計(jì)算精度都非常好.圖20 和圖21 分別展示直接配點(diǎn)法和穩(wěn)定配點(diǎn)法在三維亥姆霍茲問題上的收斂性分析和矩陣條件數(shù)分析.如前面兩個(gè)算例得到的結(jié)論相同,穩(wěn)定配點(diǎn)法有更好的精度和穩(wěn)定性.

圖17 三維亥姆霍茲問題源點(diǎn)離散圖Fig.17 Discrete points of source points for the 3D Helmholtz problem

圖18 三維亥姆霍茲問題域內(nèi)位移解Fig.18 Numerical solutions of the inner domain for the 3D Helmholtz problem

圖19 三維亥姆霍茲問題邊界位移解Fig.19 Numerical solutions of boundary for the 3D Helmholtz problem

圖19 三維亥姆霍茲問題邊界位移解 (續(xù))Fig.19 Numerical solutions of boundary for the 3D Helmholtz problem (continued)

圖20 三維亥姆霍茲問題收斂性分析Fig.20 Convergence comparisons for the 3D Helmholtz problem

圖21 三維亥姆霍茲問題剛度矩陣條件數(shù)Fig.21 Condition number comparisons for the 3D Helmholtz problem

5 總結(jié)

本文將曲線拉格朗日插值多項(xiàng)式引入無(wú)網(wǎng)格直接配點(diǎn)法和穩(wěn)定配點(diǎn)法中,構(gòu)建了離散點(diǎn)可任意布置的拉格朗日插值配點(diǎn)法和拉格朗日插值穩(wěn)定配點(diǎn)法.這兩種方法的形函數(shù)與有限元法一樣具有Kronecker delta 性質(zhì),能夠簡(jiǎn)單準(zhǔn)確地施加本質(zhì)邊界條件,相比較傳統(tǒng)無(wú)網(wǎng)格方法在邊界處理方面有較大優(yōu)勢(shì),提升了無(wú)網(wǎng)格法在邊界上的求解精度.拉格朗日插值穩(wěn)定配點(diǎn)法將問題域劃分為若干個(gè)子域,在子域內(nèi)能夠?qū)崿F(xiàn)精確積分,積分降低了剛度矩陣的條件數(shù),進(jìn)一步提高了算法的精度和穩(wěn)定性,解決了傳統(tǒng)無(wú)網(wǎng)格法難以實(shí)現(xiàn)精確積分的難題.該兩種方法由于采用曲線拉格朗日插值,離散點(diǎn)可任意布置,不再局限于結(jié)構(gòu)化網(wǎng)格布點(diǎn),使得這兩種方法可以適用于各類復(fù)雜區(qū)域問題的求解.基于拉格朗日插值的直接配點(diǎn)法和穩(wěn)定配點(diǎn)法的計(jì)算成本都比較低,邊界求解精度高,穩(wěn)定配點(diǎn)法可以通過(guò)精確積分進(jìn)一步提高精度和穩(wěn)定性,未來(lái)可以將這兩種方法應(yīng)用于更多力學(xué)問題和實(shí)際工程問題的數(shù)值分析.