不規(guī)則區(qū)域平面彈性問(wèn)題的正則區(qū)域重心插值配點(diǎn)法

， ， ，

(1.山東建筑大學(xué) 工程力學(xué)研究所，濟(jì)南 250101； 2.山東建筑大學(xué) 理學(xué)院，濟(jì)南 250101)

1 引 言

求解彈性力學(xué)平面問(wèn)題的主要方法有解析解法[1]和數(shù)值分析方法。對(duì)于不規(guī)則平面彈性問(wèn)題，很難得到其解析解，因此需要借助數(shù)值方法進(jìn)行求解。求解不規(guī)則區(qū)域平面彈性問(wèn)題的數(shù)值方法計(jì)算中，有限元法[2,3]、邊界元法[4]、有限差分法[5]和無(wú)網(wǎng)格法[6,7]是應(yīng)用最廣泛的方法。有限元法有其不足之處，需要進(jìn)行大量網(wǎng)格劃分，影響了有限單元法的計(jì)算精度和實(shí)際應(yīng)用范圍。邊界元法可只對(duì)不規(guī)則邊界進(jìn)行離散，離散化誤差僅僅來(lái)源于不規(guī)則邊界，但邊界元法依賴于所求問(wèn)題的基本解，需要進(jìn)行奇異積分計(jì)算。對(duì)于不規(guī)則問(wèn)題，采用有限差分法計(jì)算，需要構(gòu)造特殊的差分格式。近年來(lái)興起的無(wú)網(wǎng)格法有效彌補(bǔ)了有限單元法、邊界元法和有限差分法的不足之處?；趶较蚧瘮?shù)法的數(shù)值方法[7]將徑向基函數(shù)引入未知函數(shù)近似過(guò)程，該方法不需要?jiǎng)澐志W(wǎng)格,局限性在于計(jì)算精度非常依賴形狀參數(shù)的選取。

配點(diǎn)法是一種真正意義上的無(wú)網(wǎng)格方法，該方法是基于微分方程強(qiáng)形式的離散方法，程序?qū)嵤┓奖?，無(wú)需進(jìn)行網(wǎng)格劃分。常用的配點(diǎn)法有擬譜法[8]和微分求積法(DQM)[9]。擬譜法是基于譜函數(shù)近似展開(kāi)，數(shù)值求解微分方程的配點(diǎn)法。通過(guò)微分矩陣來(lái)表達(dá)未知函數(shù)的各階導(dǎo)數(shù)，微分矩陣通過(guò)對(duì)未知函數(shù)的近似函數(shù)求導(dǎo)得到。為了提高計(jì)算精度，在擬譜法中的節(jié)點(diǎn)一般采用譜多項(xiàng)式的零點(diǎn)來(lái)求解。微分求積法是由Bellman和Casti在20世紀(jì)70年代提出的求解偏微分方程的一種數(shù)值計(jì)算方法，其本質(zhì)是通過(guò)在全域節(jié)點(diǎn)的函數(shù)值加權(quán)近似得到未知函數(shù)的導(dǎo)數(shù)值。微分求積法的權(quán)一般通過(guò)Lagrange插值計(jì)算。擬譜方法和微分求積法不能直接應(yīng)用于不規(guī)則問(wèn)題的數(shù)值分析。

Lagrange插值是函數(shù)近似的重要方法，隨著插值節(jié)點(diǎn)數(shù)量的增加，Lagrange插值表現(xiàn)出極大的數(shù)值不穩(wěn)定性。將Lagrange插值公式改寫(xiě)成重心公式的形式，可以顯著降低插值的數(shù)值不穩(wěn)定性和改善計(jì)算精度。采用重心Lagrange插值近似未知函數(shù)的配點(diǎn)型方法即重心Lagrange插值配點(diǎn)法具有不需要?jiǎng)澐志W(wǎng)格、程序?qū)嵤┓奖愫陀?jì)算精度高的優(yōu)點(diǎn)，已廣泛應(yīng)用于求解線性和非線性微分方程的初值問(wèn)題、邊值問(wèn)題和初邊值問(wèn)題[10-16]。

