引言

随着航天技术的发展，轻质、高精度、柔性零部件被广泛应用于众多工程领域，如空间站柔性操作机械臂、卫星环形桁架天线、航天器太阳帆板、空间机器人等.这些系统中存在刚性部件与柔性部件的大范围运动相互耦合，而且柔性部件的大范围运动和自身的大变形相互耦合^[1]，属于刚柔耦合多体系统动力学问题.基于小变形假设的柔性多体系统动力学是多刚体系统动力学与线性有限元的组合方法，无法描述柔性部件的大变形问题.因此，非线性有限元逐渐成为描述和预测这类柔性多体系统动力学响应的主流方法^[2].

1997年，美国多体动力学专家Shabana教授为了实现柔性体大范围运动和大变形的统一描述，提出了绝对节点坐标方法（Absolute Nodal Coordinate Formulation，ANCF）^[3].该方法基于连续介质力学的构型描述，在惯性坐标系中定义单元节点坐标，并采用节点斜率矢量代替传统有限元中的节点转角坐标，得到的柔性多体系统微分-代数方程具有常数质量矩阵、不存在科氏力和离心力项等特点，被认为是多体系统动力学研究历史上的重要进展之一，受到国内外众多学者的研究与关注.但该方法采用斜率矢量作为节点坐标，导致单元自由度多，并且存在严重的闭锁问题，影响了该方法的计算效率和精度.当涉及复杂几何构型时，采用绝对节点坐标方法单元进行网格划分会使得单元节点处的斜率矢量难以确定，造成网格的生成非常困难和耗时.此外，绝对节点坐标方法单元多选用Hermite多项式作为形函数，使得离散之后的网格模型是初始几何模型的近似，这种近似会降低动力学仿真结果的精度.因此，将绝对节点坐标方法单元集成到现有的多体系统动力学分析软件时遇到不少困难.

在柔性多体系统动力学研究中，采用绝对节点坐标方法对柔性构件进行建模时，需要人工划分网格，这是一个非常耗时的过程.此外，离散后的网格模型通常是初始几何模型的近似，对于研究大应变和接触碰撞等问题，这种近似会造成较大误差.美国能源部Sandia国家实验室曾做过一项统计，在汽车、航空航天和造船行业等领域，采用传统有限元进行分析，大约80%的时间用于网格剖分以及剖分前的几何模型准备工作^[4]，如图1所示.因此，在多柔体系统建模时，一个关键问题就是如何将几何描述模型和力学仿真计算模型进行统一，简化网格生成和细化的过程，提高数值仿真计算效率和精度，实现模型设计和分析一体化.

2005年，Hughes等^[5]为了解决这一问题，提出了等几何分析（Isogeometric Analysis，IGA）的概念，直接将构造CAD几何造型的样条函数[如Bézier、B-spline、非均匀有理B样条（Non-Uniform Rational B-Splines，NURBS）和T-spline等]及其控制点坐标作为有限元分析的输入信息，从而省去了网格划分过程，形成了几何与分析的直接联系.因此可以说，等几何分析是有限元法的扩展，成熟的CAD建模技术和有限元计算为等几何分析方法提供了坚实基础.由于该方法在设计与仿真阶段对模型的几何描述形式是完全一致的，因而被称为“等几何分析”.在等几何分析中，几何模型和分析网格模型采用同一表达方式，避免了网格细化过程中与几何模型数据的频繁交互，同时避免了传统有限元中网格划分的复杂过程.此外，和绝对节点坐标方法相比，等几何分析的结果精度更高，其分析使用的多项式直接来源于几何模型数据，避免了传统有限元中采用分段多项式逼近的方式而引入的误差.采用等几何分析方法建立的多柔体系统动力学模型不仅能精确描述柔性构件的几何特征，而且能高效描述柔性构件的大转动、大变形耦合动力学，还具有单元自由度少和收敛性高的特性，成为一种颇具发展前景的建模和计算方法.

