胡明皓 王莉華

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

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

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

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

1 拉格朗日插值近似

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1.3 拉格朗日插值近似

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

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

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

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

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

于是可以得到

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

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

式中子矩陣為

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

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

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

在穩(wěn)定配點法中,對方程式(14)～式(16)的強形式在相應子域內進行積分,得到如下積分形式

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

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

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

4 數(shù)值算例

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

4.1 一維直桿問題

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

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

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

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

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

圖9 一維直桿問題計算結果及誤差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.文獻[29]中研究了均勻離散和對應結構化網格布點的兩種離散方式,本文討論如圖12 所示任意非均勻離散模式.源點和配點布置在相同位置.

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

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

圖13 二維泊松問題域內位移解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.

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

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

圖18 三維亥姆霍茲問題域內位移解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 總結

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