李樹忱，王兆清，袁超

(1. 山東大學(xué) 巖土與結(jié)構(gòu)工程研究中心，山東 濟(jì)南，250061；2. 山東建筑大學(xué) 工程力學(xué)研究所，山東 濟(jì)南，250101 )

彈性力學(xué)是研究巖石變形、斷裂破壞問題的重要基礎(chǔ)，對于復(fù)雜的巖石彈性力學(xué)問題一般很難得到問題的解析解，數(shù)值計算方法成為一種重要的研究手段。研究彈性力學(xué)問題的基本數(shù)值方法有有限元方法和有限差分法。有限元方法是一種基于低階多項式插值的Galerkin方法，其計算精度依賴于單元大小和積分點數(shù)量的選?。挥邢薏罘址ú捎貌罘炙阕咏莆⒎炙阕?，其計算精度依賴于差分步長的大小。近年來興起的各種類型無網(wǎng)格方法，在巖石力學(xué)領(lǐng)域得到廣泛關(guān)注[1-3]。能夠采用較少的計算代價，得到問題的高精度數(shù)值解，是數(shù)值分析領(lǐng)域研究者努力尋求的目標(biāo)?；谝苿幼钚《朔ǖ臒o單元Galerkin法盡管有不需要劃分計算網(wǎng)格的優(yōu)點，但在形成剛度矩陣過程中，需要進(jìn)行大量的計算，并且不能直接施加邊界條件[3-4]。基于自然鄰點插值的自然單元法，邊界條件施加方便，但在形成剛度矩陣過程中，需要尋找自然鄰點，計算工作量巨大[5]。Lagrange插值是函數(shù)近似的重要手段，但當(dāng)插值節(jié)點數(shù)量加大時，Lagrange插值表現(xiàn)為數(shù)值不穩(wěn)定，Runge現(xiàn)象說明了這一問題，因此，工程計算中較少采用 Lagrange插值作為近似未知函數(shù)的手段。將Lagrange插值公式改寫為重心公式形式，得到重心Lagrange插值[6]。重心Lagrange插值具有很好的數(shù)值穩(wěn)定性[7]和極高的插值近似精度[8-9]。采用重心Lagrange插值作為未知函數(shù)的近似函數(shù)的配點型方法 —— 重心插值配點法，不需要劃分計算網(wǎng)格，也不需要數(shù)值積分，已應(yīng)用于求解一維工程邊值問題[10]和動力學(xué)問題[11-12]中?，F(xiàn)有的大部分?jǐn)?shù)值方法是在直角坐標(biāo)系下求解彈性力學(xué)問題。對于一些圓形區(qū)域上的工程問題，例如巖石力學(xué)的巴西圓盤實驗，采用極坐標(biāo)系分析非常方便[13-14]。在此，本文采用重心Lagrange插值的張量積形式作為二維問題未知函數(shù)的近似函數(shù)，基于極坐標(biāo)系下彈性力學(xué)位移平衡方程的強(qiáng)形式，提出求解極坐標(biāo)系下彈性力學(xué)問題的重心插值配點法。

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

考慮定義在區(qū)間[a,b]上的函數(shù)u(x)，函數(shù)u(x)在節(jié)點a=x1＜x2＜…＜xn=b上的取值為ui=u(xi),i=1,2, …,n。函數(shù)u(x)的重心Lagrange插值為[8]

函數(shù)u(x)的重心Lagrange插值可表示為

則函數(shù)u(x)的m階導(dǎo)數(shù)可表示為

函數(shù)u(x)在節(jié)點x1,x2, …,xn處的m階導(dǎo)數(shù)可表示為

2 極坐標(biāo)系下彈性問題的重心插值配點法公式

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

取極坐標(biāo)系(r,θ)，某點P(r,θ)處的位移分量為ur和uθ，應(yīng)力分量為σr和σθ和τrθ=τθr，體力分量為fr和fθ，則極坐標(biāo)系下的平衡方程可表示為

位移-應(yīng)變的幾何關(guān)系為：