根据是否考虑梁截面的剪切变形，可将基于等几何分析的梁单元分为两类:只考虑轴向伸长和横向弯曲变形的Euler-Bernoulli梁单元，以及考虑截面剪切和轴向转动效应的Timoshenko梁单元.Weeger等^[6]提出了基于等几何分析的细长梁单元，将其应用于对空间细长结构的非线性振动进行分析，并和传统拉格朗日插值得到的梁单元进行对比.分析结果表明，等几何梁单元表现出更好的收敛性，且不存在转角坐标，避免了转动插值的奇异性.Bouclier等^[7]提出了基于NURBS插值的等几何Timoshenko梁单元，并采用缩减积分和B投影技术解决了短粗梁结构分析中常见的剪切和薄膜闭锁问题.随后，Sang等^[8]通过研究发现，在等几何Timoshenko梁单元中，随着NURBS插值形函数次数的升高，剪切闭锁效应逐渐减弱，并通过对梁结构自由振动的数值分析验证了该结论.此外，Luu等^[9]通过对抛物线形和椭圆形空间曲梁结构进行动力学分析，并和同阶拉格朗日插值多项式构造的梁单元进行对比，计算结果表明等几何单元的精度更高.这是由于NURBS能够对此类结构的几何构型进行精确描述，而以拉格朗日插值多项式构造的梁单元是对该结构几何构型的近似.近年来Borkovicab等^[10]提出了考虑三维空间梁厚度方向非线性应变的几何精确欧拉伯努利梁单元; 陈明飞等^[11]基于一阶剪切变形理论并采用等几何有限元方法对任意曲率的功能梯度曲梁进行自由振动分析; 夏阳等^[12]应用拟协调有限元，使用降阶基函数逼近梁内应变项，解决了复杂三维曲梁结构仿真中的闭锁问题.

在等几何分析概念提出后，基于等几何分析的计算力学研究自然延伸到板壳结构.基于Kirchhoff-Love理论的薄板壳有限元要求近似位移场具有C1连续性，这是传统有限元构造薄壳单元的难点，但该条件在基于光滑样条理论的等几何分析中通常可以自然满足.Kiendl等^[13]基于等几何分析提出了一种薄壳单元，通过构造光滑的基函数保证应变和应力的全局连续性.由于其研究仅考虑壳单元节点的3个中面位移自由度，单元的自由度数非常少.基于该薄壳单元，Chen等^[14]进一步提出了等几何薄膜单元，采用显示时间积分策略对空间薄膜结构的充气过程进行了数值仿真，并对充气过程中薄膜结构出现的动力松弛现象进行了讨论.为了描述面外剪切变形，Benson等^[15]提出了基于实体退化的Reissner-Mindlin壳单元，通过减缩积分（reduced integration）技术提高了计算效率.他们基于该单元引入弹塑性本构模型，分析了结构中出现的塑性变形，并将提出的单元集成到LS-DYNA软件中.此外，Echter^[16]基于高阶剪切板理论提出了多种等几何板壳单元，并对单元的收敛性进行了系统的对比.

在动力学与控制领域，等几何分析方法也逐渐受到众多学者的关注.王悦等^[17]提出了一种基于T样条曲面的变网格柔性系统等几何分析方法，通过对受冲击柔性薄板的动力学分析表明所提出T样条单元及局部细化算法可以只在接触碰撞等应变剧烈变化的区域实现局部网格细化，从而控制系统自由度数，提高计算效率; 刘涛等^[18]针对表面黏贴有压电层的功能梯度板的动力学及主动振动控制问题，建立了一种基于三阶剪切变形理论的等几何分析求解方法; 陈龙等^[19]基于等几何分析方法研究了齿轮接触问题，研究结果表明该方法具有计算效率高、速度快、应力场更加光滑等优势; 荣吉利等^[20]为实现大变形结构动力学问题中的网格自适应算法，将多层贝塞尔提取方法用于HB样条建模分析以实现形函数的统一; Tobias等^[21]采用等几何分析模型研究了柔性多体系统的冲击响应，模拟了弹性球体对刚性表面的冲击和长弹性杆的冲击并进行了试验验证，表明等几何分析方法在研究柔性多体系统碰撞问题上具有高精度和高效率的特点.

在等几何分析中，B样条和NURBS是两类常用的样条基函数，对于二维或三维问题，B样条与NURBS通常采用张量积形式构造，局部加密将会导致全局网格的变化，无法实现局部细化.为实现网格局部细化及自适应网格，多种可局部加密的样条函数在计算机图形学领域被提出，例如，T样条、R样条（Local Refined Splines），以及层次B样条（Hierarchical B-splines，简称HB样条）^[22].此外，三角B样条^[22]采用三角形拓扑结构，天然具有实现网格局部细化的优势.虽然B样条和NURBS在网格局部细化上存在不足，但在CAD领域，根据国际标准化组织（ISO）颁布的工业产品数据交换标准STEP，NURBS是定义工业产品几何形状的唯一数学方法.因此NURBS在实现计算机辅助设计与计算机辅助工程的统一方面具有重要的地位.

本工作将等几何分析方法的研究范围扩展至含大转动、大变形的柔性多体系统动力学研究中，基于NURBS曲线和曲面分别提出了Euler-Bernoulli梁单元和Kirchhoff-Love壳单元，通过数值算例对单元的性能进行验证，表明提出的等几何单元不仅能精确描述柔性构件的几何特征，而且能精确描述柔性构件的大转动、大变形耦合动力学特性，还具有单元自由度少和收敛性高的特性，为柔性多体系统建模提供了一种高效、精确的新单元.

1 NURBS基本理论