本文采用重心插值配點(diǎn)法求解不規(guī)則區(qū)域彈性平面問(wèn)題。將不規(guī)則區(qū)域嵌入到規(guī)則的矩形區(qū)域，實(shí)現(xiàn)區(qū)域正則化。與有限單胞法相比[17]，本文方法為全局配點(diǎn)方法，不需要?jiǎng)澐钟?jì)算單元。在規(guī)則區(qū)域上，通過(guò)重心插值離散和施加邊界條件得到規(guī)則區(qū)域上的位移數(shù)值解，利用重心插值得到不規(guī)則區(qū)域內(nèi)任意節(jié)點(diǎn)的位移解。提出了一種研究不規(guī)則區(qū)域彈性平面問(wèn)題的正則區(qū)域法，為求解不規(guī)則彈性平面問(wèn)題提供了一種高精度的數(shù)值算法，通過(guò)數(shù)值算例驗(yàn)證了方法的有效性和計(jì)算精度。

2 重心Lagrange插值及其微分矩陣

2.1 一維重心Lagrange插值函數(shù)及其微分矩陣

對(duì)于區(qū)間[a,b]上的節(jié)點(diǎn)a=x1,x2,…,xn=b,函數(shù)p(x)在節(jié)點(diǎn)處的函數(shù)值為pj=p(xj),j= 1,2,…,n,則函數(shù)p(x)的重心Lagrange型插值函數(shù)為[16]

(j= 1,2,…,n) (1)

(m= 1,2,…)

(2)

寫(xiě)成矩陣的形式為

p(m)=D(m)p

(3)

2.2 二維重心Lagrange插值函數(shù)及其偏微分矩陣

假設(shè)函數(shù)p(x,y)在定義的矩形區(qū)域Ω0= [a,b]×[c,d]上，由節(jié)點(diǎn)a=x1,x2,…,xm=b,c=y1,y2,…,yn=d構(gòu)成Ω0上張量積型節(jié)點(diǎn)(xi,yj),對(duì)應(yīng)的函數(shù)值pi j=p(xi,yj)。函數(shù)p(x,y)的重心Lagrange插值為[16]

(4)

式中Li(x)和Mj(y)是一維重心Lagrange插值基函數(shù)。由式(4)，函數(shù)p(x,y)的(l+k)階偏微分可以寫(xiě)為

(5)

未知函數(shù)p(x,y)在節(jié)點(diǎn)(xi,yj),i= 1,2,…,m,j= 1,2,…,n處的(l+k)階偏微分值為

(p= 1,2,…,m,q= 1,2,…,n)

(6)

式(8)可以用矩陣形式表示為

P(l,k)= (C(l)?D(k))P

(7)

式中C(l)為節(jié)點(diǎn){xi}處l階重心Lagrange插值微分矩陣，?為矩陣的Kronecker積，D(k) 為節(jié)點(diǎn){yj}處k階重心Lagrange插值微分矩陣。定義C(0)=Im和D(0)=In,Im和In分別為m階和n階單位矩陣。向量P和P(l,k)分別為函數(shù)值及其(l+k)階偏微分值構(gòu)成的列向量。

3 不規(guī)則區(qū)域彈性平面問(wèn)題的正則區(qū)域重心Lagrange插值配點(diǎn)法

3.1 直角坐標(biāo)系下彈性力學(xué)的基本方程

在直角坐標(biāo)系下，平面問(wèn)題的位移表達(dá)控制方程為

(8)

位移邊界條件

(9)

位移形式表達(dá)的力邊界條件為

(10)

3.2 不規(guī)則區(qū)域彈性平面的正則區(qū)域法

正則區(qū)域法是將任意不規(guī)則物理區(qū)域Ω嵌入到一個(gè)規(guī)則的矩形計(jì)算區(qū)域Ω0= [a,b]×[c,d],將規(guī)則區(qū)域Ω0的幾何邊界條件記作Γ0,不規(guī)則區(qū)域Ω的幾何邊界條件記作Γ。在Ω0的x方向布置m個(gè)節(jié)點(diǎn)x1,x2,…,xm,y方向布置n個(gè)點(diǎn)y1,y2,…,yn,從而可以布置成m×n的張量型數(shù)值計(jì)算節(jié)點(diǎn)(xi,yj),i= 1,2,…,m,j= 1,2,…,n。如圖1所示，·為區(qū)域Ω內(nèi)的節(jié)點(diǎn)，°為區(qū)域Ω0外的節(jié)點(diǎn)，構(gòu)成計(jì)算節(jié)點(diǎn)。

利用插值式(4)近似位移函數(shù)u(x,y)和v(x,y),并將其代入位移控制方程(8),利用偏導(dǎo)數(shù)的微分矩陣表達(dá)式(7),位移控制方程(8)的重心插值配點(diǎn)法離散方程可寫(xiě)作如下矩陣形式。

(11)

引進(jìn)記號(hào)

