多体系统动力学综述

1. 绝对节点坐标法传统有限元方法建立的单元为非等参数单元，其使用节点处的位移梯度来描述物体的无限小的转动，但在物体发生大变形时，节点处的位移梯度已不能准确描述物体的转动变形，从而极大影响到计算的精度。

Shabana [1]提出了绝对节点坐标法（Absolute nodal coordinate formulation, ANCF ），其理论基础主要是有限元和连续介质力学理论。

该方法将物体的单元节点坐标定义在全局坐标系下，使用节点处的斜率(slope)矢量作为节点坐标而不是节点处的无限小转动[2]，不需要另外计算刚体位移与柔性变形之间的耦合，能较精确地计算大变形的多体系统动力学问题。

其最终推导出的多体系统的微分代数方程组（DAEs ）中，质量矩阵是一个常数矩阵，但刚度矩阵将是一个非线性的时间函数。

1.1梁单元的绝对节点坐标法Shabana 首先推导出一维梁单元的绝对节点坐标法模型[1][3]。

在这种模型中，梁单元用中性轴来简化，如图1所示，其上面任意一点P 在全局坐标系下的坐标表达为：23101232320123r =Se r a a x a x a x r b b x b x b x ⎡⎤+++⎡⎤==⎢⎥⎢⎥+++⎣⎦⎣⎦图1其中，x 为沿轴线的单元局部坐标，[]0,x l ∈，l 为梁单元初始长度；S 为单元形函数；e 为含有8个单元节点坐标的广义坐标矢量。

123456781102205162e []|,|,|,|,Tx x x l x l e e e e e e e e e r e r e r e r ========= 1212304078,,,x x x l x l r r r r e e e e x x x x ====∂∂∂∂====∂∂∂∂最终，通过绝对节点坐标法得到的无约束的单元动力学方程为：k e Me+Q =Q 其中，M 为常数质量矩阵，Q k 为广义弹性力矩阵，Q e 为广义外力矩阵。

由于一维梁单元模型无法考虑到梁的剪切变形，Omar 和Shabana [4]接着又提出了一种二维的考虑剪切变形的梁单元的绝对节点坐标法模型。

此模型中，单元上任意一点在全局坐标系下的位置坐标为：231012345232012345r Se r a a x a y a xy a x a x r b b x b y b xy b x b x ⎡⎤+++++⎡⎤===⎢⎥⎢⎥+++++⎣⎦⎣⎦其中，x 为梁单元上任一点在局部坐标系下相对于中性轴的横向坐标，y 表示任意一点在局部坐标系下的纵向坐标，单元节点坐标矢量为：123456789101112121211022030405060121271829101112e []|,|,,,,,|,|,,,,Tx x x x x x x l x l x l x l x l x l e e e e e e e e e e e e r r r r e r e r e e e e xx y y r r r r e r e r e e e e x x y y=============∂∂∂∂======∂∂∂∂∂∂∂∂======∂∂∂∂ 由上式可以看出，梁上任意一点的横向斜率是二次插值多项式，而任意一点的纵向斜率坐标却是一次多项式，变形过程中，横向应力与纵向应力相互耦合，进而在建模时很容易产生剪切自锁（shear-locking ）以及泊松自锁（Possion’s locking ）等问题；此外，这种ANCF 模型还存在收敛速度低的问题[5][6]。

基于此，很多学者做出了很多研究工作，提出了一些新的求解方法和模型。

Dufva 和Mikkola 等[7]改变单元的运动学描述公式，提出了一种更精确、更简便的平面剪切梁单元。

而后，Garc´ıa -Vallejo 和Mikkola 等[5]重新定义了单元上任意一点在全局坐标系下的插值函数，去掉了单元节点与单元中心线相切的斜率坐标，同时在单元重点再增加一个节点，提出了一种三节点的二维剪切梁单元模型，如图2所示。

此模型下任意节点在全局坐标系下的插值函数以及单元坐标矢量如下所示：2210123452220123451234567891011121211022030401251/262/27/28/2r Se e []|,|,,|,|,,Tx x x x x l x l x l x l r a a x a y a xy a x a x y r b b x b y b xy b x b x y e e e e e e e e e e e e r r e r e r e e y yr r e r e r e e y y ========⎡⎤+++++⎡⎤===⎢⎥⎢⎥+++++⎣⎦⎣⎦=∂∂====∂∂∂∂====∂∂12911021112|,|,,x l x l x l x l r r e r e r e e y y ====∂∂====∂∂图2结合减缩积分的方法，计算效率和精度可以大大提高，同时能有效解决原先的二维剪切梁单元面临的一些问题。

由于传统的基于绝对节点坐标法的一维梁单元很难准确的反映扭转以及剪切变形等的影响，Shabana 等[8][9]进一步提出了三维梁单元的绝对节点坐标法模型。

这种模型能放松Euler-Bernoulli 梁理论以及Timoshenko 梁理论关于梁变形过程中截面为刚性的假设，较好的反映梁变形过程中转动惯量、剪切变形以及扭转的影响。

Shabana 提出了两种基于绝对节点坐标法的三维梁单元，分别为2节点单元（图3）和4节点单元（图4）。