NURBS是在1975年由Versprille^[23]在其博士学位论文中提出.1991年，国际标准化组织（ISO）颁布的工业产品数据交换标准STEP中，把NURBS作为定义工业产品几何形状的唯一数学方法.由于NURBS是B样条的有理形式，所以首先介绍B样条的构造.

1.1 节点矢量

B样条是在1972年由Riesenfeld^[24]提出的一种样条曲线特殊表示形式，在等几何分析中一般称为B样条基函数.B样条基函数是定义在一个节点矢量（knot vector）上的分段多项式，节点矢量的定义为

其中ξ_i（i=1，2，···，n+p+1）称为节点（knot），p为B样条基函数的次数（order），n为B样条基函数的个数.节点矢量表示的是一组非递减的非负实数，即0≤ξ_i≤ξ_i+1，节点矢量$\Xi$构成了定义B样条基函数的参数空间^[25].

节点矢量中节点的分布决定了B样条基函数的连续性.在节点矢量中相邻的两个或两个以上的节点值相等，称这些节点为重复节点，节点的重复i称为节点重复度（multiplicities），节点重复度决定了B样条基函数的连续性.没有重复节点时，B样条基函数在该节点具有C^p-1连续; 如果存在重复节点，基函数在该节点具有C^p-k连续，k为节点重复度.

在计算机图形学中，通常采用开放型（open）节点矢量定义B样条基函数，开放型节点矢量表示在节点矢量$\Xi$的两端节点ξ₁和ξ_n+p+1处，节点重复度为p+1.对于开放型节点矢量定义的B样条基函数，在参数区间端点处具有插值性，但在内部节点通常不具有插值性.下文不作特殊说明的情况下B样条和NURBS所采用的节点矢量均为开放型节点矢量.

1.2 B样条曲线和曲面

对于给定的节点矢量$\Xi$，B样条基函数由Coxde Boor递归方法^[26-27]定义，首先对于p=0，即0阶B样条可定义为

其次，对于p=1，2，···，n，B样条基函数定义为

由式（3）可知，一个p阶B样条基函数是两个p-1阶B样条基函数的线性组合，这为B样条基函数的计算提供了方法.

图2给出了定义在节点矢量$\Xi$ ={0，0，0，0.25，0.5，0.75，0.75，1，1，1}上的7个二次B样条基函数.从图中可以看出，该二阶B样条基函数在节点0.75处的节点重复度为2，因此在该节点处是C⁰连续; 而在节点0.25和0.5处的重复度为1，因而是C¹连续.

根据B样条基函数的定义，对于给定的节点矢量$\Xi$，p阶B样条基函数具有以下性质:

（1）非负性（non-negativity）.即$\forall \xi$，有N_i，p（ξ）≥0.由传统Lagrange和Hermite多项式构造的基函数则无法保证非负性，导致非协调元在界面处的几何不连续.此外，B样条基函数的非负性使得等几何单元的质量阵也是非负的.

（2）单位分解性（partition of unity）.即$\forall \xi$，均有$\sum _{i=1}^n{N}_{i，p}（\xi ）⩾0$，这是构造有限元基函数必须具备的性质.