方程組(11)可化簡(jiǎn)為

LU=F

(12)

式中u,v,fx和fy分別為位移和體力函數(shù)在計(jì)算節(jié)點(diǎn)處函數(shù)值構(gòu)成列向量。

3.3 邊界條件的離散及施加

(13)

式(13)的矩陣形式為

(14)

圖1 不規(guī)則區(qū)域Ω嵌入一個(gè)矩形區(qū)域Ω0

Fig.1 Irregular domain Ω embedded into a rectangular regionΩ0

力邊界條件(10)可采用類似的方式離散。將位移和力邊界條件的離散方程，簡(jiǎn)記為

B0U=g,B1U=h

(15)

采用附加法[15,16]施加邊界條件(15)，即將方程(15)附加到方程(12)，得到

(16)

方程組(16)為過(guò)約束線性方程組，采用最小二乘法求解。在數(shù)值算例實(shí)施過(guò)程中，采用MATLAB的反斜杠命令求解方程組(16)。為確保方程組(16)解的唯一性，選取的邊界節(jié)點(diǎn)數(shù)量應(yīng)大于正則區(qū)域邊界點(diǎn)的數(shù)量。

4 數(shù)值算例

本文數(shù)值算例，采用MATLAB進(jìn)行數(shù)值算例程序的編寫(xiě)，計(jì)算節(jié)點(diǎn)采用與Chebyshev節(jié)點(diǎn)同分布的點(diǎn)。不同數(shù)量節(jié)點(diǎn)計(jì)算的絕對(duì)誤差和相對(duì)誤差分別定義為Ea=‖Uc-Ue‖2,Er=‖Uc-Ue‖2/‖Ue‖2。Uc和Ue為計(jì)算值和解析解值列向量；‖?‖2表示向量的2范數(shù)。

算例1拱形區(qū)域Ω,三條直角邊固支，拱頂承受面力作用，各點(diǎn)坐標(biāo)分別為O(-1,0),B(1,0),C(1,1),D(-1,1),如圖2所示。拱頂曲線函數(shù)為拋物線y= 2-x2。拱頂上的面力和體力分量由位移解析解u(x,y) =y(x2-1)和v(x,y) =y(x2-1)確定。

圖2 拱形區(qū)域嵌入矩形區(qū)域

Fig.2 Arch domain embedded rectangular region

在本文的數(shù)值算例中，不規(guī)則區(qū)域內(nèi)的插值節(jié)點(diǎn)采用結(jié)構(gòu)化的布置方式，這有利于計(jì)算結(jié)果的可視化繪圖。不同數(shù)量計(jì)算節(jié)點(diǎn)的計(jì)算誤差列入 表1。由表1可知，采用5×5個(gè)計(jì)算節(jié)點(diǎn)即可達(dá)到極高的計(jì)算精度，隨著計(jì)算節(jié)點(diǎn)數(shù)量的增加，計(jì)算精度有所降低。這是因?yàn)楸舅憷慕鉃槎味囗?xiàng)式，重心Lagrange插值為多項(xiàng)式插值，采用高次多項(xiàng)式近似低階多項(xiàng)式的誤差，隨著近似多項(xiàng)式的次數(shù)增加而增大。

exp(x)(xcosy+sinx)

表1 拱形區(qū)域的位移絕對(duì)和相對(duì)誤差

Tab.1 Computational errors of displacement by barycentric interpolation collocation method under different numbers of nodes in example 1

計(jì)算節(jié)點(diǎn)數(shù)m=n絕對(duì)誤差Ea相對(duì)誤差Er52.3750×10-131.0065×10-14113.4157×10-121.4476×10-13151.5458×10-106.5509×10-12212.5137×10-91.0653×10-10

表2 花瓣型區(qū)域的位移絕對(duì)和相對(duì)誤差

Tab.2 Computational errors of displacement by barycentric interpolation collocation method under different numbers of nodes in example 2

計(jì)算節(jié)點(diǎn)數(shù)m=n絕對(duì)誤差Ea相對(duì)誤差Er91.3919×10-31.1896×10-5131.3971×10-71.1941×10-9171.1323×10-109.6777×10-13217.4544×10-106.3709×10-12

圖3 拱形區(qū)域插值節(jié)點(diǎn)分布圖

Fig.3 Interpolated nodes in arch domain

圖4 花瓣型區(qū)域嵌入矩形區(qū)域

Fig.4 Petals domain embedded rectangular region

圖5 花瓣型區(qū)域插值節(jié)點(diǎn)分布圖

