引 言

近年来发展迅速的无网格方法[1-7]可以减少或消除网格划分给数值计算带来的困难，在自适应分析、大变形和裂纹扩展等问题领域比传统的有限元法更加灵活和有效．在现有的无网格方法中，Liu等[8]于1995年提出的重构核粒子法是目前应用最为广泛的无网格方法之一．其核心思想是在光滑粒子法[9]中引入核函数修正项，使理论上精确重构有限域的近似函数成为可能．然而，重构核粒子法的形函数不满足插值特性．重构核插值法[10-11]通过简单函数引入插值特性，并利用增强函数构造重构条件，得到一个具有任意离散点插值特性的形函数．重构核插值法不仅免去了处理本质边界条件的困难，而且其形函数能够精确重构插值点多项式的真值、近似精度高．

在很多实际工程中，有一类几何形状、约束条件与外载荷均对称于某一固定轴的轴对称问题．这类构件在任一对称面内相对应点上的位移、应变以及应力大小相等且方向对称，因此可将三维空间问题转换为二维问题，从而能快速有效地求解．目前，复杂轴对称问题求解的数值方法研究已引起了许多学者[12-14]的极大关注．为充分发挥重构核插值法的优势，本文将其与弹性动力学空间轴对称问题相结合，利用虚位移原理详细推导了相应的离散方程，并编制了相应的计算程序．最后，通过典型算例的计算和对比,分析验证了该方法的有效性和合理性．

1 重构核插值法的形函数

在轴对称面Ω内，任一点x处的函数u(x)可近似为[10-11]

(1)

式中N为计算点x邻域内的节点数，ψI(x)为重构核插值法的形函数，且有

(2)

由粒子xI处的核函数及其在粒子xI的值来构造，取具有Kronecker delta函数特性的简单形式，由扩展的重构条件构造，即

(3)

(4)

式中HT(x-xI)为m维基函数向量，且

HT(x-xI)=[1, x-xI, y-yI, …, (y-yI)m],

(5)

(6)

(7)

为便于形函数偏导数的计算，将式(4)改写为

(8)

式中

AT(x)=HT(x-xI)Q-1(x)，

(9)

(10)

因此，的偏导数为

(11)

若取支持域为圆形，中心为xI，则两个核函数分别为

(12)

(13)

式中Φ为四次样条函数，和本文分别取为和和分别为节点xI到距其最近的第1个和第9个节点之间的距离．

2 轴对称弹性动力学问题

在轴对称面Ω上，位移u=[ur,uz]T、加速度应力σ=[σr，σθ，σz,τrz]T和应变ε=[εr，εθ，εz，γrz]T的关系方程为

(14a)

ε=Lu,

(14b)

σ=Dε,

(14c)

此外还应满足相应的边界条件和初始条件．上式中，ρ为质量密度，b=[br,bz]T为体力向量， L为微分算子矩阵，D为弹性矩阵，且

(15)

(16)

式中

(17)

其中E和ν分别为材料的弹性模量和Poisson比．

由于只对空间域进行离散，轴对称面Ω内的试函数u(x,t)可由式(1)表示为

(18)

由式(14b)、(15)和(18)可得应变ε为

(19)

式中

(20)

对于轴对称弹性动力学问题，平衡方程对应的标准变分弱形式为

(21)

式中为给定的面力向量．

将式(18)和(19)代入式(21)，并注意到节点位移变分δu的任意性，最终可得控制方程的离散形式为

(22)

式中M,K和f(t)分别是系统的质量矩阵、刚度矩阵和节点载荷向量．它们的各元素可具体表示为

(23a)

(23b)

(23c)

其中

(24)

3 时间积分方案

式(22)是二阶常微分方程组，其求解主要有直接积分法和振型叠加法两类．文献[15]对直接积分法和振型叠加法的实用范围进行了讨论，本文采用直接积分法中的Newmark-β法．

在t～t+Δt的时间区域内，Newmark-β积分方法采用下列的假设，即

(25)

(26)

式中α和δ代表Newmark-β积分参数，本文取α=1/4和δ=1/2．将式(25)和(26)代入式(22)可得

(27)

式中

(28)

(29)

4 数 值 算 例

4.1 受突加内压的厚壁圆筒

