平面细长梁基于无网格径向基点插值的绝对节点坐标法
陈渊钊 章定国 黎亮
摘要： 为了消除或减弱传统绝对节点坐标法（Absolute Nodal Coordinate Formulation，ANCF）中缩减梁单元的“失真现象”，构造了一种适用于描述柔性梁绝对位形的无网格径向基点插值 （Radial Point Interpolation Method， RPIM）形函数，提出了柔性梁基于无网格RPIM的ANCF法。传统ANCF梁单元在描述纯弯曲悬臂梁的位形（一段圆弧）时，即便获得精确的单元节点坐标，通过梁单元插值得到的位形与悬臂梁的实际位形存在差异，即失真现象，悬臂梁越弯曲该差异越明显，失真越大。失真导致伪应变的产生，极大地影响数值求解的精度。而RPIM法采用一组场节点离散问题域，通过计算点支持域内的场节点构造形函数，计算点一般位于支持域的中心区域，不同计算点之间的支持域有较多重合的部分，加强了节点之间的联系，能更合理、准确地描述绝对位形，能有效减小失真。研究表明：基于RPIM的ANCF法较传统ANCF法精度更高、计算效率更快、对不等距分布节点的适应性更强，在大变形柔性多体系统动力学领域内具有推广性。
关键词： 多体系统动力学； 柔性梁； 绝对节点坐标法； 径向基点插值法； 失真
中图分类号： O313.7； O322 文献标志码： A 文章编号： 1004-4523（2018）02-0245-10
DOI：10.16385/j.cnki.issn.1004-4523.2018.02.007
引 言
ANCF法最早由Shabana[1]于1996年提出，该方法直接从连续介质力学出发，推导得到的动力学方程具有常数质量阵、不存在离心力项和科氏力项等特点，是一种适用于求解大转动、大变形动力学问题的方法[2-5]。一般地，ANCF法基于有限单元法，许多ANCF梁、板单元是通过多项式构造得到的[6-11]。然而Sanborn[12]发现三次多项式曲线的弯曲程度会影响参数点的分布，在曲线弯曲程度较大时会出现伪应变及曲线失真。张越[13]指出纯弯曲的ANCF缩减梁单元在弯曲程度较大时存在轴向拉伸变形，即存在伪拉伸应变，并认为其原因为弯曲、拉伸应变相互耦合，提出相应的解耦方法以消除伪应变，但其形函数仍与ANCF缩减梁单元一致。而Hyldahl [14]使用ANCF矩形壳单元离散一个四分之一空心圆盘时，发现一些单元之间存在“空白”，即相邻两个单元弯曲的共用边不重合，存在明显的缝隙。空心圆盘和ANCF壳单元可以分别看成是矩形板、ANCF矩形单元在变形后的形状，ANCF矩形单元基于精确的节点坐标插值得到的绝对位形与变形后的矩形板存在误差，即失真现象，说明基于多项式的ANCF单元并不能很好地描述变形体的绝对位形。事实上，若给定一段圆心角较大的圆弧曲线（纯弯曲悬臂梁的位形），将两端端点的位置及梯度精确赋值给ANCF缩减梁单元节点坐标，通过单元插值得到的曲线与实际曲线会存在偏差，曲线越弯曲偏差越大，因此，单个ANCF一维梁单元并不能很好地描述大曲率的曲线，单元的位形不同于纯弯曲细长梁的位形，即存在失真现象，从而使得单元中产生伪应变。
为避免伪应变的产生，应使ANCF精准地描述整个绝对位形，消除或减小失真，一般的解决方法便是大量增加单元数目，避免失真现象的产生，但这样会大大增加计算的规模及耗时。因此，构造一种精确描述绝对位形的离散方法是非常有必要的。计算力学中除了有限元法（Finite Element Method， FEM）以外，仍有不少新发展起来的离散方法，如无网格法（Meshfree Method）[15]。无网格法是一种不依赖单元网格的离散方法，该方法采用一组点来离散求解区域，并应用计算点支持域内的场节点构造近似函数，可以使用较多的场节点来提高近似函数的连续性[15]。FEM法两个相邻单元的计算点之间的联系仅仅依靠单元间的共同节点，而无网格法不同计算点之间的联系则是依靠对应支持域重合部分的场节点，比单元的共同节点多得多，使得无网格法节点之间的联系更加紧密，可以更好地描述变形场。而目前应用无网格法对绝对坐标场插值的研究尚未有报道。
现有的无网格法有多种，如再生核粒子法（Reproducing Kernel Particle Method， RKPM）[16-18]、无网格局部Perov-Galerkin法（Meshfree Local Perov-Galerkin Method， MLPGM）[19]、径向基点插值法[20]等。其中Liu GR等 [15，20] 提出的RPIM具有鲁棒性、高精度、Kronecker函数特性等特点，因此多被用于固体、流体力学问题的研究中[20-25]。杜超凡等成功将无网格法RPIM应用到柔性梁一次近似刚柔耦合模型[26]，并在柔性多体系统动力学问题上取得了良好的结果，说明无网格法在柔性多体系统动力学具有推广性[26-28]。
鉴于无网格RPIM法的优势，本文构造一种适用于描述柔性梁绝对位形的RPIM形函数，并基于连续介质力学和哈密顿原理建立平面细长梁系统的动力学方程，提出基于无网格RPIM的绝对节点坐标法。本文首先阐明ANCF缩减梁单元中失真现象的客观存在性，通过与本文方法对比，分析伪应变产生原因，同时检验本文方法描述曲线的性能。再通过大变形静力学、动力学算例，分析失真对传统ANCF缩减梁单元精度的影响，同时比较说明本文方法精确处理大变形问题的精度及效率。
1 柔性梁的ANCF模型
本文研究对象为平面梁系统，采用Euler-Bernouli假设，不考虑剪切变形，并认为变形时梁的横截面仍保持为平面且与中轴线垂直避免伪应变的产生，其实就是避免失真现象，但失真现象很难完全消除，仅能尽量减小失真的程度，所以需要增加離散的节点数。事实上，除了加密离散以外，还可以采用其他更合理的插值方法来减小失真。若取长度为3 m整圆，等长划分三段，即取四节点，应用传统ANCF法和本文方法进行插值，如图4所示，传统ANCF插值结果偏差较大，失真明显，本文方法所得结果基本与圆重合，说明本文方法较传统ANCF描述曲线更合理，失真的程度较传统ANCF法小得多。表3给出了此时基于两种方法求得的势能。在传统ANCF结果中，一个长度为1 m、曲率κ=2π/3单元的伪应变能为7.3985 J（如表2所示），三个单元的伪应变之和为22.1958 J，三个单元所得伪势能为一个单元伪势能的三倍，说明传统ANCF单元之间联系简单，伪势能较大，而本文方法所得结果较传统ANCF法小得多，说明本文方法能有效减小失真。
4 仿真算例
由上一节得知ANCF描述弯曲较大的曲线时会失真，从而产生伪应变，这种现象对ANCF求解大变形问题时有着怎样的影响，下面通过两个简单的例子进行讨论。
4.1 受集中力矩作用的悬臂梁
考虑图5所示的悬臂梁，梁末端受集中力矩τ的作用，研究其静力学问题。梁长L=1.8 m，横截面积A=2.5 cm2，截面惯性矩I=0.130 cm4，弹性模量E=68.95 GPa，密度ρ=2766.67 kg/m3，即梁参数选取与文献[30-32]一致。
图5 受集中力矩作用的悬臂梁
Fig.5 Cantilever beam under a moment
若梁末端受到力矩τ作用，力矩大小为τ=λEIL（51）其中λ为一常数。梁变形后的轴线为一段弧线，相应的曲率为κ=τLEI=λ（52） 当λ=2π时柔性梁的位形为一整圆，梁末端的水平、竖直方向上的绝对位移为0。选取不同数目的等距离分布的节点和不等距离分布的节点，对比本文方法和传统ANCF法的差异。
等距分布的节点在梁上第i个节点的水平绝对坐标为xi=（i-1）Lm， i=1，…，m（53） 不等距分布的节点选取切比雪夫插值节点，梁两端节点分布较密集，中间稀疏，梁上第i个节点的水平绝对坐标为xi=L2-L2cosi-1πm， i=1，…，m（54） 传统ANCF法在使用这些节点时为相邻两节点构成一个单元。
在计算时，需忽略动能的影响，故静力学方程为Qτ-F=0（55） 利用Newton-Raphson迭代即可求解静力学方程。
表4为梁末端的水平方向的绝对位移，其解析解为0。从表4可以看出，在节点数较少时，传统ANCF法计算结果不准确，本文方法计算结果仍有较高精度。此外，使用等距节点的传统ANCF法仿真数值结果的精度要比使用不等距节点好，甚至高一个数量级，说明传统ANCF法的计算精度受节点的分布情况影响很大，对网格依赖性高，而本文方法几乎不受节点分布影响，对节点分布的依赖性很弱，由此可说明本文方法较传统ANCF法对不等距分布的节点更具适应性，数值稳定性高。
表4 梁末端的水平方向绝对位移（单位：m）
Tab.4 Horizontal position of the free end of beam （Unit：m）节点数ANCFRPIM等距不等距等距不等距4-0.31663-0.36876-0.003530.003815-0.30480-0.514720.001210.000136-0.13772-0.428260.000080.000027-0.04940-0.29022-0.00006-0.000018-0.01829-0.16898-0.000070.000039-0.00767-0.08916-0.00003〖〗0.0000210-0.00368-0.04552-0.00001-0.00001 图6为5节点等距分布时，应用两种方法分别计算λ=0，0.5π，π，1.5π，2π时柔性梁的位形所得结果。如图所示，λ=0.5π时，两种方法基本与解析解吻合，说明求解柔性梁变形较小的静力学问题时，两种方法使用较少节点也可得到较为可信的结果，同时也验证了本文方法的正确性。当λ=2π时，梁末端点A应与悬臂端点O重合，传统ANCF法计算得到的点A位置与解析解相差较大距离，且梁的位形与解析解误差明显，本文方法与解析解的符合度较好，说明了本文方法的正确性，且精度较ANCF法高。事实上，λ=0.5π时，每个ANCF缩减梁单元的转角θi很小，最大的仅有22.5°，所以此时的失真现象非常小，可以忽略；而λ=2π时，单元转角θi最大的有74°，失真现象明显，伪应变较大，会产生较大误差，极大地影响了传统ANCF法的求解精度。
图6 受集中力矩作用时柔性梁的位形
Fig.6 Position of beam under a moment图7为λ=2π时柔性梁的位形。此时，传统ANCF法取10个节点形成单元计算，本文方法取5个节点计算。如图7所示，传统ANCF法取等距分布的节点时计算结果基本与解析解重合，此时最大的单元转角θi有39°，失真现象基本可以忽略，而使用不等距节点时误差较大，此时最大的单元转角θi有60°，失真现象明显；本文方法取等距节点和不等距节点均与解析解基本重合，与表4所得结论一致，说明传统ANCF法的网格对精度影响很大，稳定性较低，本文方法则没有这样的缺陷。
图7 λ=2π时柔性梁的位形
Fig.7 Position of beam under the moment （λ=2π）4.2 重力單摆
考虑如图8所示的重力单摆，梁的材料参数仍与文献[30-32]一致。
图8 重力单摆
Fig.8 Gravity pendulum
图9（a）和（b）分别为梁下落过程中末端的竖向、水平绝对位移，此时传统ANCF法划分四单元、本文方法取五节点。此时两种方法仿真结果基本重合，ANCF缩减梁单元转角θi最大值为0.15°，此时的梁在小变形状态，说明在求解柔性梁小变形动力学问题时，两种方法使用较少的单元（节点）数不会导致明显的失真，可得到较好的数值结果。
图9 自由单摆末端的绝对位移
Fig.9 Position of the free end of the pendulum
图10 E=68.95 MPa时自由单摆末端的竖直绝对位移
Fig.10 Vertical position of free end of the pendulum
（E=68.95 MPa）图10给出了E=68.95 MPa时柔性梁末端的竖向绝对位移。从图10看，传统ANCF 和RPIM分别在13个单元、10个节点时基本收敛，事实上，以某一点的位移作为收敛判断的依据是不够准确的，如图11所示，传统ANCF法使用30个单元和使用13、17个单元计算得到的2.5 s时柔性梁的位形有轻微的误差，而与RPIM使用11个节点的计算结果更符合，说明RPIM精度更好。值得注意的是，传统ANCF法使用13个单元时，弯曲最大的单元θi为39.9°；使用17个单元时，弯曲最大的单元θi为32.9°。说明少量增加单元个数并不能迅速减小某些单元的弯曲程度，而大量增加单元个数则大大增加了计算的耗时及规模，若采用RPIM则可有效提高计算精度。
圖11 2.5 s时柔性梁的位形
Fig.11 Deformed shape of the pendulum （t=2.5 s）
图12 梁上每点的曲率
Fig.12 Curvature of points of the beam
由于应力应变比位移更为重要，故以2.5 s时柔性梁上每点的曲率来衡量两种方法的收敛性，如图12所示，传统ANCF法和RPIM分别在使用23个单元和14个节点时收敛。表5列出了两种方法应用广义α法[33]计算时的相对耗时，传统ANCF法使用23个单元计算耗时比RPIM使用14个节点计算耗时要多得多，可表明RPIM精确处理大变形问题的高效性。图13给出了柔性单摆下落过程中的能量变化曲线，系统的动能、势能、应变能之和恒为零，即系统能量守恒，再次说明RPIM法的正确性。
表5 两种方法的相对耗时
Tab.5 The relative time consuming of two methodsANCFRPIM19
单元23
单元29
单元12
节点14
节点15
节点1.4272.25784.26311.3421.494注：该算例采用广义α法计算[33]
图13 能量平衡曲线
Fig.13 Energy balance curve5 结 论
本文应用无网格RPIM法对平面柔性梁的绝对位形进行插值，利用支持域内的场节点所对应的径向基函数及其导数来构造绝对位形的插值形函数，并通过哈密顿原理建立柔性梁系统的动力学方程，提出基于无网格RPIM的ANCF法。研究表明：
1） 传统ANCF缩减梁单元在求解大变形问题时会存在失真现象，单元越弯曲失真现象就越明显，产生的伪应变也越大，严重降低了数值求解的精度。故应控制每个单元的弯曲程度，加密网格、大量增加单元个数，但也同时大大增加了计算成本。
2） 本文方法能合理、准确地描述柔性梁变形后的绝对位形曲线，有效减小失真的程度，降低伪应变的影响。通过静力学和动力学仿真算例可以说明：本文方法较传统ANCF法精度更高、计算效率更快，对不等距分布节点的适应性更强。
参考文献：
[1] Shabana A A. An absolute nodal coordinate formulation for the large rotation and deformation analysis of flexible bodies [R].Department of Mechanical and Industrial Engineering， University of Illinois at Chicago， 1996.
[2] Yakoub R Y， Shabana A A. Three dimensional absolute nodal coordinate formulation for beam elements： implementation and applications [J]. Journal of Mechanical Design， 2001， 123（4）： 614—621.
[3] Shabana A A， Yakoub R Y. Three dimensional absolute nodal coordinate formulation for beam elements： theory [J]. Journal of Mechanical Design， 2001， 123（4）： 606—613.
[4] Omar M A， Shabana A A. A two-dimensional shear deformable beam for large rotation and deformation problems [J]. Journal of Sound and Vibration， 2001， 243（3）： 565—576.
[5] Berzeri M， Shabana A A. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation [J]. Journal of Sound and Vibration， 2000， 235（4）： 539—565.
[6] Olshevskiy A， Dmitrochenko O， Dai M D， et al. The simplest 3-， 6- and 8-noded fully-parameterized ANCF plate elements using only transverse slopes [J]. Multibody System Dynamics， 2014， 34（1）： 23—51.
[7] Yan D， Liu C， Tian Q， et al. A new curved gradient deficient shell element of absolute nodal coordinate formulation for modeling thin shell structures [J]. Nonlinear Dynamics， 2013， 74（1-2）： 153—164.
[8] Dmitrochenko O， Mikkola A. Two simple triangular plate elements based on the absolute nodal coordinate formulation[J]. Journal of Computational and Nonlinear Dynamics， 2008， 3（4）： 041012.
[9] Gerstmayr J， Shabana A A. Analysis of thin beams and cables using the absolute nodal coordinate formulation [J]. Nonlinear Dynamics， 2006， 45（1-2）： 109—130.
[10]Dufva K E， Sopanen J T， Mikkola A M. A two-dimensional shear deformable beam element based on the absolute nodal coordinate formulation[J]. Journal of Sound and Vibration， 2005， 280（3-5）： 719—738.
[11]Dufva K， Shabana A A. Analysis of thin plate structures using the absolute nodal coordinate formulation[J]. Proceedings of the Institution of Mechanical Engineers， Part K： Journal of Multi-body Dynamics， 2005， 219（4）： 345—355.
[12]Sanborn G G， Choi J， Choi J H. Curve-induced distortion of polynomial space curves， flat-mapped extension modeling， and their impact on ANCF thin-plate finite elements [J]. Multibody System Dynamics， 2011， 26（2）： 191—211.
[13]張 越， 赵 阳， 谭春林，等. ANCF索梁单元应变耦合问题与模型解耦[J]. 力学学报， 2016， 48（6）： 1406—1415.
Zhang Yue， Zhao Yan，Tan Chunlin， et al. The strain coupling problem and model decoupling of ANCF cable/beam element[J]. Chinese Journal of Theoretical and Applied Mechanics，2016，48（6）：1406—1415.
[14]Hyldahl P， Mikkola A M， Balling O， et al. Behavior of thin rectangular ANCF shell elements in various mesh configurations [J]. Nonlinear Dynamics， 2014， 78（2）： 1277—1291.
[15]刘桂荣，顾元通. 无网格法理论及程序设计[M].王建明，周学军，译. 济南：山东大学出版社，2007：45—133.
Liu G R， Gu Y T. An Introduction to Meshfree Methods and Their Programming[M]. Jinan：Shandong University Press， 2007：45—133.
[16]Chen J S， Pan C， Roque C M O L， et al. A lagrangian reproducing kernel particle method for metal forming analysis[J]. Computational Mechanics， 1998， 22（3）： 289—307.
[17]Chen J S， Pan C， Wu C T， et al. Reproducing kernel particle methods for large deformation analysis of non-linear structures[J]. Computer Methods in Applied Mechanics & Engineering， 1996， 139（1-4）： 195—227.
[18]Liu W K， Jun S， Zhang Y F. Reproducing kernel particle methods. Int J Numer Meth Eng[J]. International Journal for Numerical Methods in Fluids， 1995， 20（8-9）： 1081—1106.
[19]Xiong Y B， Long S Y. Local petrov-galerkin method for a thin plate[J]. Applied Mathematics and Mechanics， 2004， 25（2）： 210—218.
[20]Wang J G， Liu G R. A point interpolation meshless method based on radial basis functions[J]. International Journal for Numerical Methods in Engineering， 2002， 54（11）： 1623—1648.
[21]Cui X Y， Feng H， Li G Y， et al. A cell-based smoothed radial point interpolation method （CS-RPIM） for three-dimensional solids[J]. Engineering Analysis with Boundary Elements， 2015， 50： 474—485.
[22]Pilafkan R， Folkow P D， Darvizeh M， et al. Three dimensional frequency analysis of bidirectional functionally graded thick cylindrical shells using a radial point interpolation method （RPIM） [J]. European Journal of Mechanics / A Solids， 2013， 39（5）： 26—34.
[23]Gu Y T， Wang W L， Fu Q. An enriched radial point interpolation method （e-RPIM） for the analysis of crack tip[J]. Engineering Fracture Mechanics， 2010， 78（1）： 175—190.
[24]Chen S L， Li Y X. An efficient RPIM for simulating wave motions in saturated porous media[J]. International Journal of Solids & Structures， 2008， 45（25）： 6316—6332.
[25]Dai K Y， Liu G R， Han X， et al. Inelastic analysis of 2D solids using a weak-form RPIM based on deformation theory[J]. Computer Methods in Applied Mechanics & Engineering， 2006， 195（33）： 4179—4193.
[26]杜超凡， 章定国， 洪嘉振. 径向基点插值法在旋转柔性梁动力学中的应用[J]. 力学学报， 2015， 47（2）： 279—288.
Du Chaofan， Zhang Dingguo， Hong Jiazhen. A meshfree method based on radial point interpolation method for the dynamic analysis of rotating flexible beams[J]. Chinese Journal of Theoretical and Applied Mechanics， 2015，47（2）： 279—288.
[27]杜超凡， 章定国. 光滑节点插值法：计算固有频率下界值的新方法[J]. 力学学报， 2015， 47（5）： 839—847.
Du Chaofan， Zhang Dingguo. Node-based smoothed point interpolation method： a new method for computing lower bound of natural frequency[J]. Chinese Journal of Theoretical and Applied Mechanics， 2015， 47（5）： 839—847.
[28]杜超凡， 章定国. 基于无网格点插值法的旋转悬臂梁的动力学分析[J]. 物理学报， 2015，64（3）： 396—405.
Du Chaofan， Zhang Dingguo. A meshfree method based on point interpolation for dynamic analysis of rotating cantilever beams[J]. Acta Phys. Sin.， 2015， 64（3）： 396—405.
[29]章孝順，章定国，陈思佳，等. 基于绝对节点坐标法的大变形柔性梁几种动力学模型研究[J].物理学报， 2016， 65（9）：148—157.
Zhang Xiaoshun， Zhang Dingguo， Chen Sijia，et al. Several dynamic models of a large deformation flexible beam based on the absolute nodal coordinate formulation[J]. Acta Phys. Sin.， 2016， 65（9）：148—157.
[30]陈思佳， 章定国， 洪嘉振. 大变形旋转柔性梁的一种高次刚柔耦合动力学模型[J]. 力学学报， 2013， 45（2）：251—256.
Chen Sijia， Zhang Dingguo， Hong Jiazhen. A high-order rigid-flexible coupling model of a rotating flexible beam under large deformation[J]. Chinese Journal of Theoretical and Applied Mechanics， 2013， 45（2）：251—256.
[31]李 彬， 刘锦阳. 大变形柔性梁系统的绝对坐标方法[J]. 上海交通大学学报， 2005， 39（5）： 827—831.
Li Bin， Liu Jinyang. Application of absolute nodal coordination formulation in flexible beams with large deformation[J]. Journal of Shanghai Jiaotong University， 2005， 39（5）： 827—831.
[32]章孝順， 章定国， 洪嘉振. 考虑曲率纵向变形效应的大变形柔性梁刚柔耦合动力学建模与仿真[J]. 力学学报，2016，48（3）：692—701.
Zhang X H， Zhang D G， Hong J Z. Rigid-flexible couple dynamic modeling and simulation with the longitudinal deformation induced curvature effect for a rotating flexible beam under larger deformation [J]. Chinese Journal of Theoretical and Applied Mechanics，2016，48（3）：692—701.
[33]Amold M， Brüls O. Convergence of the generalized-α scheme for constrained mechanical systems[J]. Multibody System Dynamics， 2007，18（2）：185—202.
An absolute nodal coordinate formulation based on
radial point interpolation method for planar slender beams
CHEN Yuan-zhao， ZHANG Ding-guo， LI Liang
（School of Science， Nanjing University of Science and Technology， Nanjing 210094， China）
Abstract： In order to alleviate or eliminate the ‘distortion phenomenon of the deficient beam elements in the traditional absolute nodal coordinate formulation （ANCF）， an ANCF based on radial point interpolation method （RPIM） for flexible beams is proposed in which a new RPIM shape functions are constructed to describe the absolute configuration of flexible beams. For a pure bending cantilever beam （an arc）， there is always difference between configuration of the beams by using the gradient deficient beam elements and the actual configuration of the beam in the traditional ANCF. The difference， namely ‘distortion phenomenon， becomes more obvious with the increase of the bending deformation of the beam， which may cause the pseudo strain and have serious influence on accuracy of numerical solution. In the present method， the RPIM is used to discretize the deformation field through a set of field nodes and the shape functions are generally established based on field nodes within a support domain of the calculating point. The calculating points are generally located in the central region of the support domain， and the support domain of different calculating points can have more coincident parts. Thus， the connection between field nodes is strengthened， which makes the method describe the configuration in more reasonable and effectively alleviate influence of the distortion and pseudo strain. The simulation results show that the proposed method has higher calculation accuracy and efficiency and is more adaptive for the non-equidistant nodes compared with the traditional ANCF， which can be further extended in the dynamic field of flexible multi-body system.
Key words： multibody dynamics； flexible beam； absolute nodal coordinate formulation； radial point interpolation method； distortion
摘要： 为了消除或减弱传统绝对节点坐标法（Absolute Nodal Coordinate Formulation，ANCF）中缩减梁单元的“失真现象”，构造了一种适用于描述柔性梁绝对位形的无网格径向基点插值 （Radial Point Interpolation Method， RPIM）形函数，提出了柔性梁基于无网格RPIM的ANCF法。传统ANCF梁单元在描述纯弯曲悬臂梁的位形（一段圆弧）时，即便获得精确的单元节点坐标，通过梁单元插值得到的位形与悬臂梁的实际位形存在差异，即失真现象，悬臂梁越弯曲该差异越明显，失真越大。失真导致伪应变的产生，极大地影响数值求解的精度。而RPIM法采用一组场节点离散问题域，通过计算点支持域内的场节点构造形函数，计算点一般位于支持域的中心区域，不同计算点之间的支持域有较多重合的部分，加强了节点之间的联系，能更合理、准确地描述绝对位形，能有效减小失真。研究表明：基于RPIM的ANCF法较传统ANCF法精度更高、计算效率更快、对不等距分布节点的适应性更强，在大变形柔性多体系统动力学领域内具有推广性。
关键词： 多体系统动力学； 柔性梁； 绝对节点坐标法； 径向基点插值法； 失真
中图分类号： O313.7； O322 文献标志码： A 文章编号： 1004-4523（2018）02-0245-10
DOI：10.16385/j.cnki.issn.1004-4523.2018.02.007
引 言
ANCF法最早由Shabana[1]于1996年提出，该方法直接从连续介质力学出发，推导得到的动力学方程具有常数质量阵、不存在离心力项和科氏力项等特点，是一种适用于求解大转动、大变形动力学问题的方法[2-5]。一般地，ANCF法基于有限单元法，许多ANCF梁、板单元是通过多项式构造得到的[6-11]。然而Sanborn[12]发现三次多项式曲线的弯曲程度会影响参数点的分布，在曲线弯曲程度较大时会出现伪应变及曲线失真。张越[13]指出纯弯曲的ANCF缩减梁单元在弯曲程度较大时存在轴向拉伸变形，即存在伪拉伸应变，并认为其原因为弯曲、拉伸应变相互耦合，提出相应的解耦方法以消除伪应变，但其形函数仍与ANCF缩减梁单元一致。而Hyldahl [14]使用ANCF矩形壳单元离散一个四分之一空心圆盘时，发现一些单元之间存在“空白”，即相邻两个单元弯曲的共用边不重合，存在明显的缝隙。空心圆盘和ANCF壳单元可以分别看成是矩形板、ANCF矩形单元在变形后的形状，ANCF矩形单元基于精确的节点坐标插值得到的绝对位形与变形后的矩形板存在误差，即失真现象，说明基于多项式的ANCF单元并不能很好地描述变形体的绝对位形。事实上，若给定一段圆心角较大的圆弧曲线（纯弯曲悬臂梁的位形），将两端端点的位置及梯度精确赋值给ANCF缩减梁单元节点坐标，通过单元插值得到的曲线与实际曲线会存在偏差，曲线越弯曲偏差越大，因此，单个ANCF一维梁单元并不能很好地描述大曲率的曲线，单元的位形不同于纯弯曲细长梁的位形，即存在失真现象，从而使得单元中产生伪应变。
为避免伪应变的产生，应使ANCF精准地描述整个绝对位形，消除或减小失真，一般的解决方法便是大量增加单元数目，避免失真现象的产生，但这样会大大增加计算的规模及耗时。因此，构造一种精确描述绝对位形的离散方法是非常有必要的。计算力学中除了有限元法（Finite Element Method， FEM）以外，仍有不少新发展起来的离散方法，如无网格法（Meshfree Method）[15]。无网格法是一种不依赖单元网格的离散方法，该方法采用一组点来离散求解区域，并应用计算点支持域内的场节点构造近似函数，可以使用较多的场节点来提高近似函数的连续性[15]。FEM法两个相邻单元的计算点之间的联系仅仅依靠单元间的共同节点，而无网格法不同计算点之间的联系则是依靠对应支持域重合部分的场节点，比单元的共同节点多得多，使得无网格法节点之间的联系更加紧密，可以更好地描述变形场。而目前应用无网格法对绝对坐标场插值的研究尚未有报道。
现有的无网格法有多种，如再生核粒子法（Reproducing Kernel Particle Method， RKPM）[16-18]、无网格局部Perov-Galerkin法（Meshfree Local Perov-Galerkin Method， MLPGM）[19]、径向基点插值法[20]等。其中Liu GR等 [15，20] 提出的RPIM具有鲁棒性、高精度、Kronecker函数特性等特点，因此多被用于固体、流体力学问题的研究中[20-25]。杜超凡等成功将无网格法RPIM应用到柔性梁一次近似刚柔耦合模型[26]，并在柔性多体系统动力学问题上取得了良好的结果，说明无网格法在柔性多体系统动力学具有推广性[26-28]。
鉴于无网格RPIM法的优势，本文构造一种适用于描述柔性梁绝对位形的RPIM形函数，并基于连续介质力学和哈密顿原理建立平面细长梁系统的动力学方程，提出基于无网格RPIM的绝对节点坐标法。本文首先阐明ANCF缩减梁单元中失真现象的客观存在性，通过与本文方法对比，分析伪应变产生原因，同时检验本文方法描述曲线的性能。再通过大变形静力学、动力学算例，分析失真对传统ANCF缩减梁单元精度的影响，同时比较说明本文方法精确处理大变形问题的精度及效率。
1 柔性梁的ANCF模型
本文研究对象为平面梁系统，采用Euler-Bernouli假设，不考虑剪切变形，并认为变形时梁的横截面仍保持为平面且与中轴线垂直避免伪应变的产生，其实就是避免失真现象，但失真现象很难完全消除，仅能尽量减小失真的程度，所以需要增加離散的节点数。事实上，除了加密离散以外，还可以采用其他更合理的插值方法来减小失真。若取长度为3 m整圆，等长划分三段，即取四节点，应用传统ANCF法和本文方法进行插值，如图4所示，传统ANCF插值结果偏差较大，失真明显，本文方法所得结果基本与圆重合，说明本文方法较传统ANCF描述曲线更合理，失真的程度较传统ANCF法小得多。表3给出了此时基于两种方法求得的势能。在传统ANCF结果中，一个长度为1 m、曲率κ=2π/3单元的伪应变能为7.3985 J（如表2所示），三个单元的伪应变之和为22.1958 J，三个单元所得伪势能为一个单元伪势能的三倍，说明传统ANCF单元之间联系简单，伪势能较大，而本文方法所得结果较传统ANCF法小得多，说明本文方法能有效减小失真。
4 仿真算例
由上一节得知ANCF描述弯曲较大的曲线时会失真，从而产生伪应变，这种现象对ANCF求解大变形问题时有着怎样的影响，下面通过两个简单的例子进行讨论。
4.1 受集中力矩作用的悬臂梁
考虑图5所示的悬臂梁，梁末端受集中力矩τ的作用，研究其静力学问题。梁长L=1.8 m，横截面积A=2.5 cm2，截面惯性矩I=0.130 cm4，弹性模量E=68.95 GPa，密度ρ=2766.67 kg/m3，即梁参数选取与文献[30-32]一致。
图5 受集中力矩作用的悬臂梁
Fig.5 Cantilever beam under a moment
若梁末端受到力矩τ作用，力矩大小为τ=λEIL（51）其中λ为一常数。梁变形后的轴线为一段弧线，相应的曲率为κ=τLEI=λ（52） 当λ=2π时柔性梁的位形为一整圆，梁末端的水平、竖直方向上的绝对位移为0。选取不同数目的等距离分布的节点和不等距离分布的节点，对比本文方法和传统ANCF法的差异。
等距分布的节点在梁上第i个节点的水平绝对坐标为xi=（i-1）Lm， i=1，…，m（53） 不等距分布的节点选取切比雪夫插值节点，梁两端节点分布较密集，中间稀疏，梁上第i个节点的水平绝对坐标为xi=L2-L2cosi-1πm， i=1，…，m（54） 传统ANCF法在使用这些节点时为相邻两节点构成一个单元。
在计算时，需忽略动能的影响，故静力学方程为Qτ-F=0（55） 利用Newton-Raphson迭代即可求解静力学方程。
表4为梁末端的水平方向的绝对位移，其解析解为0。从表4可以看出，在节点数较少时，传统ANCF法计算结果不准确，本文方法计算结果仍有较高精度。此外，使用等距节点的传统ANCF法仿真数值结果的精度要比使用不等距节点好，甚至高一个数量级，说明传统ANCF法的计算精度受节点的分布情况影响很大，对网格依赖性高，而本文方法几乎不受节点分布影响，对节点分布的依赖性很弱，由此可说明本文方法较传统ANCF法对不等距分布的节点更具适应性，数值稳定性高。
表4 梁末端的水平方向绝对位移（单位：m）
Tab.4 Horizontal position of the free end of beam （Unit：m）节点数ANCFRPIM等距不等距等距不等距4-0.31663-0.36876-0.003530.003815-0.30480-0.514720.001210.000136-0.13772-0.428260.000080.000027-0.04940-0.29022-0.00006-0.000018-0.01829-0.16898-0.000070.000039-0.00767-0.08916-0.00003〖〗0.0000210-0.00368-0.04552-0.00001-0.00001 图6为5节点等距分布时，应用两种方法分别计算λ=0，0.5π，π，1.5π，2π时柔性梁的位形所得结果。如图所示，λ=0.5π时，两种方法基本与解析解吻合，说明求解柔性梁变形较小的静力学问题时，两种方法使用较少节点也可得到较为可信的结果，同时也验证了本文方法的正确性。当λ=2π时，梁末端点A应与悬臂端点O重合，传统ANCF法计算得到的点A位置与解析解相差较大距离，且梁的位形与解析解误差明显，本文方法与解析解的符合度较好，说明了本文方法的正确性，且精度较ANCF法高。事实上，λ=0.5π时，每个ANCF缩减梁单元的转角θi很小，最大的仅有22.5°，所以此时的失真现象非常小，可以忽略；而λ=2π时，单元转角θi最大的有74°，失真现象明显，伪应变较大，会产生较大误差，极大地影响了传统ANCF法的求解精度。
图6 受集中力矩作用时柔性梁的位形
Fig.6 Position of beam under a moment图7为λ=2π时柔性梁的位形。此时，传统ANCF法取10个节点形成单元计算，本文方法取5个节点计算。如图7所示，传统ANCF法取等距分布的节点时计算结果基本与解析解重合，此时最大的单元转角θi有39°，失真现象基本可以忽略，而使用不等距节点时误差较大，此时最大的单元转角θi有60°，失真现象明显；本文方法取等距节点和不等距节点均与解析解基本重合，与表4所得结论一致，说明传统ANCF法的网格对精度影响很大，稳定性较低，本文方法则没有这样的缺陷。
图7 λ=2π时柔性梁的位形
Fig.7 Position of beam under the moment （λ=2π）4.2 重力單摆
考虑如图8所示的重力单摆，梁的材料参数仍与文献[30-32]一致。
图8 重力单摆
Fig.8 Gravity pendulum
图9（a）和（b）分别为梁下落过程中末端的竖向、水平绝对位移，此时传统ANCF法划分四单元、本文方法取五节点。此时两种方法仿真结果基本重合，ANCF缩减梁单元转角θi最大值为0.15°，此时的梁在小变形状态，说明在求解柔性梁小变形动力学问题时，两种方法使用较少的单元（节点）数不会导致明显的失真，可得到较好的数值结果。
图9 自由单摆末端的绝对位移
Fig.9 Position of the free end of the pendulum
图10 E=68.95 MPa时自由单摆末端的竖直绝对位移
Fig.10 Vertical position of free end of the pendulum
（E=68.95 MPa）图10给出了E=68.95 MPa时柔性梁末端的竖向绝对位移。从图10看，传统ANCF 和RPIM分别在13个单元、10个节点时基本收敛，事实上，以某一点的位移作为收敛判断的依据是不够准确的，如图11所示，传统ANCF法使用30个单元和使用13、17个单元计算得到的2.5 s时柔性梁的位形有轻微的误差，而与RPIM使用11个节点的计算结果更符合，说明RPIM精度更好。值得注意的是，传统ANCF法使用13个单元时，弯曲最大的单元θi为39.9°；使用17个单元时，弯曲最大的单元θi为32.9°。说明少量增加单元个数并不能迅速减小某些单元的弯曲程度，而大量增加单元个数则大大增加了计算的耗时及规模，若采用RPIM则可有效提高计算精度。
圖11 2.5 s时柔性梁的位形
Fig.11 Deformed shape of the pendulum （t=2.5 s）
图12 梁上每点的曲率
Fig.12 Curvature of points of the beam
由于应力应变比位移更为重要，故以2.5 s时柔性梁上每点的曲率来衡量两种方法的收敛性，如图12所示，传统ANCF法和RPIM分别在使用23个单元和14个节点时收敛。表5列出了两种方法应用广义α法[33]计算时的相对耗时，传统ANCF法使用23个单元计算耗时比RPIM使用14个节点计算耗时要多得多，可表明RPIM精确处理大变形问题的高效性。图13给出了柔性单摆下落过程中的能量变化曲线，系统的动能、势能、应变能之和恒为零，即系统能量守恒，再次说明RPIM法的正确性。
表5 两种方法的相对耗时
Tab.5 The relative time consuming of two methodsANCFRPIM19
单元23
单元29
单元12
节点14
节点15
节点1.4272.25784.26311.3421.494注：该算例采用广义α法计算[33]
图13 能量平衡曲线
Fig.13 Energy balance curve5 结 论
本文应用无网格RPIM法对平面柔性梁的绝对位形进行插值，利用支持域内的场节点所对应的径向基函数及其导数来构造绝对位形的插值形函数，并通过哈密顿原理建立柔性梁系统的动力学方程，提出基于无网格RPIM的ANCF法。研究表明：
1） 传统ANCF缩减梁单元在求解大变形问题时会存在失真现象，单元越弯曲失真现象就越明显，产生的伪应变也越大，严重降低了数值求解的精度。故应控制每个单元的弯曲程度，加密网格、大量增加单元个数，但也同时大大增加了计算成本。
2） 本文方法能合理、准确地描述柔性梁变形后的绝对位形曲线，有效减小失真的程度，降低伪应变的影响。通过静力学和动力学仿真算例可以说明：本文方法较传统ANCF法精度更高、计算效率更快，对不等距分布节点的适应性更强。
参考文献：
[1] Shabana A A. An absolute nodal coordinate formulation for the large rotation and deformation analysis of flexible bodies [R].Department of Mechanical and Industrial Engineering， University of Illinois at Chicago， 1996.
[2] Yakoub R Y， Shabana A A. Three dimensional absolute nodal coordinate formulation for beam elements： implementation and applications [J]. Journal of Mechanical Design， 2001， 123（4）： 614—621.
[3] Shabana A A， Yakoub R Y. Three dimensional absolute nodal coordinate formulation for beam elements： theory [J]. Journal of Mechanical Design， 2001， 123（4）： 606—613.
[4] Omar M A， Shabana A A. A two-dimensional shear deformable beam for large rotation and deformation problems [J]. Journal of Sound and Vibration， 2001， 243（3）： 565—576.
[5] Berzeri M， Shabana A A. Development of simple models for the elastic forces in the absolute nodal co-ordinate formulation [J]. Journal of Sound and Vibration， 2000， 235（4）： 539—565.
[6] Olshevskiy A， Dmitrochenko O， Dai M D， et al. The simplest 3-， 6- and 8-noded fully-parameterized ANCF plate elements using only transverse slopes [J]. Multibody System Dynamics， 2014， 34（1）： 23—51.
[7] Yan D， Liu C， Tian Q， et al. A new curved gradient deficient shell element of absolute nodal coordinate formulation for modeling thin shell structures [J]. Nonlinear Dynamics， 2013， 74（1-2）： 153—164.
[8] Dmitrochenko O， Mikkola A. Two simple triangular plate elements based on the absolute nodal coordinate formulation[J]. Journal of Computational and Nonlinear Dynamics， 2008， 3（4）： 041012.
[9] Gerstmayr J， Shabana A A. Analysis of thin beams and cables using the absolute nodal coordinate formulation [J]. Nonlinear Dynamics， 2006， 45（1-2）： 109—130.
[10]Dufva K E， Sopanen J T， Mikkola A M. A two-dimensional shear deformable beam element based on the absolute nodal coordinate formulation[J]. Journal of Sound and Vibration， 2005， 280（3-5）： 719—738.
[11]Dufva K， Shabana A A. Analysis of thin plate structures using the absolute nodal coordinate formulation[J]. Proceedings of the Institution of Mechanical Engineers， Part K： Journal of Multi-body Dynamics， 2005， 219（4）： 345—355.
[12]Sanborn G G， Choi J， Choi J H. Curve-induced distortion of polynomial space curves， flat-mapped extension modeling， and their impact on ANCF thin-plate finite elements [J]. Multibody System Dynamics， 2011， 26（2）： 191—211.
[13]張 越， 赵 阳， 谭春林，等. ANCF索梁单元应变耦合问题与模型解耦[J]. 力学学报， 2016， 48（6）： 1406—1415.
Zhang Yue， Zhao Yan，Tan Chunlin， et al. The strain coupling problem and model decoupling of ANCF cable/beam element[J]. Chinese Journal of Theoretical and Applied Mechanics，2016，48（6）：1406—1415.
[14]Hyldahl P， Mikkola A M， Balling O， et al. Behavior of thin rectangular ANCF shell elements in various mesh configurations [J]. Nonlinear Dynamics， 2014， 78（2）： 1277—1291.
[15]刘桂荣，顾元通. 无网格法理论及程序设计[M].王建明，周学军，译. 济南：山东大学出版社，2007：45—133.
Liu G R， Gu Y T. An Introduction to Meshfree Methods and Their Programming[M]. Jinan：Shandong University Press， 2007：45—133.
[16]Chen J S， Pan C， Roque C M O L， et al. A lagrangian reproducing kernel particle method for metal forming analysis[J]. Computational Mechanics， 1998， 22（3）： 289—307.
[17]Chen J S， Pan C， Wu C T， et al. Reproducing kernel particle methods for large deformation analysis of non-linear structures[J]. Computer Methods in Applied Mechanics & Engineering， 1996， 139（1-4）： 195—227.
[18]Liu W K， Jun S， Zhang Y F. Reproducing kernel particle methods. Int J Numer Meth Eng[J]. International Journal for Numerical Methods in Fluids， 1995， 20（8-9）： 1081—1106.
[19]Xiong Y B， Long S Y. Local petrov-galerkin method for a thin plate[J]. Applied Mathematics and Mechanics， 2004， 25（2）： 210—218.
[20]Wang J G， Liu G R. A point interpolation meshless method based on radial basis functions[J]. International Journal for Numerical Methods in Engineering， 2002， 54（11）： 1623—1648.
[21]Cui X Y， Feng H， Li G Y， et al. A cell-based smoothed radial point interpolation method （CS-RPIM） for three-dimensional solids[J]. Engineering Analysis with Boundary Elements， 2015， 50： 474—485.
[22]Pilafkan R， Folkow P D， Darvizeh M， et al. Three dimensional frequency analysis of bidirectional functionally graded thick cylindrical shells using a radial point interpolation method （RPIM） [J]. European Journal of Mechanics / A Solids， 2013， 39（5）： 26—34.
[23]Gu Y T， Wang W L， Fu Q. An enriched radial point interpolation method （e-RPIM） for the analysis of crack tip[J]. Engineering Fracture Mechanics， 2010， 78（1）： 175—190.
[24]Chen S L， Li Y X. An efficient RPIM for simulating wave motions in saturated porous media[J]. International Journal of Solids & Structures， 2008， 45（25）： 6316—6332.
[25]Dai K Y， Liu G R， Han X， et al. Inelastic analysis of 2D solids using a weak-form RPIM based on deformation theory[J]. Computer Methods in Applied Mechanics & Engineering， 2006， 195（33）： 4179—4193.
[26]杜超凡， 章定国， 洪嘉振. 径向基点插值法在旋转柔性梁动力学中的应用[J]. 力学学报， 2015， 47（2）： 279—288.
Du Chaofan， Zhang Dingguo， Hong Jiazhen. A meshfree method based on radial point interpolation method for the dynamic analysis of rotating flexible beams[J]. Chinese Journal of Theoretical and Applied Mechanics， 2015，47（2）： 279—288.
[27]杜超凡， 章定国. 光滑节点插值法：计算固有频率下界值的新方法[J]. 力学学报， 2015， 47（5）： 839—847.
Du Chaofan， Zhang Dingguo. Node-based smoothed point interpolation method： a new method for computing lower bound of natural frequency[J]. Chinese Journal of Theoretical and Applied Mechanics， 2015， 47（5）： 839—847.
[28]杜超凡， 章定国. 基于无网格点插值法的旋转悬臂梁的动力学分析[J]. 物理学报， 2015，64（3）： 396—405.
Du Chaofan， Zhang Dingguo. A meshfree method based on point interpolation for dynamic analysis of rotating cantilever beams[J]. Acta Phys. Sin.， 2015， 64（3）： 396—405.
[29]章孝順，章定国，陈思佳，等. 基于绝对节点坐标法的大变形柔性梁几种动力学模型研究[J].物理学报， 2016， 65（9）：148—157.
Zhang Xiaoshun， Zhang Dingguo， Chen Sijia，et al. Several dynamic models of a large deformation flexible beam based on the absolute nodal coordinate formulation[J]. Acta Phys. Sin.， 2016， 65（9）：148—157.
[30]陈思佳， 章定国， 洪嘉振. 大变形旋转柔性梁的一种高次刚柔耦合动力学模型[J]. 力学学报， 2013， 45（2）：251—256.
Chen Sijia， Zhang Dingguo， Hong Jiazhen. A high-order rigid-flexible coupling model of a rotating flexible beam under large deformation[J]. Chinese Journal of Theoretical and Applied Mechanics， 2013， 45（2）：251—256.
[31]李 彬， 刘锦阳. 大变形柔性梁系统的绝对坐标方法[J]. 上海交通大学学报， 2005， 39（5）： 827—831.
Li Bin， Liu Jinyang. Application of absolute nodal coordination formulation in flexible beams with large deformation[J]. Journal of Shanghai Jiaotong University， 2005， 39（5）： 827—831.
[32]章孝順， 章定国， 洪嘉振. 考虑曲率纵向变形效应的大变形柔性梁刚柔耦合动力学建模与仿真[J]. 力学学报，2016，48（3）：692—701.
Zhang X H， Zhang D G， Hong J Z. Rigid-flexible couple dynamic modeling and simulation with the longitudinal deformation induced curvature effect for a rotating flexible beam under larger deformation [J]. Chinese Journal of Theoretical and Applied Mechanics，2016，48（3）：692—701.
[33]Amold M， Brüls O. Convergence of the generalized-α scheme for constrained mechanical systems[J]. Multibody System Dynamics， 2007，18（2）：185—202.
An absolute nodal coordinate formulation based on
radial point interpolation method for planar slender beams
CHEN Yuan-zhao， ZHANG Ding-guo， LI Liang
（School of Science， Nanjing University of Science and Technology， Nanjing 210094， China）
Abstract： In order to alleviate or eliminate the ‘distortion phenomenon of the deficient beam elements in the traditional absolute nodal coordinate formulation （ANCF）， an ANCF based on radial point interpolation method （RPIM） for flexible beams is proposed in which a new RPIM shape functions are constructed to describe the absolute configuration of flexible beams. For a pure bending cantilever beam （an arc）， there is always difference between configuration of the beams by using the gradient deficient beam elements and the actual configuration of the beam in the traditional ANCF. The difference， namely ‘distortion phenomenon， becomes more obvious with the increase of the bending deformation of the beam， which may cause the pseudo strain and have serious influence on accuracy of numerical solution. In the present method， the RPIM is used to discretize the deformation field through a set of field nodes and the shape functions are generally established based on field nodes within a support domain of the calculating point. The calculating points are generally located in the central region of the support domain， and the support domain of different calculating points can have more coincident parts. Thus， the connection between field nodes is strengthened， which makes the method describe the configuration in more reasonable and effectively alleviate influence of the distortion and pseudo strain. The simulation results show that the proposed method has higher calculation accuracy and efficiency and is more adaptive for the non-equidistant nodes compared with the traditional ANCF， which can be further extended in the dynamic field of flexible multi-body system.
Key words： multibody dynamics； flexible beam； absolute nodal coordinate formulation； radial point interpolation method； distortion