（3）紧支性（compact support）.如果$\xi \notin \left[{\xi }_i，\right.\left.{\xi }_{i+p+1}\right)$，则N_i，p（ξ）=0，即B样条基函数N_i，p（ξ）仅在[ξ_i，ξ_i+p+1）上非零.

（4）可微性（differentiability）.B样条基函数N_i，p（ξ）在节点ξ_i处是p-m_i次可微的，其中m_i为节点ξ_i的重复度.B样条基函数的可微性保证了构造的样条曲线是全局光滑的，并且保证基于B样条的等几何分析单元在内部节点处应力和应变也是连续的.

由式（3）可知，给定B样条基函数的次数p和合适的节点矢量$\Xi$，就可以构造n个B样条基函数.同样，在空间中给定n个控制点，则p阶分段多项式B样条曲线可通过B样条基函数和相应控制点的线性组合表示

其中r为B样条曲线上任意一点的位置矢量，N_i，p是定义在节点矢量$\Xi$上的p次B样条基函数，B_i是空间中定义B样条曲线的控制点，由控制点B_i在空间中按下标顺序连接构成的多边形称为B样条曲线的控制多边形.图3给出了平面空间内的一个二次B样条曲线及其控制点和控制多变形示意图，对应的基函数如图2所示.

B样条曲面是张量积曲面，通常由两个节点矢量$\Xi$={ξ₁，ξ₂，···，ξ_n+p+1}和Γ={η₁，η₂，···，η_m+q+1}以及n×m个控制点B_i，j构成.N_i，p（ξ），i=1，2，···，n+p+1和M_j，q（η），j=1，2，···，m+q+1分别为定义在节点矢量$\Xi$和Γ上的一维B样条基函数.则B样条曲面可表示为

其中r（ξ，η）为B样条曲面上任意一点，其中p和q分别为ξ和η方向上B样条基函数的次数，由控制点B_i，j在空间中构成的多边形网格称为控制网格，定义B样条曲面的参数域为[ξ₁，ξ_n+p+1]×[η₁，η_m+q+1].

图4给出了一个双二次（p=2，q=2）B样条曲面，其中图4（a）为该B样条曲面的基函数示意图，图4（b）给出了该B样条曲面及其控制点和控制网格，其中两个参数方向的节点矢量分别为$\Xi$ ={0，0，0，0.5，1，1，1}和Γ={0，0，0，1/3，2/3，1，1，1}.

B样条曲面的基函数由单变量B样条基函数的张量积形式构成，则B样条曲面的基函数同样具有非负性、单位分解性和紧支性等性质.此外，B样条曲线和曲面还具有仿射变换（affine transformation）几何不变性，仿射变换是指平移、旋转和投影等线性变换.仿射变换几何不变性保证了B样条单元在刚性位移下的零应力特征.

1.3 B样条基函数的导数

在等几何分析中需要通过基函数的导数推导弹性力及其雅可比矩阵，因此本节简单介绍B样条基函数的求导方法.由式（3）可知，一个p阶B样条基函数是两个p-1阶B样条基函数的线性组合，则B样条基函数的导数可由低阶基函数组合递推得到.已知B样条的阶次p和节点矢量$\Xi$ ={ξ₁，ξ₂，···，ξ_n+p+1}，则p阶B样条基函数的一阶导数为

由式（7）依次递推得

其中系数α_i，j满足以下关系

式（7）给出了对B样条基函数求导的通用公式.需要注意的是:求导次数k不应超过p，对于B样条基函数，所有高于p阶的导数均为零; 在求导过程中，某些项的分母中节点之差可能为零，规定这种情况商为零，即0/0=0.

1.4 NURBS曲线和曲面

B样条虽然有强大的曲线曲面表示能力，但是对于许多在工程设计中常见的几何形状，如圆锥曲线、旋转面和圆锥曲面等解析曲面却无法精确表示.为解决这一问题，工程界通常采用NURBS曲线和曲面，下面对NURBS进行简单介绍.

NURBS是B样条的有理形式，则R^d空间的NURBS可以由R^d+1空间的B样条进行投影变换（projective transformation）得到.例如，R²空间内的圆弧曲线是R³空间的二次B样条曲线在平面内投影所得，如图5所示.

通过B样条的投影变换可以得到有理多项式C_R（ξ）=f（ξ）/g（ξ），其中f（ξ）和g（ξ）是分段多项式，下面首先介绍在R^d空间NURBS曲线的构造策略.设{${\mathit{B}}_i^w$}为R^d+1空间的B样条曲线的控制点集合，节点矢量为$\Xi$，则R^d内的NURBS控制点可由R^d+1空间中的B样条控制点进行投影得到，即

其中（B_i）_j为控制点B_i的第j个分量，w_i为控制点B_i对应的权重因子.在图5中，权重即为中B样条曲线控制顶点的纵坐标.由此可得NURBS基函数和NURBS曲线的定义为

其中，$R_i^p（\xi ）$为NURBS基函数.和B样条曲面类似，NURBS曲面同样由张量积的NURBS基函数构造，可表示为

其中r为NURBS曲面上的任意一点，为NURBS曲面的基函数.图6给出了由张量积NURBS构造的三维曲面，如球壳、圆环壳和旋转曲面等.对于这些圆锥曲面和旋转曲面，可由双二次NURBS曲面精确构造.NURBS是B样条基函数的有理形式，NURBS除了具有B样条的特性以外，还具有以下特殊性质:

（1）当所有权重因子相等时，NURBS将退化为B样条，也就是说B样条是NURBS的一种特例.

（2）NURBS基函数是有理多项式，与传统Lagrange和Hermite多项式一样，NURBS基函数的计算也是稳定的.

（3）NURBS可以精确表示标准的解析（如圆锥曲线曲面、旋转面）和自由曲线、曲面，为产品几何模型的定义提供了标准.

（4）为了修改曲线曲面的几何形状，既可通过调整控制点的位置，又可通过改变权重因子的值来实现，具有很大灵活性.

NURBS和B样条可通过h细分（节点插入），p细分（基函数升阶）和k细分（基函数升阶+节点插入）实现控制点个数的增加和基函数次数的升高，并且这些细分方式都不会改变样条的几何性质.这些细分方法为基于NURBS和B样条的等几何分析提供了丰富的网格细化方法.

2 基于NURBS的等几何单元

在结构力学中，广泛应用的梁理论有两种:Euler-Bernoulli梁理论和Timoshenko梁理论，Euler-Bernoulli梁理论也称为工程梁理论.其中Euler-Bernoulli梁理论的运动学假设是:梁单元在变形过程中，梁的中线法平面一直保持平面，并与梁的中线直垂.这个假定认为弯曲变形是梁的主要变形，剪切变形是次要的变形，因而可以忽略不计，这对于截面尺度远小于轴向尺寸的实腹梁来说，不会引起显著的误差，因而在工程界得到广泛应用.

2.1 Euler-Bernoulli梁单元

如图7所示，采用NURBS曲线对细长梁的中线进行建模，由NURBS曲线的定义可知，梁单元的中线上任意一点P的全局位置矢量可表示为

其中S为参数域中任意一点ξ∈[ξ_i，ξ_i+1）的非零基函数，构成该等几何梁单元的形函数.根据NURBS基函数的紧支性，对于p阶NURBS基函数，点处有且仅有p+1个非零基函数，可表示为

其中I₃为3×3的单位矩阵，e为p+1个非零基函数所对应的控制点构成的单元节点坐标，可表示为

在完全拉格朗日格式（Total Lagrangian Formulation）框架下，不考虑梁单元的截面剪切变形，则梁单元上任意一点的Green应变可解耦为拉伸应变ε^l与弯曲应变ε^κ两部分

其中拉伸应变ε^l可表示为

弯曲应变可表示为

弯曲应变中的表示为材料曲率（material measure of curvature）$\stackrel{-}{\kappa }（\xi ）=\left|{\mathit{r}}_{\xi }\times {\mathit{r}}_{\xi \xi }\right|/\left|{\mathit{r}}_{\xi }\right|，{\stackrel{-}{\kappa }}_0（\xi ）=\left|{\left({\mathit{r}}_0\right)}_{\xi }\times {\left({\mathit{r}}_0\right)}_{\xi \xi }\right|/\mid {\left({\mathit{r}}_0\right)}_{\xi }$.

在得到梁单元上任一点的应变后，在梁单元的初始构形上进行积分，则单元应变能可表示为

其中E为材料的杨氏模量，V^e为初始构型中曲梁单元的体积.从式（21）可知，单元应变能也可分为拉伸应变能与弯曲应变能两部分，其中拉伸应变能可表示为

弯曲应变能可进一步展开为

单元弹性力可表示为应变能对广义坐标的一次偏导数

其中拉伸弹性力与弯曲弹性力分别为

单元弹性力雅可比矩阵可表示为应变能对广义坐标的二次偏导数，具体推导过程参考文献^[28].根据虚功原理，梁单元的质量矩阵可表示为

上式表明，NURBS描述的梁单元质量矩阵也为常数矩阵.此外，若忽略梁单元中的弯曲应变能，只考虑拉伸应变能，则该梁单元退化成索单元.

2.2 Kirchhoff-Love壳单元

在许多工程构件的模拟中，壳单元都是极为重要的.其中薄壳出现在许多产品中，诸如汽车上的金属薄板，飞机的机舱、机翼和方向舵，以及某些产品的外壳，如火箭、手机、洗衣机和计算机等.用实体单元模拟这些构件通常需要大量单元，同时存在闭锁问题，从而导致非常昂贵的计算费用.因此采用壳单元模拟薄壁结构，能够极大地改善运算效率.本节根据Kirchhoff-Love假设，采用NURBS曲面对薄壳的中面进行描述，提出了适合对薄壳结构进行建模的等几何分析壳单元.

如图8所示，在壳单元的初始构型中，基于Kirchhoff-Love假设，壳体上任一点$P_0^z（\xi ，\eta ，\zeta ）$的位置矢量可通过壳的中面上对应点P₀（ξ，η）的位置矢量r₀以及该点的法线方向n₀确定

其中‘0’表示初始构型，ξ和η为薄壳的中面对流曲线坐标（convective curvilinear coordinates），ξ（-h/2≤ξ≤h/2）是点与点之间的距离，h是壳的厚度.r₀是由NURBS曲面描述的壳中面上点P₀（ξ，η）的位置矢量，即

其中$R_{ij}^{pq}（\xi ，\eta ）$和B_i，j分别为NURBS曲面基函数和控制点.

为了精确描述薄壳单元的变形，在薄壳单元初始构形中的任一点${\mathrm{P}}_0^{\mathrm{z}}（\mathrm{\xi }，\mathrm{\eta }，\mathrm{\zeta }）$处，建立两个局部坐标系${\left({\mathbf{g}}_0^z\right)}_1-{\left({\mathbf{g}}_0^z\right)}_2-{\left({\mathbf{g}}_0^z\right)}_3$和${\left({\mathbf{e}}_0^z\right)}_1-{\left({\mathbf{e}}_0^z\right)}_2-{\left({\mathbf{e}}_0^z\right)}_3$.类似地，在当前构型中的点$P^z（\xi ，\eta ，\zeta ）$处也可建立两个局部坐标系${\left({\mathbf{g}}^z\right)}_1-{\left({\mathbf{g}}^z\right)}_2-{\left({\mathbf{g}}^z\right)}_3$和${\left({\mathbf{e}}_{}^z\right)}_1-{\left({\mathbf{e}}_{}^z\right)}_2-{\left({\mathbf{e}}_{}^z\right)}_3$.这些坐标系的详细定义读者可参考文献^[29].

基于连续介质力学理论^[30]，连续体壳单元中任意一点的变形梯度矩阵F表示为

其中dX₀和dx₀分别是在参考构型和当前构型中，点$P_0^z（\xi ，\eta ，\zeta ）$在切平面上的微元.变形梯度矩阵F能够通过一个正交矩阵R和一个非奇异对称矩阵U进行分解^[21]，即

其中R表示微元的纯刚体转动，R矩阵能够通过右伸长张量U从变形梯度矩阵F中进行提取.矩阵R和U可进一步表示为

其中，R为局部笛卡儿坐标${\left({\mathbf{e}}^z\right)}_1-{\left({\mathbf{e}}^z\right)}_2-{\left({\mathbf{e}}^z\right)}_3$到坐标系${\left({\mathbf{e}}_0^z\right)}_1-{\left({\mathbf{e}}_0^z\right)}_2-{\left({\mathbf{e}}_0^z\right)}_3$的转换.T₀表示初始构型中局部笛卡儿坐标${\left({\mathbf{e}}_0^z\right)}_1-{\left({\mathbf{e}}_0^z\right)}_2-{\left({\mathbf{e}}_0^z\right)}_3$到坐标系${\left({\mathbf{g}}_0^z\right)}_1-{\left({\mathbf{g}}_0^z\right)}_2-{\left({\mathbf{g}}_0^z\right)}_3$的转换.由此可知矩阵T₀可进一步表示为

类似地，在当前构型中，矩阵T也可进一步表示为

由向量${\left({\mathbf{e}}_0^z\right)}_1$和${\left({\mathbf{e}}_0^z\right)}_2$之间的正交性可知，薄壳单元的应变张量可表示为

其中${\left(g_0^z\right)}_{\alpha \beta }={\left({\mathbf{g}}_0^z\right)}_{\alpha }\cdot {\left({\mathbf{g}}_0^z\right)}_{\beta }$是壳单元中面的第一基本量^[32].因此，Green应变张量ε能够进一步分解为

其中ε^mem是壳单元的薄膜应变，ε^bend是弯曲应变.弯曲应变可表示为

其中κ₀和κ分别为壳单元的初始曲率和变形之后的曲率，可分别表示为

基于St. Venant-Kirchhoff假设，采用Voigt符号，第二类Piola-Kirchhoff应力张量σ_S可表示为

其中E是材料的杨氏模量，ν是泊松比.ε₁₁，ε₁₂和ε₂₂是Green应变张量的三个分量.

将式（36）代入式（39），并且在厚度方向ζ（-h/2≤ζ≤h/2）进行积分，则第二类Piola-Kirchhoff应力可解耦成薄膜应力σ_n和弯曲应力σ_m两部分，分别表示为

和

式（40）和式（41）中的${\epsilon }_{11}^{\mathrm{m}\mathrm{e}\mathrm{m}}$，${\epsilon }_{12}^{\mathrm{m}\mathrm{e}\mathrm{m}}$，${\epsilon }_{22}^{\mathrm{m}\mathrm{e}\mathrm{m}}$，${\epsilon }_{11}^{\text{bend}}$，${\epsilon }_{12}^{\text{bend}}$，${\epsilon }_{22}^{\text{bend}}$分别为薄膜应力σ_n和弯曲应力σ_m的分量.

针对线弹性材料，由广义胡克定律可知，单元应变能可表示为

单元弹性力可表示为应变能对广义坐标的一次偏导数

对式（43）进一步求导可得单元弹性力的雅各比矩阵

3 刚柔耦合系统动力学方程建立与求解

按照有限元方法中单元的组装过程，将柔性多体系统中各单元的质量矩阵、弹性力以及外力向量等进行组装，同时考虑系统约束条件，通过第一类拉格朗日方程可以建立描述刚柔耦合系统动力学特性的指标-3（Index-3）微分-代数方程组（DAEs）

其中M为系统的常数质量矩阵，q为系统广义坐标，F（q）、$\mathit{q}（\mathit{q}，\dot{\mathit{q}}）$分别为系统弹性力以及外广义外力向量，Φ（q，t）、Φ_q分别为系统约束方程及其对广义坐标q的雅可比矩阵，λ为拉格朗日乘子向量.

对于上述柔性多体系统动力学方程，通常可采用隐式时间积分算法或者显示时间积分算法求解.鉴于显式时间积分算法的条件收敛特性会限制时间积分步长，因此对于柔性多体系统动力学问题，特别是具有大变形和大范围运动的多柔体系统动力学问题，多采用隐式时间积分算法进行求解.其中在Newmark方法基础上发展而来的广义-α方法受到广泛关注.

广义-α方法将指标-3的方程（45）经过差分直接离散成代数方程进行求解，其迭代过程如下

其中，h为时间积分步长，β和γ为决定计算精度与效率的算法参数.式（46）中矢量参数a满足如下关系

上式中各参数的选取方法如下

其中$\tilde{\rho }$∈[0，1]为算法的谱半径，决定着算法能量耗散分布的频率范围. $\tilde{\rho }$取值越大，则算法产生的能量耗散越小.当$\tilde{\rho }$=1时，广义-α方法将保持系统能量而不产生数值耗散; 当$\tilde{\rho }$=0时，则产生最大的能量耗散.

通过广义-α方法对动力学方程进行离散后，需要在每个时间步长内采用Newton-Rapson迭代求解如下具有稀疏矩阵的大规模线性代数方程组

式中

$\widehat{\beta }$，$\widehat{\gamma }$为满足如下关系的算法参数

4 数值算例

4.1 半圆形悬臂梁弯曲静力学分析

对图9所示的半圆形悬臂梁进行弯曲静力学分析.为了对比研究，本算例采用和文献^[34]一致的参数进行建模.在初始构型中，该悬臂梁曲率半径R=0.5m，截面半径r=0.0173m，材料的杨氏模量E=2.1×10¹⁰Pa.悬臂梁在O点固支，若在悬臂梁的自由端A点施加如图方向所示的外力矩M=λπEI/l，由材料力学可知，当仅发生纯弯曲变形，且λ=1时，该悬臂梁最终将形成一个整圆.本算例分别采用二阶和三阶NURBS梁单元以及ANCF缩减梁单元对该悬臂梁进行离散.

图10中给出了采用128个三阶NURBS梁单元进行离散，悬臂曲梁在不同外力矩作用下的弯曲构形和Mises应力分布图.从图中可看出，随着外力矩的增大，悬臂梁所受到的Mises应力也逐渐增大，并且应力在悬臂梁上的分布是很均匀的.当λ=1时，悬臂梁几乎弯曲变形为一个整圆.

为了和ANCF缩减梁单元的计算结果进行对比，以悬臂曲梁的固支点O和末端点A之间的距离来衡量单元的收敛性，得到如图11所示的误差收敛曲线图.横坐标为总自由度，纵坐标为平衡之后OA两点之间的距离.从图中可以看出，NURBS梁单元的收敛性要优于ANCF缩减梁单元，这是因为NURBS采用有理多项式进行插值，而ANCF是采用三次Hermite多项式插值，因而在对半圆形悬臂梁的初始几何描述上NURBS单元更加精确.此外，在自由度相同的情况下，和二阶NURBS梁相比，三阶NURBS单元可以得到更加精确的结果.

4.2 平面四连杆机构动力学分析

如图12所示，本节通过由一个曲柄，一个摇杆和一个连杆构成的平面四连杆机构检验提出的NURBS梁单元的动力学特性.

连杆和曲柄以及摇杆之间采用球铰连接，曲柄和摇杆的另一个端点通过转动副固定在地面上.该四连杆机构通过作用在曲柄上的驱动力矩M进行转动，驱动力矩随时间的变化如下式所示

该模型是一个典型的多体系统动力学模型.其中各构件的几何和材料参数如表1所示.

为了进行单元的收敛性测试，分别采用68个（曲柄:12，连杆:32，摇杆:24）三阶NURBS梁单元对该四连杆机构进行离散，同时采用68个和136个ANCF缩减梁单元进行动力学计算，得到的连杆中点在空间中的位置随时间变化曲线如图13所示.从图中可以看出68个三阶NURBS梁单元计算结果和136个ANCF缩减梁单元计算结果吻合良好，当采用68个ANCF单元离散时计算结果并没有收敛.

为了对比NURBS梁单元的计算效率，本算例的计算在一台16核（主频为2.90GHz）、64GB内存的工作站上进行，其中计算步长取为5.0×10^-4s，容许误差设为1.0×10^-6.表2给出了划分不同单元数目采用ANCF单元和NURBS单元进行离散的模型所消耗的CPU时间.从表中可看出，计算所需时间随着单元数目增多而增加.当采用相同单元数目进行离散，NURBS消耗的CPU时间比ANCF消耗的CPU时间要少，当单元足够多时，这种优势更加明显.

通过以上的对比可以发现NURBS单元的收敛性明显优于ANCF缩减梁单元，同时NURBS单元的计算效率也更高.这是因为采用三阶NURBS离散，单元之间是C2连续的，而ANCF缩减梁单元之间只能保证C1连续，造成单元节点处弯曲应变的不连续.此外，NURBS梁单元的节点坐标只有控制点，而ANCF缩减梁单元的节点坐标除了节点位置以外还有节点处的斜率矢量，这导致在相同单元总数情况下，采用NURBS离散系统自由度更少，因而计算效率更高.该系统采用136个单元离散时，NURBS单元建模求解所消耗的CPU时间仅为ANCF建模的48.7%.在NURBS单元建模中，NURBS插值、单元组装和动力学方程求解三部分消耗的CPU时间占比约为13.5%、18.7%和67.8%.

最后给出了四连杆机构的在驱动力矩作用下8个不同时刻的构型图，如图14所示.从图中可以看出，四连杆结构在转动过程中连杆发生了很大变形，在t＞0.5s之后，驱动力矩逐渐衰减到0，可以观察到曲柄发生了回弹.

4.3 圆柱壳在集中力作用下的静力学分析

本节采用文献^[18]中圆柱壳的拉伸静力学算例来验证提出的NURBS薄壳单元的正确性，并将结果与该文献进行比较.如图15所示，一个无约束的圆柱壳在A、E两点受到一对大小相等、方向相反的集中力F=40kN的作用，其中A点为圆柱体上母线的中点，E点为下母线的中点.圆柱壳的长度L为10.35mm，内径R与厚度h分别为4.953mm与0.094mm，材料杨氏模量E=10.5×10⁶MPa，泊松比υ=0.3125.

考虑到圆柱壳模型的对称性，本节仅对该模型上半部分的四分之一（ABCD）进行建模，同时在边AB、AD和BC上施加对称边界条件.现将在圆柱壳上施加外力F的过程平均分成10个增量步，并设静力学迭代容许误差为1.0×10^-6.图16中给出了分别采用8×12、16×24以及32×48的三次NURBS对八分之一圆柱壳进行离散时，壳体上A、B和C三点的位移随外力增加的变化曲线.从图中可看出，采用16×24的网格离散可得到收敛的结果.

当集中力F=40kN时，壳的变形达到最大，为了进行收敛性的对比，给出采用相同单元数进行离散，三阶NURBS和ANCF缩减壳单元计算结果，如表3所示.将有限元软件ABAQUS采用64×96的S4R单元的计算结果作为参考值，从该表可以看出NURBS单元的收敛性明显优于ANCF.

图17给出了圆柱壳在不同外力作用下的变形构形图.从该图可以看出，圆柱壳的变形可明显地分为两部分.当拉力小于20kN时，圆柱壳发生以弯曲变形为主的大变形.当拉力继续增大，圆柱壳几乎不再增加Z方向的位移，而以发生拉伸变形为主，并发生了一定程度的屈曲变形.

4.4 矩形单摆的动力学分析

本节通过由矩形薄板的单摆模型检验NURBS薄壳单元的动力学特性，如图18所示，O点采用球铰进行固定.其中材料杨氏模量E=1.0×109Pa，泊松比ν=0.3，密度ρ=1800 kg/m³.薄板的几何尺寸为:a×b×h=0.3×0.3×0.001m³.本节研究该单摆模型在重力作用下的动力学响应，其中重力加速度为G=9.81 m/s².

为了进行单元的收敛性测试，分别采用三阶NURBS单元将薄板离散为4×4、8×8、16×16和32×32个单元.图19中给出了不同单元数目下单摆模型中的A点沿Z方向的位移时间历程曲线.从图中可看出，前0.3s内不同单元离散模型计算得到的动力学响应类似.随后，单摆变形进一步增大，需要离散为16×16个单元才能得到收敛的数值结果.图20给出了4个不同时刻双摆的构形图.此外给出了三阶NURBS单元和ANCF缩减壳单元的计算效率的对比，如表4所示.当单元总数相同时，采用三阶NURBS壳单元离散时系统自由度大大小于ANCF缩减壳单元，因而计算效率得以提高.

5 结论

本文基于NURBS曲线和张量积NURBS曲面在等几何分析框架下分别提出了Euler-Bernoulli梁单元和Kirchhoff-Love壳单元，为柔性多体系统细长梁和薄板壳结构提供了两类精确、高效的建模方法.同时在完全拉格朗日格式下，根据Green应变张量对单元变形进行描述，结合第二类Piola-Kirchhoff应力张量给出单元应变能公式，从而推导了单元的弹性力和弹性力雅可比矩阵表达式.最后通过静力学及动力学数值算例对提出的两类单元的性能进行对比和验证.表明提出的等几何单元不仅能精确描述柔性构件的几何特征，而且能精确描述柔性构件的大转动、大变形耦合动力学特性，还具有单元自由度少和收敛性高的特性，为柔性多体系统建模提供了一种高效、精确的新单元.

参考文献