受突加内压力作用的无限长厚壁圆筒如图1所示．内半径r2=0.1 m，外半径r1=0.2 m，弹性模量E=2.1×1011 Pa，Poisson比ν=0.3，质量密度ρ=7.85×103 kg/m3,内压力p=75 MPa．截取长L=0.1 m的厚壁圆筒作为计算模型，并在其轴对称面上布置16×5个计算节点，如图2所示．为了进行数值积分，轴对称面上划分15×4个积分背景网格，且每一网格中采用3×3个Gauss点．图3为时间歩长取Δt=4.0×10-7 s时计算得到的内表面径向位移随时间的变化曲线．为了进行对比，图3还给出了ABAQUS和重构核粒子法的计算结果．显然，本文方法的计算精度要好于重构核粒子法，与ABAQUS的结果吻合得更好．

4.2 受突加内压的厚壁球壳

考虑一受突加内压力p=75 MPa作用的厚壁球壳．球壳内壁半径r2=0.16 m，外壁半径r1=0.2 m，弹性模量E=2.1×1011 Pa，Poisson比ν=0.3，密度ρ=7.85×103 kg/m3．此问题是球对称问题，现将它作为轴对称问题进行求解，相应的节点布置方案如图4所示．为了进行数值积分，轴对称面上划分6×18个积分背景网格，且每一网格中采用3×3个Gauss点．计算总时间为3.0×10-4 s，时间歩长取为Δt=3.0×10-7 s．图5给出了本文方法和ABAQUS软件计算得到的球壳内壁的径向位移随时间的变化曲线．从图5可以看出，两者吻合得很好，表明本文所推导的轴对称弹性动力学问题的重构核插值法及其离散形式是正确的，此方法是可行的．

5 结 论

重构核插值法作为一种具有代表性的无网格方法，由于具有易于建立高阶近似函数和形函数满足Kronecker delta函数特性等优势而被计算力学界广泛关注．本文将重构核插值法进一歩推广应用于轴对称弹性动力学问题．基于虚位移原理，本文详细推导了轴对称弹性动力学问题的离散方程，并编制了相应的计算程序．最后，数值算例的分析表明，本文数值解与有限元软件ABAQUS的计算结果吻合得很好，说明本文提出的轴对称弹性动力学问题的重构核插值法及其离散形式是正确的，并且具有前处理方便、精度高和可直接施加本质边界条件等优点．

参考文献

[1] LI X L, LI S L. A meshless projection iterative method for nonlinear Signorini problems using the moving Kriging interpolation[J]. Engineering Analysis With Boundary Elements, 2019, 98: 243-252.

[2] 杨建军, 郑健龙. 无网格局部强弱法求解不规则域问题[J]. 力学学报, 2017, 49(3): 659-666.(YANG Jianjun, ZHENG Jianlong. Meshless local strong-weak (MLSW) method for irregular domain problems[J]. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(3): 659-666.(in Chinese))

[3] 肖毅华, 张浩锋, 平学成. 无网格对称粒子法中两类热边界条件的处理[J]. 华东交通大学学报, 2014, 31(4): 65-70.(XIAO Yihua, ZHANG Haofeng, PING Xuecheng. Treatment of two kinds of thermal boundary conditions in meshless symmetric particle method[J]. Journal of East China Jiaotong University, 2014, 31(4): 65-70.(in Chinese))

[4] FU Z J, XI Q, CHEN W, et al. A boundary-type meshless solver for transient heat conduction analysis of slender functionally graded materials with exponential variations[J]. Computers and Mathematics With Applications, 2018, 76(4): 760-773.

[5] LI M, DOU F F, KORAKIANITIS T, et al. Boundary node Petrov-Galerkin method in solid structures[J]. Computational and Applied Mathematics, 2018, 37(1): 135-159.

[6] 王峰, 周宜红, 郑保敬, 等. 基于滑动Kriging插值的MLPG法求解结构非耦合热应力问题[J]. 应用数学和力学, 2016, 37(11): 1217-1227.(WANG Feng, ZHOU Yihong, ZHENG Baojing, et al. A meshless local Petrov-Galerkin method based on the moving Kriging interpolation for structural uncoupled thermal stress analysis[J]. Applied Mathematics and Mechanics, 2016, 37(11): 1217-1227.(in Chinese))

[7] 孙新志, 李小林. 复变量移动最小二乘近似在Sobolev空间中的误差估计[J]. 应用数学和力学, 2016, 37(4): 416-425.(SUN Xinzhi, LI Xiaolin. Error estimates for the complex variable moving least square approximation in Sobolev spaces[J]. Applied Mathematics and Mechanics, 2016, 37(4): 416-425.(in Chinese))

[8] LIU W K, JUN S, ZHANG Y F. Reproducing kernel particle methods[J]. International Journal for Numerical Methods in Engineering, 1995, 20(8/9): 1081-1106.