但无论哪种单元，梁上面任意一点在全局坐标系下的插值多项式均为：231012345672320123456723301234567r Se r a a x a y a z a xy a xz a x a x r b b x b y b z b xy b xz b x b x r c c x c y c z c xy c xz c x c x ⎡⎤+++++++⎡⎤⎢⎥⎢⎥==+++++++=⎢⎥⎢⎥⎢⎥+++++++⎢⎥⎣⎦⎣⎦且单元节点坐标都为24个。

图3图4其中，2节点单元的单元节点坐标矢量为：,,,,,,e=[e e ]e =[r r r r ],,r r r r =,r =,r =T TT T T T T T A A j x y z x y z j A Bx y z=∂∂∂∂∂∂ 4节点单元的单元节点坐标矢量为：3,2,2,1,3,3,2,1,2,e=[r r v r v r ],v TA A C C D DB B y z y z x z x y z r r r r r r r r r ϑϑϑ⎡⎤-⎡⎤⎢⎥⎢⎥=-=⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦通过数值结果模拟可以发现，这两种单元的结果是相同的。

然而，上述的三维梁单元模型中，在用虚功原理计算梁的广义弹性力时，是基于Green-Lagrange 应变张量与第二Piola-Kirchhoff 应力张量想结合的表达方式，最终计算效率很低。

Dufva 和Mikkola [10]等引入梁截面坐标系(Cross-section frame)和切线坐标系(Tangent frame)的概念，将梁上任意一点在全局坐标系中的坐标矢量表达为：z yy y z z 0t s t s r =r R A r +R A r +其中，矩阵R 表示截面坐标系与切线坐标系之间的旋转矩阵，A t 表示截面坐标系与全局坐标系之间的转移矩阵。

0r 是在全局坐标系下的梁中心轴上一点的坐标矢量，r s 为截面坐标系中的坐标矢量。

基于上述的2节点三维梁单元，最终构造出一种更高效率的基于绝对节点坐标法的三维梁单元模型。

Sugiyama 等[11]进一步对传统的基于ANCF 的三维梁单元的自锁问题作了研究，并提出了一种基于ANCF 的初始弯曲梁单元。

Sugiyama 同样指出，传统的ANCF 梁单元[4][8][9]利用Green-Lagrange 应变张量(1-2T J J I ε=（）)定义单元的变形，由于梁截面的变形，会导致高度耦合的变形模式，进而引起自锁问题。

而为了消除这种耦合的变形模式，将应变分量定义为沿着梁中心线的线性部分和弯曲/扭转两个部分。

为了避免自锁问题，利用Hellinger –Reissner 变分原理修正沿着梁中心线的剪切应力分布，同时使用假设应变场来减缓由于截面变形所引起的自锁问题。

1.2板单元和其他单元的绝对节点坐标法Shabana [12]回顾了传统的薄板的动力学研究方法，主要可分为浮动节点坐标法、增量有限元方法以及大旋转矢量方法。

提出了基于绝对节点坐标法的板单元的研究思路。

而后，Mikkola 和Shabana [13]进一步指出了传统的有限元板单元研究的缺陷和问题。

在增量有限元方法中，使用无限小的转动来定义板单元的运动，是一种非等参数单元，会导致刚体运动方程的线性化，影响运算精度；而大转动矢量方法则无法保证节点位移梯度的连续性，运算十分复杂。

并深入研究了基于绝对节点坐标法的板单元。

这种板单元是一种等参数单元，其质量矩阵是一个常数矩阵，离心力和科氏力为零。

图5为一个板单元的示意图。

基于ANCF 的板单元有四个节点，每个节点有12个坐标，因此一个板单元有48个坐标。

最终，单元上任意一点在全局坐标系下的坐标为：r S(,,)e x y z =其中，单元形函数S 可以分为两类，S A 利用不完全的四阶插值多项式推得，它不能保证单元表面节点的连续性；S B 能够更好的反映单元收敛性的要求，并能更准确的描述刚体运动和常数应变，通过数值仿真，这两种形函数的运行结果相互符合的很好。

e 为单元坐标矢量，r r r e r e e e e e TT T T j j j T j j T T T T A B C D x y z ⎡⎤∂∂∂=⎢⎥∂∂∂⎣⎦⎡⎤=⎣⎦图5然而，对于很薄的板(thin plate)，其沿着厚度方向的单元变形很小，可以忽略。

基于此，Dufva 等[14]提出了一种新的基于绝对节点坐标法的薄板单元。

这种单元省略了位置矢量对z 方向的矢量梯度坐标，最终形成了一个36个自由度的减阶的板单元，每个节点有9个坐标：r r e r e e e e e TT T j j Tj j T T T T A B C D x y ⎡⎤∂∂=⎢⎥∂∂⎣⎦⎡⎤=⎣⎦ 相比较此前的完全参数的基于ANCF 的板单元模型，这种新的减阶板单元的计算效率可提高100倍以上，且可以使用更少的单元数量。

Dmitrochenko 等[15]研究了16自由度的Hermitian 矩形板单元，在这中单元中，每个节点有四个自由度，分别为垂直位移i z ，两个斜率坐标()'/xi i z z x =∂∂和()'/yi i z z y =∂∂以及一个二阶斜率坐标()''2/xyi iz z x y =∂∂∂，如图6所示。