應(yīng)力-應(yīng)變本構(gòu)關(guān)系為：

將式(9)和(11)代入式(8)，得到位移表達(dá)的平面應(yīng)力問題平衡方程：

位移表達(dá)的應(yīng)力分量為

對于平面應(yīng)變問題，需將式(10)～(12)中的彈性常數(shù)進(jìn)行下列替換：

2.2 重心插值配點法計算公式

極坐標(biāo)系下平面問題的位移解法需要求解偏微分方程(式(12))的邊值問題。對于幾何形狀和外載荷不隨θ改變的軸對稱結(jié)構(gòu)，其應(yīng)力分布與θ無關(guān)，只是坐標(biāo)r的函數(shù)，方程(12)將變?yōu)槌Ｎ⒎址匠?。軸對稱結(jié)構(gòu)常微分方程邊值問題的重心插值配點法求解方法，參見文獻(xiàn)[10]。

在求解區(qū)域?的r方向布置m個計算節(jié)點r1,r2, …,rm，在θ方向布置n個計算節(jié)點θ1,θ2, …,θn，構(gòu)成求解區(qū)域上的m×n個張量積計算節(jié)點(ri,θj)(其中，i=1, 2, …,m;j=1, 2, …,n)。采用節(jié)點r1,r2, …,rm和θ1,θ2, …,θn上的重心 Lagrange 插值(式(1))的張量積形式插值近似徑向位移ur和環(huán)向位移uθ：

其中：urij和uθij分別為徑向位移ur和環(huán)向位移uθ在計算節(jié)點(ri,θk)處的值；uθij=uθ(ri,θj) ；i=1, 2, …,m;j=1, 2, …,n；Li(r),Mj(θ)分別為在節(jié)點r1,r2, …,rm和θ1,θ2, …,θn上的重心 Lagrange 插值基函數(shù)。

將式(15)代入位移表達(dá)的平衡方程(式(12))，得：

令式(16)在所有的計算節(jié)點(rk,θl)(k=1, 2, …,m;l=1, 2, …,n)成立，可得：

將徑向位移ur、環(huán)向位移uθ和體力fr，fθ在計算節(jié)點(rk,θl),k=1, 2, …,m;l=1, 2, …,n處的值，排列成1個向量，表示為

其中：frij=fr(ri,θj) ；fθij=fθ(ri,θj) ；i=1, 2, …,m;j=1, 2, …,n。利用微分矩陣記號和矩陣張量積(Kronecker積)運算符號“?”[15]，微分方程(式(17))的矩陣表達(dá)式為：

其中：C(m)和D(m)分別為在節(jié)點r1,r2, …,rm和θ1,θ2, …,θn上的m階重心Lagrange插值微分矩陣；Im和In分別為m階和n階單位矩陣；diag(1/r)，diag(1/r2)分別為函數(shù)1/r和1/r2在節(jié)點r1,r2, …,rm處函數(shù)值構(gòu)成的對角矩陣。引進(jìn)記號：

將式(19)簡記為：

則式(21)的解耦形式為

式(22)可簡記為

按照平衡方程離散的同樣方式，位移表達(dá)的應(yīng)力分量(式(13))的重心插值離散形式表示為

由式(24)抽取邊界節(jié)點的應(yīng)力分量表達(dá)式，得到應(yīng)力邊界條件的重心插值離散計算公式。邊界條件的施加可以采用置換法或附加法[10-12]。求解代數(shù)方程方程組得到計算節(jié)點的位移值后，采用式(24)可以方便地計算出所有節(jié)點的應(yīng)力。

3 數(shù)值算例

數(shù)值算例采用 MATLAB編寫，以充分發(fā)揮MATLAB強(qiáng)大的矩陣運算能力。數(shù)值計算的絕對誤差和相對誤差分別定義為：

其中：Uc和Ue分別為數(shù)值計算解和解析解列向量；表示向量的 2范數(shù)。計算節(jié)點采用特殊分布的Chebyshev 節(jié)點[15]。