[9] GINGOLD R A, MONAGHAN J J. Smoothed particle hydrodynamics: theory and applications to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society, 1977, 181(3): 275-389.

[10] CHEN J S, HAN W, YOU Y, et al. A reproducing kernel method with nodal interpolation property[J]. International Journal for Numerical Methods in Engineering, 2003, 56(7): 935-960.

[11] 李中华, 秦义校, 崔小朝. 弹性力学的插值型重构核粒子法[J]. 物理学报, 2012, 61(8): 25-31.(LI Zhonghua, QIN Yixiao, CUI Xiaochao. Interpolating reproducing kernel particle method for elastic mechanics[J]. Acta Physica Sinica, 2012, 61(8): 25-31.(in Chinese))

[12] 韩治, 杨海天, 王斌. 无网格伽辽金法求解轴对称问题[J]. 工程力学, 2005, 22(5): 64-68.(HAN Zhi, YANG Haitian, WANG Bin. Solving axisymmetric problems via EFGM[J]. Engineering Mechanics, 2005, 22(5): 64-68.(in Chinese))

[13] 陈建桥, 梁元博, 丁亮. 无网格局部Petrov-Galerkin法求解轴对称问题[J]. 华中科技大学学报(城市科学版), 2007, 24(4): 9-12.(CHEN Jianqiao, LIANG Yuanbo, DING Liang. Numerical analysis of axisymmetric problems by MLPG[J]. Journal of Huazhong University of Science and Technology(Urban Science Edition), 2007, 24(4): 9-12.(in Chinese))

[14] 陈莘莘, 李庆华, 刘永胜. 轴对称动力学问题的无网格自然邻接点Petrov-Galerkin法[J]. 振动与冲击, 2015, 34(3): 61-65.(CHEN Shenshen, LI Qinghua, LIU Yongsheng. Meshless natural neighbour Petrov-Galerkin method for axisymmetric dynamic problems[J]. Journal of Vibration and Shock, 2015, 34(3): 61-65.(in Chinese))

[15] 姜清辉, 周创兵, 漆祖芳. 基于Newmark积分方案的DDA方法[J]. 岩石力学与工程学报, 2009, 28(1): 2778-2783.(JIANG Qinghui, ZHOU Chuangbing, QI Zufang. Discontinuous deformation analysis method based on Newmark integration algorithm[J]. Chinese Journal of Rock Mechanics and Engineering, 2009, 28(1): 2778-2783.(in Chinese))

A Reproducing Kernel Interpolation Method for Axisymmetric Elastodynamic Problems

CHEN Shenshen, ZENG Jiawei

(School of Civil Engineering and Architecture, East China Jiaotong University, Nanchang 330013, P.R.China)

Abstract: The reproducing kernel interpolation method (RKIM) is a novel type of meshless method emerging in recent years. Because the shape functions of the RKIM have point interpolation property and high-order smoothness, the essential boundary conditions can be imposed directly and high computational accuracy is ensured as well. In order to solve the elastodynamic problems for 3D axisymmetric solids more effectively, a novel numerical method based on the RKIM was presented and discussed. Due to axial symmetry of geometry and boundary conditions, only a set of discrete nodes on a cross section are required in the computation and therefore the preprocessing of this method is very simple. The Newmark-β algorithm was employed for time integration. Numerical examples show that, the proposed method for solving axisymmetric elastodynamic problems possesses the advantages of meshless methods and high accuracy.

Key words: elastodynamics; axisymmetric; meshless method; reproducing kernel interpolation method

收稿日期： 2018-09-13； 修订日期： 2018-11-15

基金项目： 国家自然科学基金(11462006;11772129)

作者简介： 陈莘莘(1975—)，男，教授，博士(通讯作者. E-mail: chenshenshen@tsinghua.org.cn).

中图分类号： O241; O343

文献标志码： A

文章编号：1000-0887(2019)08-0938-07

ⓒ 应用数学和力学编委会，ISSN 1000-0887

DOI: 10.21656/1000-0887.390242

Foundation item: The National Natural Science Foundation of China(11462006;11772129)

引用本文/Cite this paper: 陈莘莘， 曾佳伟. 轴对称弹性动力学问题的重构核插值法[J]. 应用数学和力学, 2019, 40(8): 938-944.CHEN Shenshen, ZENG Jiawei. A reproducing kernel interpolation method for axisymmetric elastodynamic problems[J]. Applied Mathematics and Mechanics, 2019, 40(8): 938-944.