弹性力学和弹性断裂问题的准凸重构核粒子法

弹性力学和弹性断裂问题的准凸重构核粒子法 Quasi-Convex Reproducing Kernel Particle Method for Elasticity and Elastic Fracture Problems first-author 田利瑞，男，硕士生，现从事计算力学方面的研究.E-mail：tlr923@163.com 田利瑞，男，硕士生，现从事计算力学方面的研究.E-mail： tlr923@163.com tlr923@163.com corresp 孔令好, domili@whut.edu.cn 孔令好, domili@whut.edu.cn domili@whut.edu.cn 武汉理工大学 土木工程与建筑学院，湖北 武汉 430070 武汉理工大学 土木工程与建筑学院，湖北 武汉 430070 School of Civil Engineering and Architecture, Wuhan University of Technology, Wuhan 430070, Hubei, China School of Civil Engineering and Architecture, Wuhan University of Technology, Wuhan 430070, Hubei, China 通过放宽重构条件, 准凸重构核近似可以构造任意高阶正性更好的形函数.本文结合控制方程的Galerkin弱式，将准凸重构核近似应用于弹性力学和弹性断裂问题的计算，推导了相关计算公式，建立了弹性力学和弹性断裂问题的准凸重构核粒子法(quasi-convex reproducing kernel particle method, QCRKPM).通过编写QCRKPM的MATLAB程序，对中心圆孔方板和边缘裂纹方板这两个算例进行了模拟计算.将QCRKPM的数值精度和计算效率与无单元Galerkin(element-free Galerkin, EFG)方法和重构核粒子法(reproducing kernel particle method，RKPM)进行对比，验证了QCRKPM的高效性和高精度. By applying the relaxed reproducing conditions, the so-called quasi-convex reproducing kernel approximation gains the ability to construct arbitrary order shape functions with enhanced positivity. In this paper, combined with the Galerkin weak form of the control equations, the quasi-convex reproducing kernel approximation was applied to analyze elasticity and elastic fracture problems. The theoretical formulas were derived, and the quasi-convex reproducing kernel particle method (QCRKPM) for elasticity and elastic fracture problems were proposed. Two examples including a square plate with a central circular hole and a square plate with an edge crack were simulated by MATLAB. By comparing the numerical accuracy and computational efficiency of the proposed QCRKPM with that of the element-free Galerkin (EFG) method and reproducing kernel particle method (RKPM), the high efficiency and precision of the QCRKPM were verified. 准凸重构核近似 准凸重构核粒子法 重构核粒子法 弹性力学问题 弹性断裂问题 quasi-convex reproducing kernel approximation quasi-convex reproducing kernel particle method reproducing kernel particle method elasticity problem elastic fracture problem 引言 1 无网格方法的形函数构造完全基于节点，降低了基于单元的传统数值方法在单元生成和计算过程中网格重划等方面的计算量，简化了数据处理过程，提高了数值精度，能够有效处理有限元法等传统数值方法不擅长的大变形和断裂问题 [ 1 1 ~ 6 6 1~6 ] .其中，重构核粒子法 [ 7 7 , 8 8 ] (reproducing kernel particle method，RKPM)是目前应用最为广泛的无网格方法之一. RKPM是在光滑粒子法 [ 9 9 ] 的基础上发展起来的，其核心是构造场变量逼近函数的重构核近似 [ 7 7 , 8 8 ] .通过引入重构条件，重构核近似获得了比传统核近似更好的光滑性和更高的精度，并完全消除了光滑粒子法的不稳定性.但是，重构核近似构造的形函数不是非负函数，即不具备凸近似性质.为了构造正性更好的形函数，Wang等 [ 10 10 ] 放宽了传统重构核近似的重构条件，提出了具备准凸近似性质的准凸重构核近似(quasi-convex reproducing kernel approximation)，使其能够在保证数值精度以及计算效率的情况下，构造任意高阶正性更好的形函数.文献[ 10 10 ]将准凸重构核近似应用于弹性问题，数值结果表明准凸性质可以提高重构核近似的计算精度.文献[ 10 10 ]还指出，在进行动力分析时，准凸重构核近似能够构造正性更好的质量矩阵，对计算精度的提高将更为显著.目前，准凸重构核粒子法(quasi-convex reproducing kernel particle method, QCRKPM)已被应用于求解各向异性弹性力学问题 [ 11 11 ] 、界面问题 [ 12 12 ] 和带裂纹板的自由振动问题 [ 13 13 ] 等. 本文将准凸重构核粒子法应用于弹性力学和弹性断裂问题，采用准凸重构核近似构造位移场的逼近函数，结合弹性力学基本方程的等效Galerkin弱式，采用罚函数法施加本质边界条件，推导了相关计算公式，建立了弹性力学和弹性断裂问题的准凸重构核粒子法，编写相应的MATLAB程序，对中心圆孔方板和边缘裂纹方板这两个算例进行计算分析.将本文中的准凸重构核粒子法的数值结果与相应解析解和分别采用无单元Galerkin(EFG)方法和RKPM的数值解进行对比，验证本文方法的高效性和高精度. 准凸重构核近似 1 在不增加计算量和保证数值精度的前提下，Wang等 [ 10 10 ] 通过放宽传统重构核近似的重构条件，提出了具有准凸近似特性的准凸重构核近似. 在问题域 Ω 上满布一组总数为 N 的节点{ x I } I =1 N ， x I 为点 x 支持域内第 I 个节点.域内任意点 x 处场变量 u ( x )的逼近函数 u h ( x )可表示为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458739&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458739&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458739&type=middle 其中， N P 为点 x 支持域内节点的总数， u I 为节点 x I 处的场变量参数. Φ I ( x )为点 x 的形函数在节点 x I 处的值，其表达式为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458743&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458743&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458743&type=middle 其中， q ( x I - x , ω I )为准凸重构核近似的基函数列向量，T为转置符号， ω I = β h 2 为节点松弛参数， β 为松弛系数， h 为域内节点间距； b ( x )为待定系数列向量； w ( x I - x )为权函数；Δ V I 为节点区域度量.对于二维问题，节点区域度量为面积，本文采用Voronoi图 [ 11 11 ] 来确定. 准凸重构核近似的基函数 q ( x I - x , ω I )可由标准多项式基函数 p ( x I - x )得到 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458746&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458746&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458746&type=middle 其中， q 0 ( x I - x , ω I )= p 0 ( x I - x )； q 1 ( x I - x , ω I )= p 1 ( x I - x )； q i ( x I - x , ω I )= p i ( x I - x )- g i ( x I - x , ω I ), 2≤ i ≤ n ， g i ( x I - x , ω I )为准凸重构核近似的节点松弛函数，满足 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458749&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458749&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458749&type=middle 式中，C i i -2 = i !/[( i -2)!2!]. 待定系数向量 b ( x )可通过如下放宽后的重构条件来确定 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458753&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458753&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458753&type=middle 其中， p ( 0 )=(1, 0, …, 0) T . 由(2)式和(5)式可知 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458758&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458758&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458758&type=middle http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458762&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458762&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458762&type=middle 由(2)式可知，点 x 的准凸重构核近似形函数在其支持域内节点 x I 处的值为： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458765&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458765&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458765&type=middle 对于二维问题，准凸重构核近似的基函数如下： 线性基： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458768&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458768&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458768&type=middle 二次基： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458771&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458771&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458771&type=middle 式中，( x , y )和( x I , y I )分别为点 x 和其支持域内节点 x I 在平面直角坐标系下的坐标， ω x I 和 ω y I 分别为平面直角坐标系 x 方向和 y 方向的节点松弛函数, 其计算式分别为： ω x I = βh x 2 , ω y I = βh y 2 , 其中 h x 和 h y 分别为平面直角坐标系 x 方向和 y 方向的节点间距.从(3)式和(9)式可以看出，对于二维问题，准凸重构核近似和重构核近似的线性基函数(即线性标准多项式基函数)完全相同，由(8)式可知，相应形函数也完全相同. 准凸重构核粒子法 1 对于线弹性静力问题，考察带力边界 Γ t 和位移边界 Γ u 的二维问题域 Ω ，利用应力-应变关系和应变-位移关系，采用罚函数法 [ 14 14 ] 施加位移边界条件，可得控制方程的Galerkin弱式 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458775&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458775&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458775&type=middle 式中， δ 为变分符号， L 为微分算子矩阵， D 为弹性本构矩阵， b 和 t 分别为体力和面力向量， α 为惩罚因子， u 为节点位移向量， u 为边界位移向量. 将位移场的准凸重构核近似代入(11)式中，得到线弹性静力问题的总体离散系统方程 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458779&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458779&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458779&type=middle 式中： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458783&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458783&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458783&type=middle 其中， http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458790&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458790&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458790&type=middle 对于线弹性断裂问题，基本方程与线弹性静力问题一样，惟一的区别在于线弹性断裂分析在计算出应力解后还需要进一步计算应力强度因子. 数值算例 1 为了验证本文弹性力学和弹性断裂问题的准凸重构核粒子法的正确性和有效性，采用MATLAB(R2012a)编制了计算程序，分别对中心圆孔方板和边缘裂纹方板这两个算例进行计算分析.同时，将本文方法的数值解和解析解与由EFG方法 [ 15 15 ] 和RKPM [ 7 7 , 8 8 ] 得到的数值解进行对比验证.由于准凸重构核近似和重构核近似的一次基函数和相应形函数完全相同，因此以下算例中的3种方法均采用二次基函数.采用三次样条权函数和矩形支持域，每个积分单元采用4×4阶Gauss积分，按平面应力进行计算. 算例1为受单向拉伸作用的中心圆孔方板.该算例中，背景积分网格为16×13，计算精度采用系统能量范数误差来衡量.数值解的系统能量范数误差越小，计算精度越高，反之越低.系统能量范数误差 e e 采用下式计算 [ 11 11 ] ： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458794&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458794&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458794&type=middle 式中， ε num 和 ε exact 分别为应变数值解和解析解. 算例2为边界上施加Ⅰ型裂纹位移场的方板.该算例中背景积分网格为20×20，计算精度以Ⅰ型裂纹应力强度因子数值解的相对误差来衡量.应力强度因子数值解的相对误差越小，计算精度越高，反之越低.应力强度因子解析解为 [ 16 16 ] ： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458799&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458799&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458799&type=middle 式中， http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458880&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458880&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458880&type=middle http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458884&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458884&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458884&type=middle 为应力， l (m)为裂纹长度， L (m)为方板边长. 算例1 2 位移解与应力解 3 考虑受单向拉伸作用的中心圆孔方板的对称性，只取1/4结构进行计算，如 图 1 图 1 所示.圆孔方板的相关几何尺寸为 a =1 m， b =6 m，材料的弹性模量 E =3.0×10 5 Pa，泊松比 ν =0.25，均布荷载为 p =10 N/m，不计自重. 对于中心圆孔有限平板，当平板宽度远远大于圆孔半径时，径向位移 u r ( r , θ )，环向位移 u θ ( r , θ )，径向正应力 σ r ( r , θ )，环向正应力 σ θ ( r , θ )以及剪应力 τ rθ ( r , θ )的解析解 [ 17 17 ] 如下： http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458808&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458808&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458808&type=middle 其中， r 为待求点到坐标原点 O 的距离， θ 为待求点与坐标原点 O 连线所在直线与 x 轴的夹角.水平位移 u x 和 x 方向正应力 σ xx 可根据弹性力学的坐标转换公式，由(15)式推导得到. 如 图 2 图 2 所示，将求解域离散为15×14个场节点，采用Voronoi图 [ 11 11 ] 计算各场节点对应的Δ V I .控制影响域大小的无量纲因子 d max =4.0，惩罚因子 α 的值取1×10 4 E .参考文献[ 10 10 ]，松弛系数 β 取0.6. 采用EFG方法 [ 15 15 ] 、RKPM [ 7 7 , 8 8 ] 、QCRKPM [ 10 10 ] 计算得到1/4中心圆孔方板 y =0截面节点的水平位移 u x 和 x 方向正应力 σ xx 的数值解，如 图 3 图 3 和 图 4 图 4 所示. 表 1 表 1 给出了3种方法计算得到的1/4中心圆孔方板 y =0截面节点的水平位移相对误差.由 图 3 图 3 及 表 1 表 1 可知，3种方法计算得到的1/4中心圆孔方板 y =0截面节点的水平位移解基本重合，具有较高的精度，相对误差均小于0.05%. 表 2 表 2 列出了3种方法的形函数计算时间.由 表 2 表 2 可知，EFG方法的形函数计算时间最短，计算效率最高，QCRKPM与RKPM的形函数计算时间接近且QCRKPM略短. 表 3 表 3 给出了3种方法计算得到的1/4中心圆孔方板 y =0截面节点的 x 方向正应力相对误差.由 图 4 图 4 和 表 3 表 3 可知，3种方法得到的方板 y =0截面节点的 x 向正应力 σ xx 解也基本重合，且精度较高，相对误差均小于7%. 4 % x /m EFG RKPM QCRKPM 1.000 000 0.013 062 0.037 128 0.040 738 1.384 615 0.015 174 0.011 001 0.010 773 1.769 231 0.000 120 0.001 471 0.000 928 2.153 846 0.002 979 0.000 462 0.000 583 2.538 462 0.001 675 0.000 907 0.001 335 2.923 077 0.002 624 0.000 545 0.001 270 3.307 692 0.001 596 0.000 784 0.001 451 3.692 308 0.000 415 0.000 271 0.000 312 4.076 923 0.000 711 0.001 268 0.001 032 4.461 538 0.003 459 0.001 930 0.001 911 4.846 154 0.001 281 0.001 700 0.001 593 5.230 769 0.001 490 0.001 902 0.002 268 5.615 385 0.000 420 0.000 263 0.000 071 6.000 000 0.008 071 0.006 639 0.008 718 方法 形函数计算时间/s EFG 0.165 42 RKPM 0.628 90 QCRKPM 0.616 29 4 % x /m EFG RKPM QCRKPM 1.000 000 0.629 168 1.212 668 1.317 205 1.384 615 6.106 910 0.042 266 3.341 068 1.769 231 6.622 388 3.673 654 3.206 940 2.153 846 0.026 409 1.768 002 2.043 297 2.538 462 0.205 487 0.553 011 0.496 765 2.923 077 0.066 833 1.409 502 1.271 846 3.307 692 1.336 021 1.527 615 1.560 863 3.692 308 2.542 100 1.811 785 1.855 989 4.076 923 0.958 320 0.309 754 0.310 792 4.461 538 1.678 834 0.714 014 0.660 921 4.846 154 2.143 748 2.745 789 2.716 082 5.230 769 2.090 766 1.438 862 1.434 171 5.615 385 2.047 951 1.871 291 1.910 288 6.000 000 3.742 573 3.736 974 4.393 766 计算参数分析 3 采用EFG方法、RKPM和QCRKPM，计算不同影响域放大因子、惩罚因子、环向节点数和径向节点数情况下的系统能量范数误差和形函数计算时间. 图 5 图 5 及 表 4 表 4 分别给出了节点分布为15×14，惩罚因子 α 取1×10 4 E ，影响域放大因子(无量纲) d max 在2.5~8范围内均匀选取12个值(间隔为0.5)，采用3种方法计算得到的系统能量范数误差和形函数计算时间.由 图 5 图 5 可知，EFG方法的系统能量范数误差随 d max 的增大而波动增大，RKPM及QCRKPM的系统能量范数误差随 d max 的增大而逐渐增大.QCRKPM和RKPM的系统能量范数误差大小相近，且均明显小于EFG方法.由 表 4 表 4 可知，由于影响域内节点数量越多，形函数计算时间越长，3种方法的形函数计算时间随 d max 的增大而不断增大，QCRKPM和RKPM的形函数计算时间相近且均大于EFG方法. 4 s d max EFG RKPM QCRKPM 2.5 0.097 84 0.336 04 0.334 09 3 0.124 99 0.420 71 0.422 12 3.5 0.147 16 0.517 43 0.515 92 4 0.169 83 0.608 84 0.612 44 4.5 0.191 47 0.704 79 0.702 80 5 0.215 05 0.803 89 0.803 15 5.5 0.238 04 0.890 96 0.891 23 6 0.255 86 0.982 22 0.971 35 6.5 0.272 62 1.036 41 1.035 23 7 0.290 52 1.109 89 1.116 76 7.5 0.303 53 1.157 14 1.161 47 8 0.313 57 1.204 68 1.209 75 图 6 图 6 及 表 5 表 5 分别给出了节点分布为15×14， d max 取4.0，惩罚因子 α 分别取10 E 、1×10 2 E 、1×10 3 E 、1×10 4 E 、1×10 5 E 时，采用3种方法计算得到的系统能量范数误差和形函数计算时间.由 图 6 图 6 可知，随惩罚因子的增大，由上述3种方法计算得到的系统能量范数误差均先减小后增大，当惩罚因子取1×10 2 E 时，各方法的系统能量范数误差均最小.由 表 5 表 5 可知，随惩罚因子的增大，3种方法的形函数计算时间变化较小，QCRKPM的形函数计算时间略少于RKPM，但两者均大于EFG方法. 4 s lg( α /E) EFG RKPM QCRKPM 1 0.161 47 0.662 13 0.657 37 2 0.160 11 0.632 58 0.612 03 3 0.159 31 0.615 43 0.611 96 4 0.156 88 0.615 52 0.613 85 5 0.158 07 0.615 61 0.611 27 图 7 图 7 及 表 6 表 6 分别给出了 d max 取4.0，惩罚因子 α 取1×10 4 E ，径向节点均取14，环向节点分别为11、13、15、17、19时，采用3种方法计算得到的系统能量范数误差和形函数计算时间.由 图 7 图 7 可知，当径向节点数为14时，EFG方法的系统能量范数误差随着环向节点数的增大呈波动变化，而RKPM和QCRKPM的系统能量范数误差随着环向节点数的增大而逐渐减小.由 表 6 表 6 可知，3种方法的形函数计算时间随环向节点数的增大而逐渐增大，QCRKPM和RKPM的形函数计算时间相近且均大于EFG方法. 4 s 环向节点数 EFG RKPM QCRKPM 11 0.127 19 0.472 71 0.498 67 13 0.146 11 0.531 95 0.526 27 15 0.150 50 0.581 61 0.576 92 17 0.164 19 0.629 62 0.623 24 19 0.180 15 0.674 91 0.677 14 图 8 图 8 及 表 7 表 7 分别给出了 d max 取4.0，惩罚因子 α 取值为1×10 4 E ，环向节点均取15，径向节点分别为12、13、14、15、16时，采用3种方法计算得到的系统能量范数误差和形函数计算时间.由 图 8 图 8 可知，当环向节点数为15时，RKPM和QCRKPM的系统能量范数误差随径向节点数的增大而逐渐减小，而EFG方法的系统能量范数误差随径向节点数的增大波动变化.由 表 7 表 7 可知，3种方法的形函数计算时间随径向节点数的增大而逐渐增大，QCRKPM和RKPM的形函数计算时间相近且均大于EFG方法. 4 s 径向节点数 EFG RKPM QCRKPM 12 0.112 83 0.422 14 0.422 36 13 0.133 37 0.502 56 0.525 42 14 0.150 70 0.577 90 0.576 80 15 0.172 02 0.664 80 0.680 23 16 0.195 00 0.755 68 0.756 72 从 图 5 图 5 ~ 8 8 可以明显地观察到，QCRKPM和RKPM的系统能量范数误差均显著小于EFG方法.除 图 5 图 5 ( 图 5 图 5 中QCRKPM和RKPM的系统能量范数误差相近)外，QCRKPM的系统能量范数误差均最小，表明其计算精度最高. 随机布点 3 问题域内采用随机布点方式，4种随机布点方案的节点总数分别为182、210、238、266，如 图 9 图 9 所示. 图 10 图 10 及 表 8 表 8 分别给出了 d max 取4.0，惩罚因子 α 为1×10 4 E 时，3种方法在4种随机布点方式下的系统能量范数误差及形函数计算时间.由 图 10 图 10 可知，3种方法的系统能量范数误差都随随机节点总数的增大而逐渐降低，且与RKPM和EFG方法相比，QCRKPM的系统能量范数误差最小，计算精度最高.由 表 8 表 8 可知，3种方法的形函数计算时间均随随机节点总数的增大而逐渐增大，QCRKPM的形函数计算时间比RKPM略小，但两者均大于EFG方法. 4 s 随机节点总数 EFG RKPM QCRKPM 182 0.136 70 0.527 51 0.516 25 210 0.153 27 0.577 93 0.564 57 238 0.159 90 0.632 37 0.592 35 266 0.170 01 0.673 54 0.623 98 算例2 2 图 11 图 11 所示为一个长 L =16 mm，宽 H =16 mm的边缘裂纹方板，取单位厚度，裂纹长度 l =8 mm.在其外部边界施加Ⅰ型裂纹 [ 16 16 ] 对应的线性位移场(由(16)式和(17)式计算得到).材料的弹性模量 E =2×10 5 MPa，泊松比 μ =0.286，按平面应力进行计算.裂纹尖端场的 x 方向位移解析解 u 为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458851&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458851&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458851&type=middle y 方向位移解析解 v 为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458854&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458854&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458854&type=middle x 方向正应力解析解 σ x 为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458857&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458857&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458857&type=middle y 方向正应力解析解 σ y 为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458860&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458860&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458860&type=middle 切向应力解析解 σ xy 为 http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458864&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458864&type=small http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=4458864&type=middle 式中， G 为剪切模量， K Ⅰ 为应力强度因子. 在问题域上均布31×30的节点，支持域的无量纲因子 d max =3；采用罚函数法施加本质边界条件，惩罚因子 α =1×10 7 E .采用 J 积分 [ 15 15 ] 来计算应力强度因子 K Ⅰ ，选取的 J 积分域如 图 12 图 12 所示.采用可视性准则 [ 15 15 ] ，分别采用RKPM、QCRKPM和EFG方法3种方法进行计算. 计算以裂尖为圆心，0.2 l 为半径的半圆环上点的径向应力 σ rr 、环向应力 σ θθ 和切向应力 σ rθ ，数值结果及相对误差分别如 图 13 图 13 ~ 15 15 所示.由 图 13(a) 图 13(a) , 14(a) 14(a) , 15(a) 15(a) 可知，由QCRKPM计算得到的距离裂尖0.2 l 圆环上点的径向应力 σ rr 、环向应力 σ θθ 和切向应力 σ rθ 数值解的精度较高.由 图 13(b) 图 13(b) 可知，在计算径向应力 σ rr 时，3种形函数得到的相对误差都在8%以下，其中 θ 在0~0.5 rad范围内，EFG方法的精度最高，而QCRKPM的精度最差； θ 在0.5~1.5 rad范围内，QCRKPM的精度明显高于EFG方法； θ 在1.5~2.1 rad范围内，QCRKPM的计算精度明显高于RKPM.由 图 14(b) 图 14(b) 可知，在计算环向应力 σ θθ 时，QCRKPM和RKPM的相对误差均在15%以下，而 θ 在1.25~2.25 rad范围内EFG方法的相对误差最大，接近20%.由 图 15(b) 图 15(b) 可知，在计算切向应力 σ rθ 时， θ 在0~1.5 rad范围内，EFG方法的精度非常差，RKPM和QCRKPM的精度较高，而 θ 在1.8~2.5 rad范围内，QCRKPM的精度较差.EFG方法的计算时间最短，为0.322 s；RKPM和QCRKPM的计算时间相差不大，分别为1.259和1.225 s, 且QCRKPM的计算时间较RKPM略短. 表 9 表 9 给出了分别采用RKPM、QCRKPM及EFG方法时的应力强度因子 K Ⅰ 的计算值及相对误差.由 表 9 表 9 可知，这3种方法的计算精度均很高，均小于5%，其中RKPM的计算精度高于QCRKPM，QCRKPM的计算精度高于EFG方法. 2 K Ⅰ 解析解 2 RKPM 2 2 QCRKPM 2 2 EFG K Ⅰ 相对误差/% K Ⅰ 相对误差/% K I 相对误差/% 0.014 169 0.014 022 1.037 5 0.013 751 2.950 1 0.013 579 4.164 0 结论 1 本文将准凸重构核粒子法应用于弹性力学问题和弹性断裂问题，推导了相关计算公式，编写了相应MATLAB程序，分别对中心圆孔方板和边缘裂纹方板这两个算例进行了计算分析.数值结果表明：1)本文建立的弹性力学和弹性断裂问题的QCRKPM可以获得高精度的位移解、应力解和应力强度因子解，且不同的计算参数取值得到的弹性力学问题的系统能量范数误差略小于RKPM方法或与之相近，并显著小于EFG方法，说明本文方法的计算精度是可靠的；2)本文方法计算得到的裂尖应力解的相对误差平均在10%以下，总体计算精度高于EFG方法，且应力强度因子解的计算精度显著高于EFG方法，说明本文方法在弹性断裂问题分析中也具有较高准确性；3)本文方法的形函数计算时间与RKPM基本相当，说明准凸重构核近似不会显著增加形函数计算量和计算时间，在计算效率方面与传统重构核近似基本相同. 参考文献 LI D M, LIEW K M, CHENG Y M. Analyzing elastoplastic large deformation problems with the complex variable element-free Galerkin method [J]. Computational Mechanics , 2014, 53 (6):1149-1162. DOI: 10.1007/s00466-013-0954-4. LI D M, LIEW K M, CHENG Y M. An improved complex variable element-free Galerkin method for two-dimensional large deformation elastoplasticity problems [J]. Computer Methods in Applied Mecha-nics and Engineering , 2014, 269 (2):72-86. DOI: 10.1016/j.cma.2013.10.018. LI D M, BAI F N, CHENG Y M, et al . A novel complex variable element-free Galerkin method for two-dimensional large deformation problems [J]. Computer Methods in Applied Mechanics and Engineering , 2012, 233-236 (3):1-10. DOI: 10.1016/j.cma.2012.03.015. 李冬明, 彭妙娟, 程玉民.弹性大变形问题的复变量无单元Galerkin方法[J].中国科学:物理学 力学 天文学, 2011, 41 (8):1003-1014. LI D M, PENG M J, CHENG Y M. The complex variable element-free Galerkin (CVEFG) mothod for elastic large deformation problems [J]. Science China Physics, Mechanics & Astronomy , 2011, 41 (8):1003-1014 (Ch). LI D M, ZHANG Z, LIEW K M. A numerical framework for two-dimensional large deformation of inhomogeneous swelling of gels using the improved complex variable element-free Galerkin method [J]. Computer Methods in Applied Mechanics and Engineering , 2014, 274 (6):84-102. DOI: 10.1016/j.cma.2014.02.008. LI D M, ZHANG L W, LIEW K M. A three-dimensional element-free framework for coupled mechanical-diffusion induced nonlinear deformation of polymeric gels using the IMLS-Ritz method [J]. Computer Methods in Applied Mechanics and Engineering , 2015, 296: 232-247. DOI: 10.1016/j.cma.2015.05.011. LIU W K, JUN S, ZHANG Y F. Reproducing kernel particle methods [J]. International Journal for Numerical Methods in Fluids , 1995, 20 (8-9):1081-1106. DOI: 10.1002/fld.1650200824. LIU W K, CHEN Y, JUN S, et al . Overview and applications of the reproducing Kernel Particle methods [J]. Archives of Computational Methods in Engineering , 1996, 3 (1):3-80. DOI: 10.1007/BF02736130. MONAGHAN J J. An introduction to SPH [J]. Computer Physics Communications , 1988, 48 (1):89-96. DOI: 10.1016/0010-4655(88)90026-4. WANG D D, CHEN P J. Quasi-convex reproducing kernel meshfree method [J]. Computational Mecha-nics , 2014, 54 (3):689-709. DOI: 10.1007/s00466-014-1022-4. 陈波, 田利瑞, 李冬明.正交各向异性弹性力学的准凸重构核粒子法[J].华中科技大学学报(自然科学版), 2018, 46 (1):110-114. DOI:10.13245/j.hust.180121. CHEN B, TIAN L R, LI D M. Quasi-convex reproducing kernel particle method for orthotropic elasticity [J]. Journal of Huazhong University of Science and Technology (Natural Science Edition) , 2018, 46 (1):110-114. DOI:10.13245/j.hust.180121 (Ch). ZHANG H J, WANG D D. An isogeometric enriched quasi-convex meshfree formulation with application to material interface modeling [J]. Engineering Analysis with Boundary Elements . 2015, 60:37-50. DOI: 10.1016/j.enganabound.2015.03.016. ZHANG H J, WU J Z, WANG D D. Free vibration analysis of cracked thin plates by quasi-convex coupled isogeometric-meshfree method [J]. Frontiers of Structural and Civil Engineering , 2015, 9 (4):405-419. DOI: 10.1007/s11709-015-0310-1. LI D M, TIAN L R. Large deformation analysis of a gel using the complex variable element-free Galerkin method [J]. Applied Mathematical Modelling , 2018, 61:484-497. DOI: 10.1016/j.apm.2018.04.004. BELYTSCHKO T, LU Y Y, GU L. Element-free Galerkin methods [J]. International Journal for Numerical Methods in Engineering , 1994, 37 (2):229-256. DOI: 10.1002/nme.1620370205. 郦正能, 关志东.应用断裂力学[M].北京:北京航空航天大学出版社, 2012:10-12. LI Z N, GUAN Z D. Applied Fracture Mechanics [M]. Beijing: Beihang University Press, 2012:10-12 (Ch). 闫晓军, 胡殿印, 张小勇, 等.弹性力学[M].北京:清华大学出版社, 2015. YAN X J, HU D Y, ZHANG X Y, et al . Theory of elasticity [M]. Beijing: Tsinghua University Press, 2015 (Ch). 10.14188/j.1671-8836.2018.05.008 O343.1;O346.1;TU318 whdxxblxb-64-5-431 湖北省自然科学基金青年项目(2018CFB129);中央高校基本科研业务费专项资金(2016IVA022) 2018-04-15 2018-10-24 2020-12-09 版权所有© 2018 《武汉大学学报(理学版)》编辑部 Copyright © 2018 Journal of Wuhan University(Natural Science Edition). All rights reserved. 2018 武汉大学学报(理学版) 5 64