算例1曲梁的一般彎曲問題。考慮圖1所示矩形等截面曲梁受端部載荷作用的平面應(yīng)力問題。

曲梁的位移分量解析解為[16]

圖1 端部受集中力作用的曲梁Fig.1 Curved beam subjected to concentrated force on end

曲梁的應(yīng)力分量解析解為

曲梁的邊界條件的Saint-Venant提法為：

為得到問題的高精度解答，邊界條件(式(29)和(30))按照應(yīng)力解析解(式(27))和位移解析解(式(26))施加。

數(shù)值計算時，曲梁的幾何參數(shù)和力學(xué)參數(shù)分別取a=13，b=17，E=107，μ=0.25，P=104，單位一律采用國際單位制。不同數(shù)量Chebyshev節(jié)點重心Lagrange插值配點法計算的位移、應(yīng)力誤差見表1。由表1可知：采用較少數(shù)量的計算節(jié)點重心插值配點法計算的位移和應(yīng)力具有很高的計算精度，并且隨著計算節(jié)點數(shù)量的增加，位移和應(yīng)力的計算精度隨著提高；重心插值配點法計算的應(yīng)力誤差大于位移誤差，這是由于計算應(yīng)力需要數(shù)值求導(dǎo)造成計算精度損失。重心插值配點法計算的曲梁位移分量和應(yīng)力分量分布分別見圖 2～6。

算例 2含有孔洞無限大板的拉伸，應(yīng)力集中因子問題?？紤]如圖7所示含有孔洞受單向拉伸的無限大板問題，孔洞的半徑為a，板邊均布拉力為P。設(shè)板的厚度為1，按照平面應(yīng)力計算。

建立圖7所示坐標(biāo)系，極坐標(biāo)系下板內(nèi)應(yīng)力分量的解析解為[17]：

表1 不同數(shù)量計算節(jié)點重心插值配點法計算的位移和應(yīng)力誤差Table 1 Computational error of displacement and stress by barycentric interpolation collocation method under different types and numbers of nodes in Example 1

圖2 曲梁的徑向位移分布Fig.2 Distribution of radial displacement in curved beam

圖3 曲梁的環(huán)向位移分布Fig.3 Distribution of hoop displacement in curved beam

圖4 曲梁的徑向正應(yīng)力分布Fig.4 Distribution of radial normal stress in curved beam

圖5 曲梁的環(huán)向正應(yīng)力分布Fig.5 Distribution of hoop normal stress in curved beam

圖6 曲梁的剪應(yīng)力分布Fig.6 Distribution of shear stress in curved beam

圖7 含孔洞無限大板Fig.7 A infinite plate with hole under uniaxial tension

位移解析解為：

其中：G=E/ [2(1 +μ)]為剪切模量；κ= ( 3 -μ)/(1 +μ)；μ為泊松比。

圖8 重心插值配點法的計算區(qū)域Fig.8 Computational domain in barycentric interpolation collocation method

r=b邊界上應(yīng)力邊界條件按照式(33)確定：

最大環(huán)向應(yīng)力 (σθ)max在孔洞邊界，應(yīng)力集中因子Scf=(σθ)max/P的解析值等于3，重心插值配點法計算的應(yīng)力集中因子等于3.000 2。

在r方向和θ方向取不同數(shù)量計算節(jié)點計算位移和應(yīng)力的相對誤差，448種不同節(jié)點數(shù)量組合和相對誤差常用對數(shù)的關(guān)系分別見圖9～13。由圖9～13可以看出：隨著計算節(jié)點數(shù)量的增加，位移和應(yīng)力的計算誤差顯著降低；在相同的徑向節(jié)點數(shù)量條件下，增加環(huán)向節(jié)點數(shù)量不能明顯提高計算精度；在相同的環(huán)向節(jié)點數(shù)量下，增加徑向節(jié)點數(shù)量可以顯著提高計算精度；在環(huán)向采用10個計算節(jié)點、徑向采用33個計算節(jié)點的計算精度可達(dá)10-6量級。