Fig.5 Interpolated nodes in petals domain

圖6 含有花瓣形孔洞的圓形域嵌入矩形區(qū)域

Fig.6 Circular domain with a petal-shaped holes embedded rectangular region

算例3圓域含有花瓣型孔洞的多連通區(qū)域Ω=Ω2Ω1,如圖6所示，花瓣型孔洞區(qū)域Ω1的邊界曲線參數(shù)方程為

x= 1+(1+0.3sin(4θ))cosθ

y=-1+(1+0.3sin(4θ))cosθ(0≤θ≤2π)

圓形區(qū)域Ω2的邊界曲線方程為x2+y2= 3.52。位移邊界條件和體力分量由位移解析解確定[18]

u(x,y) = -y(2x+y)-exp(y)sin(x)+

v(x,y) = exp(y)cos(x)-xexp(x)cos(y)+

表3 圓域含有花瓣型孔洞的位移絕對(duì)和相對(duì)誤差

Tab.3 Computational error of displacement by barycentric interpolation collocation method under different types and numbers of nodes in example 3

計(jì)算節(jié)點(diǎn)數(shù)m=n絕對(duì)誤差Ea相對(duì)誤差Er174.6793×10-75.8865×10-10212.5344×10-103.1883×10-13274.9887×10-106.2758×10-13295.2532×10-106.6085×10-13

圖7 帶有花瓣形孔洞的圓形域計(jì)算節(jié)點(diǎn)分布圖

Fig.7 Interpolated nodes in circular domain with a petal-shaped hole

圖8 含有花瓣形空洞的圓形域x方向位移數(shù)值計(jì)算等值線圖

Fig.8 Contours of displacement inxdirection for circular domain with a petal-shaped hole

圖9 含有花瓣形空洞的圓形域y方向位移數(shù)值計(jì)算等值線圖

Fig.9 Contours of displacement in y direction for circular domain with a petal-shaped hole

5 結(jié) 論

重心插值配點(diǎn)法求解不規(guī)則區(qū)域上的彈性問(wèn)題，具有類似譜方法的高精度。采用較少數(shù)量的節(jié)點(diǎn)進(jìn)行計(jì)算，即可實(shí)現(xiàn)極高的計(jì)算精度。重心插值配點(diǎn)法是一種真正意義下的無(wú)網(wǎng)格方法，計(jì)算過(guò)程不需要背景網(wǎng)格的數(shù)值積分。重心插值配點(diǎn)法既可以求解簡(jiǎn)單的單連通區(qū)域問(wèn)題，也可以求解多連通區(qū)域問(wèn)題。

矩陣-向量形式表達(dá)的計(jì)算公式，極大地簡(jiǎn)化了編寫(xiě)計(jì)算程序的復(fù)雜性。附加法施加邊界條件，就是組合離散控制方程和邊界條件為過(guò)約束代數(shù)方程，該方法可以實(shí)現(xiàn)任意邊界條件的施加。

采用正則區(qū)域法計(jì)算時(shí)，應(yīng)選取包含物理區(qū)域的最小的正則區(qū)域作為計(jì)算區(qū)域。數(shù)值試驗(yàn)表明，采用最小的正則區(qū)域，可以提高計(jì)算精度。對(duì)于多連通域上的彈性問(wèn)題，也可以采用極坐標(biāo)下的正則區(qū)域法求解[13,14]。應(yīng)當(dāng)注意的是，在極坐標(biāo)系下的正則區(qū)域?yàn)閳A、扇形或環(huán)形區(qū)域。

參考文獻(xiàn)(References):

[1] 侯祥林,李 琦,鄭夕健.按位移求解彈性力學(xué)平面問(wèn)題的解析構(gòu)造解研究[J].計(jì)算力學(xué)學(xué)報(bào),2015,32(3):411-417.(HOU Xiang-lin,LI Qi,ZHENG Xi-jian.Study on analytic construction solutions about solving elasticity plane problems by displacement method[J].ChineseJournalofComputationalMechanics,2015,32(3):411-417.(in Chinese))

[2] 卞正寧,杜世森,羅建輝.彈性力學(xué)問(wèn)題的一階有限元解法[J].計(jì)算力學(xué)學(xué)報(bào),2015,32(4):544-547.(BIAN Zheng-ning,DU Shi-sen,LUO Jian-hui.First-order finite element method for elasticity[J].ChineseJournalofComputationalMechanics,2015,32(4):544-547.(in Chinese))

