, , ,
(1.山東建筑大學(xué) 工程力學(xué)研究所,濟(jì)南 250101; 2.山東建筑大學(xué) 理學(xué)院,濟(jì)南 250101)
求解彈性力學(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ì)算精度。
對(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)
假設(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)成的列向量。
在直角坐標(biāo)系下,平面問(wèn)題的位移表達(dá)控制方程為
(8)
位移邊界條件
(9)
位移形式表達(dá)的力邊界條件為
(10)
正則區(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)成列向量。
(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ù)量。
本文數(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
重心插值配點(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.