上述2個算例說明：對于幾何區(qū)域為圓形、圓環(huán)或扇形的彈性力學(xué)問題，采用極坐標(biāo)解答可以有效地簡化問題的復(fù)雜性。對于含有圓形孔洞的無限大板拉伸問題，在板中截取1個與圓孔同心的大圓，將求解區(qū)域轉(zhuǎn)化為圓環(huán)區(qū)域，利用直角坐標(biāo)系和極坐標(biāo)系下應(yīng)力的轉(zhuǎn)換關(guān)系，可以方便地得到圓環(huán)求解區(qū)域的應(yīng)力邊界條件。

對于一般形狀或不規(guī)則形狀的物理區(qū)域，不能直接采用重心插值配點法計算，可以將物理區(qū)域嵌入 1個圓形或扇形的正則區(qū)域，在正則區(qū)域上采用極坐標(biāo)系的重心插值配點法計算。計算得到正則區(qū)域上點的位移和應(yīng)力值后，通過重心插值式(15)可以得到物理區(qū)域上任意一點的位移和應(yīng)力。

表2 重心插值配點法計算的位移和應(yīng)力誤差Table 2 Computational error of displacement and stress in barycentric interpolation collocation method

圖9 徑向位移相對誤差與計算節(jié)點數(shù)量的關(guān)系Fig.9 Relationship between relative error of radial displacement and number of computed nodes

圖10 環(huán)向位移相對誤差與計算節(jié)點數(shù)量的關(guān)系Fig.10 Relationship between relative error of hoop displacement and number of computed nodes

圖11 徑向應(yīng)力相對誤差與計算節(jié)點數(shù)量的關(guān)系Fig.11 Relationship between relative error of radial stress and number of computed nodes

圖13 剪應(yīng)力相對誤差與計算節(jié)點數(shù)量的關(guān)系Fig.13 Relationship between relative error of shear stress and number of computed nodes

4 結(jié)論

(1) 采用重心Lagrange插值作為近似函數(shù)的重心插值配點法，在極坐標(biāo)系下求解彈性力學(xué)問題，利用重心插值微分矩陣和矩陣張量積運算符號，可以將位移表達(dá)的極坐標(biāo)系平衡方程的離散代數(shù)方程表示為簡潔的矩陣形式。

(2) 將原來徑向位移和環(huán)向位移耦合的偏微分方程組，通過矩陣操作實現(xiàn)解耦，這為計算程序的實施提供了極大地便利。微分矩陣的引進(jìn)，提供了一種直接計算節(jié)點處應(yīng)力值的簡便快速方法。

(3) 重心插值配點法不需要劃分計算網(wǎng)格，不需要數(shù)值積分，是一種真正意義上的無網(wǎng)格方法。

(4) 重心插值配點法具有計算格式簡單、程序?qū)嵤┓奖愫陀嬎憔雀叩奶攸c。

[1] 盧波, 丁秀麗, 鄔愛清. 無網(wǎng)格法對巖體不連續(xù)面的模擬[J].巖石力學(xué)與工程學(xué)報, 2008, 27(10): 2108-2117.