[3] 方錫武,賈丙輝.非匹配結(jié)點(diǎn)的有限單元等參插值方法及其應(yīng)用研究[J].計(jì)算力學(xué)學(xué)報(bào),2017,34(1):43-48.(FANG Xi-wu,JIA Bing-hui.Research on isoparametric interpolation for finite element with non-matching nodes and its application[J].ChineseJournalofComputationalMechanics,2017,34(1):43-48.(in Chinese))

[4] Zhang J,Qin X,Han X,et al.A boundary face method for potential problems in three dimensions[J].InternationalJournalforNumericalMethodsinEngineering,2009,80(3):320-337.

[5] Opr?al I,Zahradník J.Elastic finite -difference method for irregular grids[J].Geophysics,1999,64(1):240-250.

[6] 張 雄,劉 巖,馬 上.無(wú)網(wǎng)格法的理論及應(yīng)用[J].力學(xué)進(jìn)展,2009,39(1):1-36.(ZHANG Xiong,LIU Yan,MA Shang.Meshfree methods and their applications [J].AdvancesinMechanics,2009,39(1):1-36.(in Chinese))

[7] Kee B B T,Liu G R,Lu C.A least-square radial point collocation method for adaptive analysis in linear elasticity[J].EngineeringAnalysiswithBoundaryElements,2008,32(6):440-460.

[8] Weideman J A C,Reddy S C.A MATLAB differentiation matrix suite[J].ACMTransactionsonMathematicalSoftware,2000,26(4):465-519.

[9] Ma H,Qin Q H.An interpolation-based local differential quadrature method to solve partial differential equations using irregularly distributed nodes[J].CommunicationsinNumericalMethodsinEngineering,2008,24(7):573-584.

[10] Jiang J,Wang Z Q,Wang J H,et al.Barycentric rational interpolation iteration collocation method for solving nonlinear vibration problems[J].JournalofComputationalandNonlinearDynamics,2016,11(2),021001(1-13).

[11] 王兆清,莊美玲,姜 劍.非線性MEMS微梁的重心有理插值迭代配點(diǎn)法分析[J].固體力學(xué)學(xué)報(bào),2015,36(5):453-459.(WANG Zhao-qing,ZHUANG Mei-ling,JIANG Jian.Nonlinear MEMS microbeam analysis by barycentric rational interpolation iteration collocation method[J].ChineseJournalofSolidMechanics,2015,36(5):453-459.(in Chinese))

[12] 王兆清,綦甲帥,唐炳濤.奇異源項(xiàng)問(wèn)題的重心插值數(shù)值解[J].計(jì)算物理,2011,28(6):883-888.(WANG Zhao-qing,QI Jia-shuai,TANG Bing-tao.Numerical solution of singular source problems with barycentric interpolation[J].ChineseJournalofComputationalPhysics,2011,28(6):883-888.(in Chinese))

[13] 李樹(shù)忱,王兆清,袁 超.極坐標(biāo)系下彈性問(wèn)題的重心插值配點(diǎn)法[J].中南大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,44(5):2031-2040.(LI Shu-chen,WANG Zhao-qing，YUAN Chao.Barycentric interpolation collocation method for solving elastic problems[J].JournalofCentralSouthUniversity(ScienceandTechnology),2013,44(5):2031-2040.(in Chinese))

[14] Wang Z Q,Li S C,Ping Y,et al.A highly accurate regular domain collocation method for solving potential problems in the irregular doubly connected domains[J].MathematicalProblemsinEngineering,2014:397327.

[15] 王兆清，李淑萍.非線性問(wèn)題的重心插值配點(diǎn)法[M].北京：國(guó)防工業(yè)出版社，2015.(WANG Zhao -qing,LI Shu-ping.BarycentricInterpolationCollocationMethodforNonlinearProblems[M].Beijing:National Defense Industry Press,2015.(in Chinese))

[16] 李樹(shù)忱，王兆清.高精度無(wú)網(wǎng)格重心插值配點(diǎn)法[M].北京:科學(xué)出版社，2012.(LI Shu-chen,WANG Zhao -qing.HighAccuracyMeshlessCenterofAtta-chingMethod[M].Beijing:Science Press,2012.(in Chinese))

[17] Schillinger D,Ruess,M,Zander,N,et al.Small and large deformation analysis with the p-and B-spline versions of the finite cell method[J].ComputationalMechanics,2012,50(4),445-478.

[18] Martins N F,Rebelo M.A meshfree method for elasticity problems with interfaces[J].AppliedMathematicsandComputation,2013,219(22):10732-10745.