多体系统动力学综述

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 []|,|,|,|,

T

x x x l x l e e e e e e e e e r e r e r e r ========= 1

2

1

2

304078,,,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 表示任意一点在局部坐标系下的纵向坐标，单元节点坐标矢量为：

1234567891011121

212110220304050601

21

2

71829101112e []|,|,,,,,|,|,,,,T

x 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 x

x 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所示。此模型下任意节点在全局坐标系下的插值函数以及单元坐标矢量如下所示：

2210123452220123451234567891011121211022030401

2

51/262/27/28/2

r Se e []|,|,,|,|,,T

x 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 y

r r e r e r e e y y ========⎡⎤+++++⎡⎤===⎢⎥⎢⎥+++++⎣⎦⎣⎦

=∂∂====∂∂∂∂====∂∂

1

2

911021112|,|,,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节点单元的单元节点坐标矢量为：