LU Bo, DING Xiuli, WU Aiqing. Modeling of rock discontinuity with meshless method[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(10): 2108-2117.

[2] 樊成, 欒茂田, 黎勇, 等. Kriging插值無網(wǎng)格方法及其在力學(xué)邊值問題中的應(yīng)用[J]. 巖石力學(xué)與工程學(xué)報, 2007, 26(1):195-200.

FAN Cheng, LUAN Maotian, LI Yong, et al. A new type meshless method based on Kriging interpolation scheme and its application to solving boundary-value problem of mechanics[J].Chinese Journal of Rock Mechanics and Engineering, 2007,26(1): 195-200.

[3] 李樹忱, 李術(shù)才, 隋斌, 等. 節(jié)理巖體滲流的無網(wǎng)格與有限元耦合方法[J]. 巖石力學(xué)與工程學(xué)報, 2007, 26(1): 75-80.

LI Shuchen, LI Shucai, SHU Bin, et al. Coupling meshless method and finite element method for seepage of jointed rock mass[J]. Chinese Journal of Rock Mechanics and Engineering,2007, 26(1): 75-80.

[4] Belystchko T, Krongauz, Organ Y D, et al. Meshless methods:An overview and recent developments[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 185(1): 3-47.

[5] 王兆清, 馮偉. 自然單元法研究進(jìn)展[J]. 力學(xué)進(jìn)展, 2004,34(4): 437-445.

WANG Zhao-qing, FENG Wei. Advances in natural element method[J]. Advances in Mechanics, 2004, 34(4): 437-445.

[6] Berrut J P, Trefethen L N. Barycentric Lagrange interpolation[J].SIAM Review, 2004, 46(3): 501-517.

[7] Nicholas J H. The numerical stability of barycentric Lagrange interpolation[J]. IMA Journal of Numerical Analysis, 2004, 24(4):547-556.

[8] 王兆清, 李淑萍, 唐炳濤. 一維重心型插值: 公式、算法和應(yīng)用[J]. 山東建筑大學(xué)學(xué)報, 2007, 22(5): 448-453.

WANG Zhaoqing, LI Shuping, TANG Bingtao. Formulations,algorithms and applications on barycentric interpolation in 1D[J].Journal of Shandong Jianzhu University, 2007, 22(5): 448-453.

[9] 王兆清, 李淑萍, 唐炳濤. 任意連續(xù)函數(shù)的多項式插值逼近[J]. 山東建筑大學(xué)學(xué)報, 2007, 22(2): 158-162.

WANG Zhaoqing, LI Shuping, TANG Bingtao. Polynomial interpolation approximations of arbitrary continuous functions[J].Journal of Shandong Jianzhu University, 2007, 22(2): 158-162.

[10] 王兆清, 李淑萍, 唐炳濤. 圓環(huán)變形及屈曲的重心插值配點法分析[J]. 機(jī)械強(qiáng)度, 2009, 31(2): 245-249.

WANG Zhaoqing, Li Shuping,Tang Bingtao. Deformation and buckling analysis of ring by barycentric interpolation collocation method[J]. Journal of Mechanical Strength, 2009, 31(2):245-249.

[11] 李淑萍, 王兆清. 非線性振動分析的重心插值配點法[J]. 噪聲與振動控制, 2008, 28(4): 49-52.

LI Shuping, WANG Zhaoqing. Barycentric Interpolation Collocation Method for Solving Nonlinear Vibration Problems[J].Noise and Vibration Control, 2008, 28(4): 49-52.

[12] 王兆清, 李淑萍, 唐炳濤, 等. 脈沖激勵振動問題的高精度數(shù)值分析[J]. 機(jī)械工程學(xué)報, 2009, 45(1): 288-292.

WANG Zhaoqing, LI Shuping, TANG Bingtao, et al. High precision numerical analysis of vibration problems under pulse excition[J]. Journal of Mechanical Engineering, 2009, 45(1):288-292.

[13] Hung K M, Ma C C. Theoretical analysis and digital photoelastic measurement of circular disks subjected to partially distributed compressions[J]. Experimental Mechanics, 2003, 43(2): 216-24.

[14] Ma C C, Hung K M. Exact full-field analysis of strain and displacement for circular disks subjected to partially distributed compressions[J]. International Journal of Mechanical Sciences,2008, 50(2): 275-292.

[15] Trefethen L N. Spectral methods in Matlab[M]. Philadelphia:SIAM, 2000: 211-215.

[16] Timoshenko S P, Goodier J N. Theory of elasticity[M]. New York: McGraw-Hill, 1951: 312-317.

[17] Atluri S N, Liu H T, Han Z D. Meshless local Petrov-Galerkin(MLPG) mixed collocation method for elasticity problems[J].CMES: Computer Modeling in Engineering & Sciences, 2006,14(3): 141-152.