极坐标系下弹性问题的重心插值配点法
李树忱1，王兆清2，袁超1
(1. 山东大学 岩土与结构工程研究中心，山东 济南，250061；
2. 山东建筑大学 工程力学研究所，山东 济南，250101 )
摘要：针对岩土工程中的孔洞及曲梁问题，提出一种在极坐标系下求解二维弹性问题的重心插值配点法。该方法分别在r和θ方向分别布置m和n个节点，生成求解区域上的节点。以一维重心Lagrange插值的张量积插值形式近似二维弹性问题的位移函数，代入位移表达的平衡方程和边界条件，平衡方程和边界条件分别在所有的计算节点和边界节点上精确成立，得到极坐标下弹性力学平衡方程和边界条件的离散代数表达式。利用一维重心Lagrange插值微分矩阵，将离散的平衡方程和边界条件表达为矩阵形式。利用置换法施加边界条件，求得在计算节点处的位移，进而通过微分矩阵直接求得计算节点处的应力。数值算例表明：极坐标下重心插值配点法具有计算格式简单、程序实施容易和计算精度高的特点。
关键词：弹性问题；极坐标系；重心Lagrange插值；微分矩阵；重心插值配点法；无网格方法
中图分类号：TU 443 文献标志码：A 文章编号：1672-7207(2013)05-2031-10
Barycentric interpolation collocation method for solving elastic problems
LI Shuchen1, WANG Zhaoqing2, YUAN Chao1
(1. Geotechnical and Structural Engineering Research Center, Shandong University, Jinan 250061, China;
2. Institute of Engineering Mechanics, Shandong Jianzhu University, Jinan 250101, China)
Abstract: Aiming at the hole and curved beam problem of geotechnical engineering, barycentric interpolation collocation method, a high accuracy meshless method for solving elastic problems in polar coordinate system was presented. The computational nodes were generated by tensor product from m and n nodes located in r and q directions, respectively. The functions of displacement were approximated by tensor product of barycentric Lagrange interpolation in one dimension (1D). The collocation method was used to obtain the systems of algebraic equations of governing equations and boundary conditions. The simple matrix form of the systems of algebraic equations was given using the notation of differential matrices from barycentric Lagrange interpolation in 1D. The boundary conditions were imposed by replacement method. The numerical results of displacements on nodes were obtained by solving the systems of algebraic equations of governing equations. The stresses values on nodes were directly computed using differential matrices. The numerical examples demonstrate that the proposed method has advantages of simple computational formulations, being easy to program and high precision.
Key words: elastic problem; polar coordinate system; barycentric Lagrange interpolation; differential matrix; barycentric interpolation collocation method; meshless method
弹性力学是研究岩石变形、断裂破坏问题的重要基础，对于复杂的岩石弹性力学问题一般很难得到问题的解析解，数值计算方法成为一种重要的研究手段。研究弹性力学问题的基本数值方法有有限元方法和有限差分法。有限元方法是一种基于低阶多项式插值的Galerkin方法，其计算精度依赖于单元大小和积分点数量的选取；有限差分法采用差分算子近似微分算子，其计算精度依赖于差分步长的大小。近年来兴起的各种类型无网格方法，在岩石力学领域得到广泛关注[1-3]。能够采用较少的计算代价，得到问题的高精度数值解，是数值分析领域研究者努力寻求的目标。基于移动最小二乘法的无单元Galerkin法尽管有不需要划分计算网格的优点，但在形成刚度矩阵过程中，需要进行大量的计算，并且不能直接施加边界条件[3-4]。基于自然邻点插值的自然单元法，边界条件施加方便，但在形成刚度矩阵过程中，需要寻找自然邻点，计算工作量巨大[5]。Lagrange插值是函数近似的重要手段，但当插值节点数量加大时，Lagrange插值表现为数值不稳定，Runge现象说明了这一问题，因此，工程计算中较少采用Lagrange插值作为近似未知函数的手段。将Lagrange插值公式改写为重心公式形式，得到重心Lagrange插值[6]。重心Lagrange插值具有很好的数值稳定性[7]和极高的插值近似精度[8-9]。采用重心Lagrange插值作为未知函数的近似函数的配点型方 法 —— 重心插值配点法，不需要划分计算网格，也不需要数值积分，已应用于求解一维工程边值问题[10]和动力学问题[11-12]中。现有的大部分数值方法是在直角坐标系下求解弹性力学问题。对于一些圆形区域上的工程问题，例如岩石力学的巴西圆盘实验，采用极坐标系分析非常方便[13-14]。在此，本文采用重心Lagrange插值的张量积形式作为二维问题未知函数的近似函数，基于极坐标系下弹性力学位移平衡方程的强形 式，提出求解极坐标系下弹性力学问题的重心插值配点法。
1 重心Lagrange插值及其微分矩阵
考虑定义在区间[a, b]上的函数u(x)，函数u(x)在节点a=x1＜x2＜…＜xn=b上的取值为ui=u(xi), i=1, 2, …, n。函数u(x)的重心Lagrange插值为[8]
(1)
其中：
；j=0, 1, …, n为重心插值权重。记重心Lagrange插值的插值基函数为
；k=1, 2, …, n (2)
函数u(x)的重心Lagrange插值可表示为
(3)
则函数u(x)的m阶导数可表示为
(4)
函数u(x)在节点x1, x2, …, xn处的m阶导数可表示为
；i=1, 2, …, n (5)
其中：
为第j个插值基函数在第i个节点处的m阶导数。将式(5)写成矩阵的形式：
(6)
其中：
为未知函数u(x)在节点处的m阶导数值列向量；矩阵
称为未知函数的m阶重心Lagrange插值微分矩阵，其元素为
；
为未知函数在节点处的函数值。一阶微分矩阵可以通过对式(2)求导直接得到，高阶微分矩阵可以通过下面递推公式得到[10-12]：
(7)
2 极坐标系下弹性问题的重心插值配点法公式
2.1 极坐标系下弹性力学的基本方程
取极坐标系
，某点
处的位移分量为
和
，应力分量为
和
和
，体力分量为
和
，则极坐标系下的平衡方程可表示为
(8)
位移-应变的几何关系为：
，
，
(9)
应力-应变本构关系为：
，
，
(10)
或者
，
，
(11)
将式(9)和(11)代入式(8)，得到位移表达的平面应力问题平衡方程：
(12)
位移表达的应力分量为
(13)
对于平面应变问题，需将式(10)~(12)中的弹性常数进行下列替换：
，
(14)
2.2 重心插值配点法计算公式
极坐标系下平面问题的位移解法需要求解偏微分方程(式(12))的边值问题。对于几何形状和外载荷不随θ改变的轴对称结构，其应力分布与θ无关，只是坐标r的函数，方程(12)将变为常微分方程。轴对称结构常微分方程边值问题的重心插值配点法求解方法，参见文献[10]。
在求解区域
的r方向布置m个计算节点r1, r2, …, rm，在θ方向布置n个计算节点θ1, θ2, …, θn，构成求解区域上的m×n个张量积计算节点(ri, θj)(其中，i=1, 2, …, m; j=1, 2, …, n)。采用节点r1, r 2, …, rm和θ1, θ2, …, θn上的重心Lagrange插值(式(1))的张量积形式插值近似径向位移ur和环向位移uθ：
(15)
其中：
和
分别为径向位移ur和环向位移uθ在计算节点
处的值；
；i=1, 2, …, m; j=1, 2, …, n；
, 
分别为在节点r1, r 2, …, rm和θ1, θ2, …, θn上的重心Lagrange插值基函数。
将式(15)代入位移表达的平衡方程(式(12))，得：
(16)
令式(16)在所有的计算节点
(k=1, 2, …, m; l=1, 2, …, n)成立，可得：
(17)
将径向位移、环向位移和体力，在计算节点, k=1, 2, …, m; l=1, 2, …, n处的值，排列成1个向量，表示为
(18)
其中：
；
；i=1, 2, …, m; j=1, 2, …, n。利用微分矩阵记号和矩阵张量积(Kronecker积)运算符号“
”[15]，微分方程(式(17))的矩阵表达式为：
(19)
其中：和
分别为在节点r1, r 2, …, rm和θ1, θ2, …, θn上的m阶重心Lagrange插值微分矩阵；Im和In分别为m阶和n阶单位矩阵；diag(1/r)，diag(1/r2)分别为函数1/r和1/r2在节点r1, r 2, …, rm处函数值构成的对角矩阵。引进记号：
(20)
将式(19)简记为：
(21)
则式(21)的解耦形式为
(22)
式(22)可简记为
(23)
按照平衡方程离散的同样方式，位移表达的应力分量(式(13))的重心插值离散形式表示为
(24)
由式(24)抽取边界节点的应力分量表达式，得到应力边界条件的重心插值离散计算公式。边界条件的施加可以采用置换法或附加法[10-12]。求解代数方程方程组得到计算节点的位移值后，采用式(24)可以方便地计算出所有节点的应力。
3 数值算例
数值算例采用MATLAB编写，以充分发挥MATLAB强大的矩阵运算能力。数值计算的绝对误差和相对误差分别定义为：
；
(25)
其中：
和
分别为数值计算解和解析解列向量；
表示向量的2范数。计算节点采用特殊分布的Chebyshev节点[15]。
算例1 曲梁的一般弯曲问题。考虑图1所示矩形等截面曲梁受端部载荷作用的平面应力问题。
曲梁的位移分量解析解为[16]
(26)
图1 端部受集中力作用的曲梁
Fig.1 Curved beam subjected to concentrated force on end
曲梁的应力分量解析解为
(27)
其中：
；
； 
；
；
；
；
；a＜r＜b；0＜θ＜
。
曲梁的边界条件的Saint-Venant提法为：
；
；
；
(28)
；
；
(29)
；
；
(30)
为得到问题的高精度解答，边界条件(式(29)和(30))按照应力解析解(式(27))和位移解析解(式(26))施加。
；
； (31)
；
；
(32)
数值计算时，曲梁的几何参数和力学参数分别取a=13，b=17，E=107，μ=0.25，P=104，单位一律采用国际单位制。不同数量Chebyshev节点重心Lagrange插值配点法计算的位移、应力误差见表1。由表1可知：采用较少数量的计算节点重心插值配点法计算的位移和应力具有很高的计算精度，并且随着计算节点数量的增加，位移和应力的计算精度随着提高；重心插值配点法计算的应力误差大于位移误差，这是由于计算应力需要数值求导造成计算精度损失。重心插值配点法计算的曲梁位移分量和应力分量分布分别见图2~6。
算例 2 含有孔洞无限大板的拉伸，应力集中因子问题。考虑如图7所示含有孔洞受单向拉伸的无限大板问题，孔洞的半径为a，板边均布拉力为P。设板的厚度为1，按照平面应力计算。
建立图7所示坐标系，极坐标系下板内应力分量的解析解为[17]：
表1 不同数量计算节点重心插值配点法计算的位移和应力误差
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 曲梁的环向位移分布
Fig.3 Distribution of hoop displacement in curved beam
图4 曲梁的径向正应力分布
Fig.4 Distribution of radial normal stress in curved beam
图5 曲梁的环向正应力分布
Fig.5 Distribution of hoop normal stress in curved beam
图6 曲梁的剪应力分布
Fig.6 Distribution of shear stress in curved beam
图7 含孔洞无限大板
Fig.7 A infinite plate with hole under uniaxial tension
(33)
位移解析解为：
(34)
其中：
为剪切模量；
；
为泊松比。
设想在孔洞外作1个半径为b的大圆，并且b ＞＞ a，大圆周和孔洞边界构成1个环形区域。由于对称性，取环形区域的1/4作为计算区域，如图8所示。采用位移法计算板的位移和应力分布时，r=a边界为自由孔洞边界，其应力边界条件为，。
图8 重心插值配点法的计算区域
Fig.8 Computational domain in barycentric interpolation collocation method
r=b边界上应力边界条件按照式(33)确定：
(35)
由于对称性，边界的边界条件为
，
；
，
。数值计算时，取a=0.5，b=10，E=107，μ=0.25，P=104，单位采用国际单位制。r方向取31个计算节点，θ方向取21个计算节点，计算节点采用Chebyshev型节点。重心插值配点法计算的位移和应力分量的计算误差如表2所示。从表2可知：重心插值配点法具有很高的计算精度。
最大环向应力
在孔洞边界，应力集中因子
的解析值等于3，重心插值配点法计算的应力集中因子等于3.000 2。
在r方向和θ方向取不同数量计算节点计算位移和应力的相对误差，448种不同节点数量组合和相对误差常用对数的关系分别见图9~13。由图9~13可以看出：随着计算节点数量的增加，位移和应力的计算误差显著降低；在相同的径向节点数量条件下，增加环向节点数量不能明显提高计算精度；在相同的环向节点数量下，增加径向节点数量可以显著提高计算精度；在环向采用10个计算节点、径向采用33个计算节点的计算精度可达10-6量级。
上述2个算例说明：对于几何区域为圆形、圆环或扇形的弹性力学问题，采用极坐标解答可以有效地简化问题的复杂性。对于含有圆形孔洞的无限大板拉伸问题，在板中截取1个与圆孔同心的大圆，将求解区域转化为圆环区域，利用直角坐标系和极坐标系下应力的转换关系，可以方便地得到圆环求解区域的应力边界条件。
对于一般形状或不规则形状的物理区域，不能直接采用重心插值配点法计算，可以将物理区域嵌入1个圆形或扇形的正则区域，在正则区域上采用极坐标系的重心插值配点法计算。计算得到正则区域上点的位移和应力值后，通过重心插值式(15)可以得到物理区域上任意一点的位移和应力。
表2 重心插值配点法计算的位移和应力误差
Table 2 Computational error of displacement and stress in barycentric interpolation collocation method
图9 径向位移相对误差与计算节点数量的关系
Fig.9 Relationship between relative error of radial displacement and number of computed nodes
图10 环向位移相对误差与计算节点数量的关系
Fig.10 Relationship between relative error of hoop displacement and number of computed nodes
图11 径向应力相对误差与计算节点数量的关系
Fig.11 Relationship between relative error of radial stress and number of computed nodes
图12 环向应力相对误差与计算节点数量的关系
Fig.12 Relationship between relative error of hoop stress and number of computed nodes
图13 剪应力相对误差与计算节点数量的关系
Fig.13 Relationship between relative error of shear stress and number of computed nodes
4 结论
(1) 采用重心Lagrange插值作为近似函数的重心插值配点法，在极坐标系下求解弹性力学问题，利用重心插值微分矩阵和矩阵张量积运算符号，可以将位移表达的极坐标系平衡方程的离散代数方程表示为简洁的矩阵形式。
(2) 将原来径向位移和环向位移耦合的偏微分方程组，通过矩阵操作实现解耦，这为计算程序的实施提供了极大地便利。微分矩阵的引进，提供了一种直接计算节点处应力值的简便快速方法。
(3) 重心插值配点法不需要划分计算网格，不需要数值积分，是一种真正意义上的无网格方法。
(4) 重心插值配点法具有计算格式简单、程序实施方便和计算精度高的特点。
参考文献：
[1] 卢波, 丁秀丽, 邬爱清. 无网格法对岩体不连续面的模拟[J]. 岩石力学与工程学报, 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插值无网格方法及其在力学边值问题中的应用[J]. 岩石力学与工程学报, 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] 李树忱, 李术才, 隋斌, 等. 节理岩体渗流的无网格与有限元耦合方法[J]. 岩石力学与工程学报, 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]. 力学进展, 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] 王兆清, 李淑萍, 唐炳涛. 一维重心型插值: 公式、算法和应用[J]. 山东建筑大学学报, 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] 王兆清, 李淑萍, 唐炳涛. 任意连续函数的多项式插值逼近[J]. 山东建筑大学学报, 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] 王兆清, 李淑萍, 唐炳涛. 圆环变形及屈曲的重心插值配点法分析[J]. 机械强度, 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] 王兆清, 李淑萍, 唐炳涛, 等. 脉冲激励振动问题的高精度数值分析[J]. 机械工程学报, 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.
(编辑 邓履翔)
收稿日期：2012-04-01；修回日期：2012-07-11
基金项目：国家重点基础研究发展计划(“973”计划)项目(2010CB732002)；国家自然科学基金面上资助项目(51179098)；山东大学自主创新基金资助项目(2010TS038)；山东建筑大学研究生教育创新计划项目(YC1005)资助项目
通信作者：王兆清(1965-)，男，山东省荣成人，博士，副教授，从事工程数值分析方法领域的教学和研究；电话：0531-88399180；Email: sdjzuwang@126.com
