刘铖^,2, 胡海岩北京理工大学宇航学院力学系, 北京 100081

DYNAMIC MODELING AND COMPUTATION FOR FLEXIBLE MULTIBODY SYSTEMS BASED ON THE LOCAL FRAME OF LIE GROUP ¹⁾

Liu Cheng^,2, Hu HaiyanDepartment of Mechanics, School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China

通讯作者: 2)刘铖，副教授，主要研究方向：柔性多体系统动力学建模与计算；非线性有限元方法等. E-mail:liucheng_bit@aliyun.com

收稿日期:2020-08-17接受日期:2020-10-13网络出版日期:2021-12-31

基金资助:

1)国家自然科学基金资助项目.11672034 国家自然科学基金资助项目.11832005 国家自然科学基金资助项目.12072026

Received:2020-08-17Accepted:2020-10-13Online:2021-12-31
作者简介 About authors


摘要
多柔体系统动力学主要研究由多个具有运动学约束、存在大范围相对运动的柔性部件构成的动力学系统的建模、计算和控制.多柔体系统不仅具有柔体大变形导致的几何非线性,更具有大范围刚体运动引起的几何非线性,其非线性程度远高于计算结构力学所研究的几何非线性问题.本文基于李群局部标架(local frame of Lie group, LFLG),讨论如何发展一套新的多柔体系统动力学建模和计算方法体系, 具体内容包括:基于局部标架的梁、板壳单元,适用于长时间历程计算的多柔体系统碰撞动力学积分算法,结合区域分解技术的大规模多柔体系统动力学并行求解器, 以及若干验证性算例.上述基于李群局部标架的方法体系可在计算中消除刚体运动带来的几何非线性问题,使柔体系统的广义惯性力、广义弹性力及其雅可比矩阵满足刚体运动的不变性,使多柔体系统动力学与大变形结构力学相互统一,有望推动新一代多柔体系统动力学建模和计算软件的发展.
关键词： 局部标架方法;SE(3)群;多体系统动力学;几何非线性;几何精确理论

Abstract
The main content of the dynamics of flexible multibody systems focuses on the dynamic modeling, computation and control of complex systems composed of flexible components, which are subjected to the relative overall motion and connected by kinematical constraints. Compared with the computational structural mechanics, the multibody dynamics issues have high geometrically nonlinear, which is not only deduced by the large rotation caused from the large deformation of flexible components, but also is deduced by the overall rigid body motion. Under the concept of the local frame of Lie group (LFLG), the topic that how to develop a new modeling and computational method for flexible multibody dynamics is discussed. The major studies of this paper include the following aspects: the modelling methods of beam elements and plate/shell elements based on the LFLG, the long-time integration algorithm for the flexible multibody systems including collision problems, the parallel algorithm for multibody systems based on the domain decomposition method, and several numerical examples to verify the feasibility of the proposed method. The unique feature of the new method can eliminate the geometrically nonlinear of the overall rigid motion for flexible components. Therefore, the generalized inertial forces and internal forces as well as their Jacobian matrices are invariable under the arbitrary rigid body motion. The proposed method can motivate the integration of the modeling method of the flexible multibody dynamics and the computational structural dynamics with large deformation components and is expected to promote the development of the next-generation software of multibody system dynamics.
Keywords：local frame formulation;SE(3) group;multibody system dynamics;geometrically nonlinear;geometrically exact theory


PDF (19483KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
刘铖, 胡海岩. 基于李群局部标架的多柔体系统动力学建模与计算 ¹⁾. 力学学报[J], 2021, 53(1): 213-233 DOI:10.6052/0459-1879-20-292
Liu Cheng, Hu Haiyan. DYNAMIC MODELING AND COMPUTATION FOR FLEXIBLE MULTIBODY SYSTEMS BASED ON THE LOCAL FRAME OF LIE GROUP ¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2021, 53(1): 213-233 DOI:10.6052/0459-1879-20-292

引言

多体系统动力学的研究对象是由多个具有运动学约束、存在大范围相对运动的刚性或柔性部件组成的复杂系统,主要研究内容是这类系统的动力学建模、计算和控制.该学科与多个相邻学科具有交叉. 例如,计算结构力学领域的****研究柔性结构的大变形动态行为,这相当于一类特殊的多柔体系统动力学问题.

在历史上, 多体系统动力学与计算结构力学沿着各自的途径发展. 早期,多体系统动力学研究针对多刚体系统,关注如何描述刚体大范围运动以及刚体间的连接关系,即如何表征系统惯性力以及约束力等相互作用力;计算结构力学则研究无刚体运动的复杂结构,关注如何近似描述结构变形及表征结构的内力. 由此,两个学科的研究思路以及研究方法存在明显区别, 导致了一定的学科壁垒.

随着学科的拓展, 上述两个学科均涉足如何描述具有大变形的柔体动力学问题,在结构力学领域, 将这类问题称为几何非线性问题. 从严格意义上看,计算结构力学中的几何非线性是由于柔体大变形引起的转动所致.而多柔体系统动力学中的几何非线性不仅包括柔体大变形引起的转动,更包含由柔体大范围刚体运动引起的转动,其非线性程度远高于计算结构力学中的几何非线性问题. 因此,从几何非线性角度来看, 多柔体系统动力学建模和计算的难度显著高于计算结构力学.

梁、板/壳结构是多柔体系统的基本部件. 多柔体系统动力学建模的核心问题是:如何精确描述梁、板/壳结构的大范围刚体运动与弹性变形的耦合. 长期以来,为建立梁、板/壳结构发生大转动时的精确动力学模型,计算结构力学领域与多柔体动力学领域的****分别提出了多种有限元建模方法,例如几何精确法、绝对节点坐标法、共旋坐标法等. 其中,结构力学领域提出的几何精确方法(geometrically exact formulation,GEF)^[1]与多柔体动力学领域提出的绝对节点坐标法(absolute nodal coordinate formulation, ANCF)^[2]应用较为广泛. 从本质上看,这两类建模方法均属于非线性有限元方法.前者中早期几类单元属于含转动参数的有限元法, 后者属于无转动参数有限元方法.关于几何精确梁建模方法的系统性综述可参见文献[3];关于ANCF法的系统性综述可参见文献[4]. 文献[5,6,7]表明,这两类建模方法在计算精度、计算效率以及实现难易程度上有明显差异.

计算结构力学中的几何精确法源自Reissner^[8]对大变形梁的几何精确描述.随后, 经Simo等^[1]完善, 形成了 $\mathbb{R}^{3}$ $\times$ SO(3)群中的几何精确梁模型, 也被称为一维Cosserat杆模型.该方法用中线位置矢量与截面转动伪矢量作为节点参数来描述三维梁的大转动和大变形,中线位置矢量属于 $\mathbb{R}^{3}$线性空间,截面转动旋转矩阵属于SO(3)李群, 可视为9维线性空间中的三维流形. 基于截面转动,在梁中线上建立一个局部标架, 在该局部标架下描述梁单元的应变与角速度.由于该方法在局部坐标系中完成转动的更新, 从而可有效规避转动参数的奇异性问题.但由于转动参数不属于线性空间,转动参数的时空离散、线性化以及更新等算法与传统有限元存在很大差异,导致该算法比较复杂. 这在一定程度上制约了该方法的推广. 考虑到在 $\mathbb{R}^{3}$ $\times$ SO(3)群中建模的复杂性, 已有****构造了 $\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$空间的几何精确梁单元^[9-10].即认为转动参数属于 $\mathbb{R}^{3}$空间,采用线性空间运算进行转动变量的线性化及更新, 降低了几何精确梁单元的复杂性.需要指出的是, 上述建模所采用的构型空间并不影响单元收敛速度与精度,单元精度仅与应变及其一次变分的计算方法相关. 然而, $\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$空间的几何精确梁单元无法避免转动参数的奇异性问题,也难以通过局部标架方法在计算中来消除刚体运动导致的几何非线性.

最初, Simo等^[1]提出的几何精确梁模型在全局坐标系中求解位移和转动伪矢量的增量,梁的力平衡方程包含刚体运动带来的几何非线性.Cardona等^[11]将力平衡方程的线性化过程从当前构型切平面转换到前一收敛步的切平面,并在前一收敛步的构型空间的局部坐标系中求解转动伪矢量增量.该方法不仅能够得到对称的切向刚度矩阵,而且能同时降低刚体运动导致的几何非线性程度. 在数值计算中,这表现为刚度矩阵中的转动部分可近似满足刚体转动的不变性.Sonneville等^[12]进一步采用局部标架描述几何精确梁的平动,在SE(3)群内统一描述梁的平动与转动, 包括统一插值、线性化、更新等.上述研究指出,基于SE(3)群局部标架的几何精确梁单元能够在计算中消除单元刚体运动带来的几何非线性,其切线刚度矩阵满足刚体运动的不变性. 相比于 $\mathbb{R}^{3}$ $\times$ SO(3)群中的几何精确梁方法,在SE(3)群中的建模方法在局部标架中计算线速度及其增量.由于在不同时刻线速度并不属于同一坐标系,其对应的位移增量与截面转动伪矢量增量需要在SE(3)群内统一采用乘法更新.刘铖等^[13]的研究表明, 虽然在SE(3)群中统一插值能消除梁单元剪切闭锁,但对平动与转动独立插值以及采用缩减积分技术可进一步提高梁单元的收敛性,并能够兼容等几何分析^[14]思想. 此外,局部标架建模方法还可用于消除部分非完整约束. 例如, 对于经典力学中的雪橇运动,其质心速度在连体坐标系下一坐标轴方向速度投影为零. 若选择合适的局部标架,可直接消除该自由度, 无需施加非完整约束.上述性质也可用于构造可精确描述大转动大变形的五自由度板壳单元,消除了一自转(drilling)自由度. 总体来看,基于局部标架的几何精确梁研究尚处于初步阶段, 例如,不同运动参数化方法对单元计算精度与效率的影响、如何构造局部标架Kirchhoff几何精确梁单元以及如何构造高频耗散的时间积分方法等方面尚无相关研究.

几何精确板壳单元的发展略滞后于几何精确梁单元, 但发展历程类似.Simo等^[15-17]发表系列论文, 在 $\mathbb{R}^{3}$ $\times$ S$^{2}$群内基于Reissner-Mindlin中厚板理论提出了考虑剪切的五自由度几何精确板壳模型,包含3个平动自由度与2个旋转自由度.该单元采用中面位置矢量与沿厚度方向的单位矢量作为广义坐标描述板壳运动,其中位置矢量属于线性空间 $\mathbb{R}^{3}$,而单位矢量属于二维球面流形S$^{2}$. 为了定义厚度方向单位矢量唯一的旋转矩阵,将单元厚度方向的单位矢量约束在S$^{2}$中并沿测地线进行旋转,约束节点角速度矢量不含沿厚度方向分量, 得到相应的旋转伪矢量. 在局部标架中,可自然消去自转自由度. 并可保证沿厚度方向的矢量是单位矢量,不需要增加额外的约束方程. 但五自由度板壳单元难以直接处理含约束的多体系统,例如, 板壳结构与梁结构的连接问题. 为此, Simo等^[18-19]进一步通过单元中面变形梯度张量的极分解,建立了自转自由度与中面运动之间的关系,基于约束变分原理提出了含自转自由度的六自由度几何精确板壳单元. 同时,引入自转自由度后, 还可用于改进单元的收敛性,克服对非规则网格的敏感性问题^[20]. 本文研究表明,在SE(3)群中构造的六自由度局部标架几何精确板壳单元可自然地消除几何非线性.然而, 由于五自由度板壳单元局部标架无法描述中面自转运动,导致只能消除部分几何非线性. 如何进一步改进五自由度板壳单元,使其同样能够完全消除几何非线性值得进一步研究.

此外, 在几何精确板壳单元中, 可选择不同的转动参数进行离散. 例如,选择转动伪矢量增量、全局坐标系下沿厚度方向的单位矢量增量、或是局部标架下沿厚度方向的单位矢量增量等.当然, 不同的转动参数化方法导致算法复杂性及收敛性不同.Sonneville^[21]对SE(3)几何精确板单元建模方法进行了初步研究,但所构造的板单元在收敛性方面存在一定问题, 同时单元适用性方面尚不完善. 目前,关于SE3群中的几何精确板壳单元研究也处于初步阶段.

在多柔体系统动力学领域,Shabana^[2]于1996年提出了描述柔体大范围运动与大变形耦合的绝对节点坐标方法.该方法在处理柔体转动时, 摒弃复杂的转动参数,采用节点位置矢量以及沿物质标架的斜率矢量来描述柔体运动.由于上述矢量均属于线性空间, 绝对节点坐标法通俗易懂,在多柔体系统动力学领域受到广泛关注. 然而,该方法在计算效率和收敛性方面存在一定的缺陷,在非线性有限元领域得到的关注较少. 在该方法中, 根据对单元变形的假设,可将绝对节点坐标单元分为全参数(fully parameterized)单元^[22]与缩减(slope deficiency)单元^[23]两类. 例如, 全参数梁单元的每节点含12个广义坐标,通过单元中线/面斜率矢量以及截面/厚度方向斜率矢量描述单元的运动与变形,采用三阶Hermite插值对中线/面位置矢量进行插值,同时采用一阶Lagrange插值对截面/厚度斜率矢量进行插值, 能够描述梁的截面变形,可认为属于一类实体梁单元;缩减梁、板壳单元则分别采用Euler-Bernoulli梁和Kirchhoff板假设,描述细长梁与薄板壳. 此类单元忽略横向剪切变形,采用中线/面位置矢量与斜率矢量描述单元运动.与上述思路类似的刚体动力学建模方法为自然坐标方法^[24],它直接采用正交矩阵元素与刚体任意点的位置矢量作为参数描述刚体运动.由于在同一惯性坐标系下描述多体系统运动, 此类建模方法较为简便.自然坐标方法与绝对节点坐标方法统也因此被统称为绝对坐标方法. 然而,多柔体系统中柔体大范围刚体运动所带来的几何非线性将完全保留于系统动力学方程,会给动力学求解带来不便.

基于局部标架思想, 可在绝对节点坐标方法单元的任意物质点上构造局部标架,从而用李群方法将现有全局坐标系下的单元改造为局部标架下的绝对节点坐标单元.事实上,基于局部标架思想可以改造几乎所有考虑几何非线性的梁、板壳以及实体单元,使其在计算中消除刚体运动导致的几何非线性, 同时继承该单元的原有特征,例如单元收敛性、闭锁处理方法等.

本文以梁、板壳结构为例, 阐述如何发展一套基于李群局部标架(local frame of lie group)的多柔体系统动力学建模和计算方法体系.该体系不同于以往在惯性坐标系下(如绝对坐标系)或在柔体的某个特殊连体坐标系下(如浮动坐标系)描述柔体运动,而是在柔体任意点的李群局部标架中表征该点的运动和变形,通过李群运算完成柔体物理量在局部标架与全局标架之间的转换,进而建立和求解多柔体系统动力学方程.

上述方法体系的特点是, 对于柔体的大范围运动和大变形耦合问题,可在计算中消除物体整体运动中所包含的刚体运动,从根本上规避刚体运动导致的几何非线性. 该方法同样适用于大变形结构力学问题,可消除由大变形所带来的大转动刚体运动. 消除刚体运动后,多柔体系统的有限元模型与大变形结构力学有限元模型仅含以弯曲应变为主导的几何非线性,具有类似的数值特性. 例如,单元广义惯性力、弹性力及其雅可比矩阵均满足对任意刚体运动的不变性,单元几何刚度矩阵的贡献可弱化. 因此,具有大范围刚体运动与大变形耦合特征的多柔体系统动力学与仅考虑大变形的计算结构力学可趋于统一.这有望推动多柔体系统动力学与计算结构力学的相互融合,在此基础上形成新一代的多柔体系统动力学建模和计算软件.

1 基于SE(3)群局部标架的几类梁单元

本节以含转动参数的几何精确梁与几种无转动参数的梁单元为例,阐述基于SE(3)群局部标架的梁单元建模方法.

1.1 几何精确Timoshenko梁单元

现以圆截面梁单元为例, 给出基于SE(3)群局部标架的几何精确梁构造方式.该方法可方便地推广至构造任意截面曲梁单元或变截面梁单元. 在几何精确梁模型中,采用Timoshenko梁模型假设, 即梁具有刚性横截面,梁的构型由梁的中线位置与刚性横截面转动唯一确定.

首先, 在三维欧氏空间 $\mathbb{R}^{3}$中建立正交的惯性坐标系($O$-$e_{1}$-$e_{2}$-$e_{3})$,描述图1所示空间梁由初始构型到当前构型的运动, 包括梁的刚体运动和变形之耦合.

图1

新窗口打开|下载原图ZIP|生成PPT
图1三维几何精确梁的初始构型与当前构型

Fig.1Initial and current configurations of a geometrically exact beam element of three dimensions



不失一般性, 设梁在初始构型时的中线沿着惯性坐标系的$X$轴方向, 采用 $\mathbb{R}^{3}$中的矢量 $\varphi^{0}$描述梁的中线位置. 以矢量$\varphi^{0}$的端点为原点, 建立梁的局部连体坐标系, 其方向与惯性坐标系一致, 记其正交基矢量为($i^{0}_{1}$,$i^{0}_{2}$, $i^{0}_{3})$. 此时, 局部连体坐标系与惯性坐标系之间的旋转矩阵可记为$R^{0}=I_{3}$.

当梁抵达当前构型时, 描述梁中线的位置矢量 $\varphi ^{0}$变为矢量 $\varphi $, 上述局部连体坐标系随着梁的刚性横截面转动, 其基矢量变为正交基矢量($i_{1}$-$i_{2}$-$i_{3})$,

局部连体坐标系与惯性坐标系间的旋转矩阵为$R$. 本文简称上述局部连体坐标系为局部标架.

基于上述局部标架和SE(3)群理论, 梁单元的中线位置矢量及梁的横截面旋转矩阵合并为运动张量$H$, 可表示为

(1)$\begin{eqnarray} \label{eq1} {H}\left( {\xi_{1} } \right)=\left[ {{\begin{array}{c@{\quad }c} {{R}\left( {\xi_{1} } \right)} & {{\varphi }\left( {\xi_{1} } \right)} \\ 0 & {1} \\ \end{array} }} \right]=\left[ {{\begin{array}{c@{\quad }c} {{R}^0} & {{\varphi }^{0}} \\ 0 & {1} \\ \end{array} }} \right]\left[ {{\begin{array}{c@{\quad }c} {\Delta {R}} & {{u}} \\ 0 & {1} \\ \end{array} }} \right] \end{eqnarray}$
其中, $\xi_{1}\in [0, L_{0}]$表示梁中线上$P$点在初始构型中的弧长坐标, $L_{0}$为梁单元长度, $u$为该点在初始构型局部标架中的位移矢量, $\Delta R$为初始构型与当前构型之间的旋转变换矩阵. $H$矩阵称为欧几里得变换, 它包含了梁截面作刚体平动和转动的完整信息.需要指出的是, SE(3)群中的欧几里得变换矩阵除了具有式(1)这样的形式,还可表示为六阶方阵形式^[25].

对于式(1)中在局部标架描述的运动增量, 可通过SE(3)群的指数映射,转换为初始构型下的位移矢量与旋转矩阵, 其表达式为

(2)$\begin{eqnarray} \label{eq2} \exp_{SE(3)} \left( {{a}} \right)=\left[ {{\begin{array}{c@{\quad }c} {\Delta {R}} & {{u}} \\ {\bf{0}_{3\times 1} } & {1} \\ \end{array} }} \right],\\ {a}\triangleq\left[ {{\begin{array}{c@{\quad }c} {{\tilde{{\varTheta }}}} & {{\bar{{u}}}} \\ {\bf{0}_{3\times 1} } & 0 \\ \end{array} }} \right] \end{eqnarray}$
其中, ${\varTheta }$为描述刚体转动的伪矢量, 符号$\tilde{{\cdot}}$帽子映射代表对应$\cdot $的反对称矩阵, ${\bar{{u}}}$为在局部标架中度量的位移矢量. 在SE(3)群中, 上述指数映射可表示为

(3)$\begin{eqnarray} \label{eq3} \exp_{SE(3)} ({a})=\left[ {{\begin{array}{c@{\quad }c} {\exp_{SO(3)} ({\tilde{{\varTheta }}})} & {{ T}_{SO(3)}^{T} ({\varTheta }){\bar{{u}}}} \\ {\bf{0}_{1\times 3} } & 1 \\ \end{array} }} \right] \end{eqnarray}$
其中, $\exp_{SO(3)} ({\tilde{{\varTheta }}})$为SO(3)群的指数映射, ${T}_{SO(3)} ({\varTheta })$为SO(3)群的切线算子.

关于SO(3)群与SE(3)群基础理论可参考文献[21].

进一步, 在初始构型与当前构型中, 梁单元上任意点的位置矢量可分别表示为

(4)$\begin{eqnarray} \left.\begin{array}{l} {{X}={\varphi }^{0}\left( {\xi_{1} } \right)+{R}^{0}\left( {\xi_{1} } \right)\left( {\xi_{2} {i}_{2}^{0} +\xi_{3} { i}_{3}^{0} } \right)} \\ {{x}={\varphi }\left( {\xi_{1} } \right)+{R}\left( {\xi_{1} } \right)\left( {\xi_{2} {i}_{2}^{0} +\xi_{3} {i}_{3}^{0} } \right)} \\ \end{array}\right\} \end{eqnarray}$
其中, $\xi_{2}$与$\xi_{3}$为梁截面上P点在局部标架中的坐标分量.

对于动力学问题, 梁单元在任意时刻的中线速度和截面转动角速度, 可通过运动学方程, 也被称为Poisson方程, 表示为

(5)$\begin{eqnarray} \label{eq5} \begin{array}{*{20}c} \left[\begin{array}{c@{\\ }c} {{\dot{{R}}}\left( {\xi_{1} } \right)} & {{\dot{{\varphi }}}\left( {\xi_{1} } \right)} \\ 0 & 0 \\ \end{array} \right] \\ \end{array} =\left[\begin{array}{c@{\\ }c} {{R}\left( {\xi_{1} } \right)} & {{\varphi }\left( {\xi_{1} } \right)} \\ 0 & {1} \\ \end{array} \right]\left[ \begin{array}{c@{\\ }c} {{\tilde{{\omega }}}\left( {\xi_{1} } \right)} & {{U}\left( {\xi_{1} } \right)} \\ 0 & 0 \\ \end{array} \right] \end{eqnarray}$
其中, 顶标表示矢量函数对时间$t$的导数, $U$与 $\omega$为梁单元中线和横截面的线速度与角速度在局部标架下的投影. 此时,惯性力所做的虚功可表示为

(6)$\begin{eqnarray} \delta W_{{iner}} =-\int {\delta{h}^{T}\left( {{C\dot{{v}}}-{\hat{{v}}}^{T}{Cv}} \right){d}\xi_{1} } ,\\ {C}=\left[ {{\begin{array}{c@{\\ }c} {\rho_{A} {I}_{3} } & {\bf{0}_{3} } \\ {\bf{0}_{3} } & {{J}} \\ \end{array} }} \right] \end{eqnarray}$
其中, ${v}=\left[ {{\begin{array}{c@{\quad }c} {{U}} \\ {{\omega }} \\ \end{array} }} \right]$, ${\hat{{v}}}\left[ {{\begin{array}{*{20}c} {{\tilde{{\omega }}}} & {{\tilde{{U}}}} \\ {\bf{0}_{3\times 3} } & {{\tilde{{\omega }}}} \\ \end{array} }} \right]$, $\delta{h}=\left[ {{\begin{array}{*{20}c} {\delta{d}} \\ {\delta{\eta }} \\ \end{array} }} \right]$, $\delta d$和$\delta $ $\eta$是梁单元中线和横截面在局部标架下的虚位移及虚转角, $\rho

_{A}$表示单位长度梁的材料密度, $J$为单位长度梁的惯性矩阵.上述转动项与传统方法一致, 但由于采用局部标架描述平动,平动惯性力的形式与转动惯性力的形式一致, 存在转动与平动的耦合项,与绝对导数与相对导数运算相关.

根据连续介质力学, 梁单元内任意点的变形梯度张量可表示为

(7)$\begin{eqnarray} F=\frac{{d}{x}}{{d}{X}}=\left[\begin{array}{c@{\quad }c@{\quad }c} \varphi '+R'\left(\xi_{2}{i}_{1}^{0} +\xi_{3} {i}_{2}^{0} \right) & {Ri}_{2}^{0} & Ri_{3}^{0} \\ \end{array} \right] \end{eqnarray}$
其中, $()'$表示向量函数对$\xi_{1}$的偏导数.对于小应变问题, 可根据Timoshenko梁的变形假设,将梁单元内的应变表示为如下矢量形式

(8)$\begin{eqnarray} {\varepsilon }=\left[ {{\begin{array}{*{20}c} {{\varGamma }} \\ {{\kappa }} \\ \end{array} }} \right]={f}-{f}_{0} ,\quad {f}=\left[ {{\begin{array}{*{20}c} {{R}^{T}\varphi'} \\ {\mbox{axial}({R}^{T}R')} \\ \end{array} }} \right] \end{eqnarray}$
其中, axial()表示反对称矩阵对应的轴矢量, ${\varGamma }$表示梁中线的拉伸应变和横向剪切应变, $\kappa $的分量$\kappa_{1}$ 表示梁的横截面扭转应变, 分量$\kappa_{2}$ 与$\kappa_{3}$分别表示横截面绕$i_{2}$与$i_{3}$轴的弯曲应变.对于线弹性材料, 通过平面应力假设以及截面积分,可将梁单元内力所做的虚功表示为

(9)$\begin{eqnarray} \delta W_{int} =\int {\delta{\varepsilon }\cdot { S}{d}\xi_{1} } ,\\ {S}=\left[ {{\begin{array}{*{20}c} {{N}} \\ {{M}} \\ \end{array} }} \right]=\left[ {{\begin{array}{c@{\quad }c} {{D}_{1} } & {\bf{0}_{3\times 3} } \\ {\bf{0}_{3\times 3} } & {{D}_{2} } \\ \end{array} }} \right]{\varepsilon } \end{eqnarray}$
其中, $N=D_{1}\varGamma$表示局部标架下的梁截面内力, $M=D_{2}\kappa $表示局部标架下的梁截面内力矩. $D_{1}$与$D_{2}$为弹性矩阵, 可分别表示为

(10)$\begin{eqnarray} {D}_{1} =\left[ {{\begin{array}{*{20}c} {EA} & 0 & 0 \\ 0 & {k_{s1} GA} & 0 \\ 0 & 0 & {k_{s2} GA} \\ \end{array} }} \right],\quad {D}_{2} =\left[ {{\begin{array}{*{20}c} {GJ_{s} } & 0 & 0 \\ 0 & {EI_{1} } & 0 \\ 0 & 0 & {EI_{2} } \\ \end{array} }} \right] \end{eqnarray}$
其中, $E$和$G$分别为材料的拉伸模量和剪切模量,$A$, $I_{1}$, $I_{2}$分别表示梁截面面积和截面惯性矩, $J_{s}$表示圣维南扭转常数, $k_{s1}$和$k_{s2}$为相应的剪切修正系数. 相应的, 外力所做的虚功表示为

(11)$\begin{eqnarray} \hspace{-3mm}\delta W_{{ext}} =\int \left( {\delta{ d}^{T}{P}+\delta{\eta }^{T}{T}} \right){d}\xi_{1} +\delta{d}(\xi )^{T}{\bar{{P}}}+\delta{\eta }(\xi )^{T}{\bar{{T}}}\quad \end{eqnarray}$
其中, $P$表示在局部标架下梁中线上的分布力, $T$表示在局部标架下梁截面上的分布力矩, ${\bar{{P}}}$表示加载于梁中线的集中力, ${\bar{{T}}}$表示加载于梁截面的集中力矩.

由式(6)、式(9)以及式(11), 可构造无约束梁单元的动力学方程. 在此基础上,通过对梁单元的动力学方程进行空间域有限元离散以及时间域差分离散,可得到一个非线性代数方程组, 进而进行数值求解.

在对上述动力学方程进行空间离散时, 如何保证插值算法的客观性至关重要. 例如,Jeleni$\acute{c}$与Crisfiled^[26]指出, 若直接用有限元方法对转动伪矢量插值,将不满足客观性条件; 即单元刚体转动将产生虚假应变能,应变更新算法则会对超弹性材料带来路径相关性问题. 随后,他们借鉴共旋坐标系概念, 提出一种相对插值方法,解决了几何精确建模方法插值的客观性问题. Ibrahimbegovic 等^[27]指出,若采用更新算法计算单元应变, 直接对转动增量进行插值,同样能够满足离散应变的客观性. 近年来,Bauchau等^[28]系统地讨论了几类常用转动/运动插值算法是否满足客观性、矢量性、本征性以及它们的计算效率.其中, 插值算法的本征性概念是指插值算法是否依赖某类特殊的转动参数化方法.

此外, 在构造离散系统动力学方程时, 需要对应变进行一次变分. 在线性空间内,有限元离散与变分操作可交换顺序. 但在SE3群中, 两者交换顺序并不等价.对于先离散后线性化, 一般称为离散一致线性化; 而对于先线性化后离散,则称为连续一致线性化. 为简单起见, 下面给出连续一致线性化过程.通过SE(3)群内的变分运算规则, 将几何精确梁单元的应变矢量在局部标架下进行变分,可表示为

(13)$\begin{eqnarray} \delta W_{int} =\delta{q}^{T}\int {{\varPsi }^{T}{B}^{T}{S}{d}\xi_{1} } ,\quad \delta{q}=\left[ {{\begin{array}{*{20}c} {\delta{\varphi }^{Node}} \\ {\delta{\varTheta }^{Node}} \\ \end{array} }} \right] \end{eqnarray}$
其中, ${\varPsi }$为有限元离散形函数矩阵. 本文选用一阶Lagrange插值进行离散,并采用选择性缩减积分处理单元的剪切闭锁问题. 为了得到切线刚度矩阵,再进一步作二次变分(即线性化)可得

(14)$\begin{eqnarray} \label{eq14} \Delta \left( {\delta W_{int} } \right)=\delta{ q}^{T}\int {{\varPsi }^{T}\left( {{B}^{T}{ DB}+{\breve{{S}}B}} \right){\varPsi }{d}\xi_{1} } \Delta {q} \end{eqnarray}$
其中, 第一项包含离散形式的材料刚度矩阵、第二项包含离散形式的几何刚度矩阵,${\breve{{S}}}\triangleq\left[ {{\begin{array}{c@{\quad }c} {{\bf 0}_{3\times 3} } & {{\tilde{{N}}}} \\ {{\tilde{{N}}}} & {{\tilde{{M}}}} \\ \end{array} }} \right]$.

从材料刚度矩阵与几何刚度矩阵的表达式可见, 其仅与局部标架下的应变相关,而与刚体运动无关. 因此, 单元刚度矩阵自然满足刚体运动的不变性,也就是消除了刚体运动导致的几何非线性. 值得指出,由于SE(3)群中的运算不具有交换性, 此处的单元切线刚度矩阵是非对称的. 当然,也可通过转换至前一收敛步切平面的方式构造对称形式的刚度矩阵,其具体推导可参见文献[11]. 由于上述结果消除了刚体运动导致的几何非线性,而在多数动力学计算问题中, 几何刚度矩阵可忽略不计. 此时,可获得一个近似的对称刚度矩阵. 值得指出, 相比于 $\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$线性空间的几何精确梁单元,基于SE(3)群局部标架的空间梁单元弹性力及其刚度矩阵表达式非常简洁,可显著提高计算效率.

消除刚体运动产生的几何非线性后, 梁的几何非线性主要由其弯曲应变引起.而随着单元数目的增加, 积分形式的弯曲应变逐渐减少, 应变一次变分$B$矩阵趋于常数矩阵, 材料刚度矩阵也将趋于一常数矩阵,梁的几何非线性程度进一步减弱. 因此,对于足够稠密的有限元网格,系统刚度矩阵可近似为一常数矩阵, 与线性有限元具有类似的数值特性. 然而,在工程实际问题中, 在保证空间离散收敛前提下,为平衡刚度矩阵的更新次数与线性方程组的求解规模, 应选择合适的有限元网格尺寸.上述性质对于任意的局部标架几何非线性有限单元都是适用的.

在时间离散方面,几何精确方法最为常用的是李群a类算法与李群保能量-(角)动量算法. 为简单起见,此处给出对应SE(3)群描述的广义a方法的基本离散格式.李群广义a方法对前后两个相邻时刻运动增量进行差分离散,并通过SE(3)指数映射递推运动张量$H$, 其具体离散格式为

(15)$\begin{eqnarray} \label{eq15} \left.\begin{array}{l} {d}_{n} =\Delta t{v}_{n} +\Delta t^{2}(0.5-\beta ){a}_{n} +\Delta t^{2}\beta {a}_{n+1} \\ {v}_{n+1} ={v}_{n} +\Delta t(1-\gamma ){a}_{n} +\Delta t\gamma {a}_{n+1}\\ \left( {1-\alpha_{m} } \right){a}_{n+1} +\alpha_{m} {a}_{n} =\left( {1-\alpha_{f} } \right){\dot{{v}}}_{n+1} +\alpha_{f} {\dot{{v}}}_{n}\\ {H}_{n+1} ={H}_{n} \exp_{SE3} \left( {{\tilde{{d}}}_{n} } \right)\\ \end{array}\right\} \end{eqnarray}$
其中, $d_{n}$表示系统从$n$时刻到$n+$1时刻在局部标架下的运动增量, $\alpha_{m}$, $\alpha_{f}$, $\beta $, $\gamma $为广义a方法算法的参数, $\Delta t$为离散的时间步长, ${H}_{n} =\left[ {{\begin{array}{*{20}c} {{R}_{n} } & {{\varphi }_{n} } \\ {{0}_{1\times 3} } & 1 \\ \end{array} }} \right]$. 由于通过前后两个时间步增量进行系统构型更新,运动增量$d_{n}$的范数通常为小量, 不会出现转动参数的奇异性问题,也可避免转动参数在周期临界位置的特殊处理. 系统广义质量矩阵的离散形式可表示

(16)$\begin{eqnarray} && \Delta \left( {\delta W_{{iner}} } \right)= -\rho_{A} \delta{q}^{T}\left[\int_l {\varPsi }^{T}\left( {M}_{{m}} \frac{{d}^{2}}{{d}t^{2}}+\right.\right. \left.\left.{M}_{{gyr}} \frac{{d}}{{d}t}+{M}_{{cent}} \right){\varPsi }{d}\xi_{1} \right]\delta{q} \end{eqnarray}$
其中, 广义质量矩阵的质量项$M_{m}$、陀螺力项$M_{gyr}$和离心力项$M_{cent}$分别为

(17)$\begin{eqnarray} \label{eq17} \left.\begin{array}{l} M_{m} ={C}\\ M_{{gyr}} =C\hat{{v}}-\hat{{v}}^{T}{C}-\breve{(Cv)}\\ M_{{cent}} =C\dot{{\hat{{v}}}}-\hat{{v}}^{T}{C\hat{{v}}}-\breve{(Cv)}\hat{v}\\ \end{array} \right\} \end{eqnarray}$
其中, $v=\left[\begin{array}{c} U\\ \omega\\ \end{array}\right]$, $\breve{{v}}\triangleq\left[\begin{array}{cc} {\bf 0}_{3\times 3} & \tilde{{U}} \\ \tilde{U} & \tilde{\omega }\\ \end{array}\right]$. 由上式可知, 在动力学计算过程中,时间离散后广义质量矩阵中包含$O(h^{-2})M_{m}+O(h^{-1})M_{m}+M_{cent}$,其中主要贡献项$O(h^{-2})M_{m}$为常数矩阵,陀螺项与离心力项相对贡献较小. 对于大多工程问题,陀螺力项与离心力项可几乎保持不变. 若对于极端绕三轴高转速旋转动力学问题,需保留陀螺项与离心力项. 此外,广义质量矩阵表达式仅与系统的惯性特性、速度以及加速度相关, 与转动参数无关.这表明, 数值计算过程中不会出现奇异性问题.

值得指出,上述基于SE(3)群局部标架的几何精确梁建模方法可退化为基于SE(3)群局部标架的刚体动力学建模方法,且具有避免转动参数奇异性以及降低刚体转动非线性的优点.

1.2 无转动参数的梁单元

对于无转动参数的梁单元, 也可建立基于SE(3)群局部标架,进而消除梁单元刚体运动导致的几何非线性.

首先, 考察基于ANCF描述的全参数梁单元.采用中线位置矢量以及沿中线方向的斜率矢量作为广义坐标,则梁单元中线上任意一点的位置矢量$r^{mid}$可通过三阶Hermite插值函数表示为

(18)$\begin{eqnarray} \label{eq18} {r}^{{mid}}(X)={S}^{{mid}}(X)\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c} {{r}_{i}^{T} } & {{r}_{i,X}^{T} } & {{r}_{j}^{T} } & {{r}_{j,X}^{T} } \\ \end{array} }} \right]^{T} \end{eqnarray}$
上式中的形函数矩阵可表示为

(19)$\begin{eqnarray} \label{eq19} {S}^{{mid}}(X)=\left[ {{\begin{array}{c@{\quad }c@{\quad }c@{\quad }c} {S_{1} {I}_{3} } & {S_{2} {I}_{3} } & {S_{3} {I}_{3} } & {S_{4} {I}_{3} } \\ \end{array} }} \right] \end{eqnarray}$
其中, $S_{1}=1-3\xi^{2}+2\xi^{3}$, $S_{2}= L_{0}(\xi -2\xi^{2}+\xi^{3})$, $S_{3}=3\xi^{2}-$ $2\xi^{3}$, $S_{4}=L_{0}(-\xi^{2}+\xi^{3})$, $\xi =X/L_{0}$为单元局部坐标, $L_{0}$为梁单元的初始长度, $I_{3}$为三阶单位矩阵. 此外, 通过一阶线性Lagrange插值得到梁截面的斜率矢量$r$, $Y$与$r$, $Z$. 因此, 该梁单元上任意点的位置矢量可表示为

(20)$\begin{eqnarray} \label{eq20} {x}({X})={r}^{{mid}}(X)+Y{r}_{i,Y} (\xi )+Z{ r}_{i,Z} (\xi )={Se} \end{eqnarray}$
其中, $S=[S_{1}I_{3}\\ S_{2}I_{3}\\ S_{5}I_{3}\\ S_{7}{I}_{3}\\ S_{3}{I}_{3}\\ S_{4}{I}_{3}\\ S_{6} {I}_{3}\\ S_{8}{I}_{3}]^{T}$, $S_{5}= Y(1-\xi)$, $S_{6}=Y\xi $,$S_{7}=Z(1-\xi)$, $S_{8}= Z\xi$, $X=[X \\ Y \\ Z]^{T}$.

该单元在轴向与横向的插值阶数不同, 会产生闭锁问题, 处理方法可详见文献[29].

该单元上任意点的局部标架可通过该点的3个斜率矢量确定.虽然局部标架的构造方法不唯一, 但这些方法的数值特性类似. 例如,可采用如下局部标架

(21)$\begin{eqnarray} \label{eq21} \left.\begin{array}{l} {R}=\left[\begin{array}{c@{\quad }c@{\quad }c} i_{1} & i_{2} &i_{3}\\ \end{array}\right],\\ {i}_{1} =\dfrac{{r}_{,Y} \times {r}_{,Z} }{\left\| {{r}_{,Y} \times {r}_{,Z} }\right\|} \\ {i}_{2} =\dfrac{{r}_{,Y} }{\left\| {{r}_{,Y} } \right\|},\\ {i}_{3} ={i}_{1} \times {i}_{2} \end{array}\right\} \end{eqnarray}$
由于该局部标架完全由梁的平动位移场确定,因此局部标架所对应的旋转矢量变分可表示为

(22)$\begin{eqnarray} \label{eq22} \delta{\tilde{{\varTheta }}}={R}^{T}\delta{R} \end{eqnarray}$
类似地, 可计算出局部标架的角速度以及角加速度,进而可采用基于SE(3)群描述的广义a方法求解系统动力学方程. 考虑到篇幅限制,此处对该单元的广义质量矩阵、切向刚度矩阵等不加以推导,部分细节及切向刚度矩阵高效计算方法可见文献[30].

除上述全参数ANCF梁单元外, 其余几类梁单元均采用Euler-Bernoulli梁的变形假设.若采用位移场的二阶导数计算连续形式的弯曲应变,则单元构造需要采纳至少$C^{1}$连续的插值函数. 对于ANCF缩减梁单元,为满足上述要求, 其节点坐标包含节点位置以及节点沿轴向的偏导数矢量,采用三阶Hermite插值离散梁中线位移场, 如式(18)所示.任意点的局部标架可通过该点在参考构型中的斜率矢量以及在当前构型中的斜率矢量定义为

(23)$\begin{eqnarray} \label{eq23} {R}=\left( {{E}\cdot {t}} \right){I}_{3} +\widetilde{{ E}\times {t}}+\frac{\left( {{E}\times {t}} \right)\otimes \left( {{E}\times {t}} \right)}{1+{E}\cdot {t}} \end{eqnarray}$
其中, ${E}={{r}_{,X}^{0} }\Big/{\left\| {{r}_{,X}^{0} } \right\|}$, ${t}={{r}_{,X} }\Big/{\left\| {{r}_{,X} } \right\|}$, $\otimes $为并矢符号. 关于ANCF缩减梁单元构造可见文献[31].

若进一步需要考虑梁扭转变形, 或矩形截面梁的弯曲变形, 细节推导可见文献[32], 其中局部标架则需在式(23)基础上考虑扭转运动.

此外, 为减弱对梁中线位移场高阶连续性要求, 减少单元广义坐标个数,有****提出了离散弹性杆(discrete elastic rods, DER)模型^[33].该模型通过离散方式描述杆的轴向、弯曲以及扭转应变, 具有计算效率高的优点,在计算几何领域被广泛应用. 值得指出,该模型的弯曲变形描述与Euler-Bernoulli梁一致, 只是中文直译为离散弹性杆单元.当不考虑扭转变形且网格均匀划分时, 杆单元节点处的局部标架也表示为式(23)构造,式中的单位矢量E与t定义如下

(24)$\begin{eqnarray} {E}=\frac{{\bar{{E}}}_{m} +{\bar{{E}}}_{m} }{\left\| {{\bar{{E}}}_{m} +{\bar{{E}}}_{m} } \right\|},\\ {t}=\frac{{\bar{{t}}}_{m} +{\bar{{t}}}_{m+1} }{\left\| {{\bar{{t}}}_{m} +{\bar{{t}}}_{m+1} } \right\|} \end{eqnarray}$

其中, ${\bar{{E}}}_{m} =\dfrac{{r}_{m+1}^{0} -{r}_{m}^{0} }{\left\| {{r}_{m+1}^{0} -{r}_{m}^{0} } \right\|},\\ {\bar{{t}}}_{m} =\dfrac{{r}_{m+1} -{r}_{m} }{\left\| {{r}_{m+1} -{r}_{m} } \right\|}$.

近年来, 随着CAE与CAD的融合, 等几何分析^[14]受到广泛关注.等几何分析将构造CAD几何造型的样条函数(如B$\acute{e}$zier、B-spline、非均匀有理B样条NURBS和T-spline等)及其控制点与权重直接作为有限元建模的输入信息.考虑到三阶NURBS样条在特殊情况下可退化为三阶Hermite插值函数,针对一类特殊等几何梁单元(三阶NURBS插值, 除首尾节点重复度为四外,其余节点矢量重复度为二),可类似对上述ANCF缩减梁单元的讨论来构造形如式(23)的局部标架,对于由一般NURBS样条函数离散的梁单元, 如何构造局部标架是个难题,尚值得进一步研究.

2 基于SE(3)群局部标架的几类板壳单元

本节基于SE(3)群局部标架, 介绍含转动参数的几何精确板壳单元以及几类无转动参数的板壳单元.

2.1 三维几何精确Reissner-Mindlin板壳单元

Simo等^[15]首先提出了$\mathbb{R}^{3}$ $\times$ S$^{2}$几何精确板壳单元, 随后,并将该单元推广至可考虑变厚度板壳单元^[16]、以及考虑弹塑性本构的板壳单元^[17].本文进一步将上述单元推广至SE(3)群局部标架中, 消除刚体运动导致的几何非线性.

图2

新窗口打开|下载原图ZIP|生成PPT
图2几何精确板壳单元中面的初始和当前构型

Fig.2Initial and current configurations of a geometrically exact plate/shell element of three dimensions



在几何精确板壳理论中, 假设板的构型为不可伸展的单方向Cosserat曲面,即板的厚度不发生变化, 沿厚度方向的单位矢量属于S$^{2}$空间.在初始构型与当前构型中, 板单元上任意点的位置矢量可分别表示为

(25)$\begin{eqnarray} \label{eq25} \left.\begin{array}{*{20}c} {{X}={\varphi }^{0}\left( {\xi_{1} ,\\ \xi_{2} } \right)+{R}^{0}\left( {\xi_{1} ,\\ \xi_{2} } \right)\xi_{3} {E}} \\ {{x}={\varphi }\\ \left( {\xi_{1} ,\\ \xi_{2} } \right)+{R}\\ \left( {\xi_{1} ,\\ \xi_{2} } \right)\xi _{3} {E}} \\ \end{array}\right\} \end{eqnarray}$
其中, 局部标架对应的旋转矩阵$R$可通过厚度方向的单位矢量确定,如式(23). $E=t_{0}$为初始构型中沿厚度方向的单位矢量, $t$为当前构型中沿厚度方向的单位矢量. 由式(23)可见,局部标架在时间域的演化并非自由运动, 其转动不含沿厚度方向的分量. 换言之,沿着该方向的转动增量与角速度分量均为零. 在局部标架中,可根据该约束直接消除一个自转自由度, 即单元的节点自由度为5.

根据连续介质力学, 板壳中任意点的应变可分解为薄膜应变 $\varepsilon$、弯曲应变 $\kappa $以及横向剪切应变 $\gamma $,并可在局部标架中分别表示为

(26)$\begin{eqnarray} \label{eq26} \left.\begin{array}{l} \varepsilon_{\alpha \beta } =\dfrac{1}{2}\left( {\bar{{\varphi}}_{,\alpha } \cdot \bar{{\varphi }}_{,\beta } -\bar{{\varphi }}_{,\alpha}^0 \cdot \bar{{\varphi }}_{,\beta }^0 } \right)\\ \kappa_{\alpha \beta } =-\bar{{\varphi }}_{,\alpha } \cdot \left( {{\tilde{{E}}\kappa }_{\beta } } \right)+\bar{{\varphi }}_{,\alpha }^0 \cdot \left( {{\tilde{{E}}\kappa }_{\beta }^{0} } \right) \\ \gamma_{\alpha } =\left( {\bar{{\varphi }}_{,\alpha } -\bar{{\varphi }}_{,\alpha }^0 } \right)\cdot {E}\\ \end{array} \right\} \end{eqnarray}$
其中, $\alpha =1, 2$, $\beta =1, 2$, $\bar{{\varphi }}_{,\alpha } ={R}^{T}\bar{{\varphi }}_{,\alpha } $, ${\tilde{{\kappa }}}_{\alpha} ={R}^{T}{R}_{,\alpha } $. 若薄膜应变为小量,则弯曲应变可简化为

(27)$\begin{eqnarray} \label{eq27} \kappa_{\alpha \beta } =-\bar{{\varphi }}_{,\alpha }^0 \left[ {{\tilde{{E}}}\left( {{\kappa }_{\beta } -{\kappa }_{\beta }^{0} } \right)} \right] \end{eqnarray}$
此时, 式(26)退化为文献[21]中的应变表达式, 弯曲应变与薄膜应变解耦,适用于描述以弯曲变形为主导的板壳结构. 若采用式(27)的弯曲应变,则在板壳单元纯弯曲测试中单个单元即可达到解析解精度. 进一步,通过与1.2节类似的时空离散方法、一次变分以及两次变分等操作,可构造出基于SE(3)群局部标架的五自由度几何精确板壳单元.

若进一步通过中面变形梯度张量的极分解,建立自转自由度与中面运动之间的关系,可推导出SE(3)群局部标架的六自由度几何精确板壳单元. 此时,局部标架坐标旋转矩阵$Q$可表示为

(28)$\begin{eqnarray} \label{eq28} {Q}={R}\mbox{exp(}\theta {\tilde{{E}}}\mbox{)} \end{eqnarray}$
其中, $\theta $表示为自转角度. 对于上述六自由度板壳单元,由于所建立的局部标架包含该点的完整的转动信息,该单元能够完全消除刚体运动带来的几何非线性. 而对于五自由度板壳单元,若直接采用式(25)中旋转矩阵$R$所定义局部标架则仅能够消除部分的几何非线性. 然而,五自由度单元的缺陷也可通过构造近似的自转转动进行弥补. 考虑到篇幅限制,此处不对板壳单元给出具体推导, 将在作者后期论文中进一步展示.

对于低阶板壳单元, 单元存在严重的薄膜与剪切闭锁.本文分别采用Hellinger-Reissner变分原理和Hu-Washizu变分原理处理薄膜与剪切闭锁. 相比于文献^[21]研究工作,本文所构造的板壳单元具有更好的数值稳定性,位移场时空离散方法也可进一步提高单元收敛性.

2.2 几类无转动参数的板壳单元

现基于SE(3)群局部标架, 简要介绍几类无转动参数的板单元建模方法,包括ANCF全参数板壳单元、ANCF缩减薄板壳单元、以及等几何薄板单元.

基于局部标架的ANCF全参数板单元构造方法与ANCF全参数梁单元几乎一样,此处不再赘述. 关于全参数ANCF板单元建模与高效计算方法可见文献[34].对于ANCF缩减薄板壳单元, 作者曾基于Kirchhoff假设,引入一个局部笛卡尔标架来描述板壳单元的应变^[35]. 事实上,该局部笛卡尔标架可直接作为局部标架, 进而描述板壳单元的速度以及角速度.该局部标架的具体表示为

(29)$\begin{eqnarray} \label{eq29} \left.\begin{array}{l} R=\left[\begin{array}{c@{\quad }c@{\quad }c} i_{1} & i_{2} & i_{3}\\ \end{array}\right],\\ {i}_{1} =\dfrac{{ r}_{,X} }{\left\| {{r}_{,X} } \right\|}\\ {i}_{2} ={i}_{3} \times {i}_{1},\\ {i}_{3} =\dfrac{{r}_{,X} \times {r}_{,Y} }{\left\| {{ r}_{,X} \times {r}_{,Y} } \right\|} \end{array}\right\} \end{eqnarray}$
上述方法可推广至等几何薄板壳单元中. 与处理等几何细长梁类似,现考察一类特殊的三阶双变量NURBS离散方法, 其中首尾节点四次重复,其余节点二次重复. 该单元可等价于48自由度的ANCF缩减板壳单元,其中有限元节点自由度包含2个一次偏导数以及1个二次混合偏导数.由此可通过式(29)构造一个局部标架. 可证明,在该局部标架下建立的动力学方程能够规避刚体运动导致的几何非线性.

事实上, 基于SE(3)群局部标架的描述适用于所有描述几何非线性问题的有限单元.只要将原有有限单元拓展到SE(3)群中,即可实现多柔体系统动力学与现有非线性有限元法的融合.

3 基于局部标架的多柔体系统动力学核心算法

多体系统动力学是典型的交叉学科,其研究内容涉及多刚体系统动力学、结构静/动力学、分析力学、计算数学等学科.历史上, 多柔体系统动力学与结构静/动力学沿着各自的途径发展. 因此,在大型CAE软件中, 多柔体系统动力学模块与结构静力/动力学模块之间并不融合.例如, 在MSC软件体系中,解决大转动、大变形、刚柔耦合动力学问题时需采用多体动力学软件MSC.ADAMS与有限元软件MSC.NASTRAN进行联合仿真.其原因是, 在多体系统动力学软件开发之初, 仅有多刚体系统动力学作为理论基础,并未将其与有限元理论相融合.现有的其他几类多体系统动力学软件也存在类似的问题,这给多柔体系统的动力学计算带来极大不便.

基于上述SE(3)群局部标架, 融合计算几何力学与非线性有限元理论,有望发展一套新的多柔体系统动力学建模与计算方法体系; 结合等几何分析,可开发一类消除刚体运动非线性的刚体、梁、板、壳以及实体单元,实现多柔体系统动力学与结构大变形动力学建模与计算方法的统一. 此外,这样的多柔体系统动力学软件需具备用于长历程动态仿真且高频耗散可控的通用时间积分算法,能模拟多柔体系统的复杂碰撞问题,同时实现大规模多柔体系统动力学计算的通用并行异步仿真,以满足复杂系统的动力学仿真需求. 上述方法体系的发展涉及如下几个关键科学问题:SE(3)群中三维运动参数化及基本运算规则问题;不同构型空间建模方法之间的转换关系; 含碰撞的多柔体系统长时间计算精度问题;复杂多柔体系统负载平衡图分解及迭代并行算法高效率预条件算子构造问题.

具体来说, 上述多柔体系统动力学建模和计算方法可分解为如下6个方面.

(1) 基于SE(3)群局部标架的建模方法. 研究SE(3)群中运动的参数化描述方法、运动参数插值方法、动力学方程线性化运算规则;基于SE(3)群局部标架, 构造含转动参数的梁板壳单元,

例如1.1节提出的几何精确梁单元、2.1节提出的几何精确板壳单元、Kirchhoff几何精确梁/板壳单元、离散弹性杆单元;基于SE(3)群局部标架, 研究考虑截面变形的任意截面的空间梁建模方法,提出截面运动与中线运动解耦的维数降阶动力学建模方法; 基于SE(3)群局部标架,构造无转动参数的梁、板壳、实体单元, 例如Kirchhoff板壳单元, 实体单元等;分析不同运动参数化对计算精度与效率的影响;借鉴有限元方法中对闭锁问题的处理手段(如混合变分原理), 提高上述单元的收敛性.

(2) 长时间历程积分器. 保能量-(角)动量时间积分算法的构造依赖于单元运动的参数化描述方法、时空离散格式以及力平衡方程、几何方程、本构方程的形式.针对上述几类基于SE(3)群局部标架的单元,分别构造其保能量-(角)动量时间积分算法,并进一步提出高频耗散能量衰减动量守恒算法; 基于Hamilton原理,构造一类时空统一离散的变分积分子,包括运动参数李群变分积分子、Hamel变分积分子,探索变分积分子处理绳索与薄膜系统动力学响应的优势.

(3)多柔体系统的碰撞模拟算法. 借鉴有限元方法研究成果,研究SE(3)群局部标架下的碰撞问题罚函数方法与Lagrange乘子法;对于碰撞问题的隐式算法, 构建保能量-(角)动量的大步长数值计算方法;对碰撞问题的显式算法, 开发异步变分积分算法,并能够兼容显隐式混合时间积分方法.

(4) 多柔体系统的分布式通用并行异步算法.多柔体系统模型通常存在多种不同类型的单元以及运动副,提出多柔体系统加权图的表征方法、基于图论的负载平衡区域分解方法;提出基于区域分解的多柔体系统并行迭代算法, 研究界面粗网格自适应构造方法;当界面材料存在明显差异时, 研究界面问题高效的预条件处理算子构造方法,分析区域分解对界面问题条件数的影响; 针对上述区域分解算法,探索多体系统并行异步算法, 不同子区域可采用独立的时间积分步长.

(5) 基于时空后验误差估计的自适应算法.针对基于SE(3)群描述的 $\alpha$类时间积分算法以及保能量-(角)动量时间积分算法,提出时间离散后验误差估计方法, 并实现变步长时间积分算法;借鉴已有的有限元方法, 采用超收敛修复类或残差类方法进行空间离散后验误差估计,借鉴等几何分析中分层基函数概念实现空间网格自适应;针对时空统一离散的变分积分子, 借鉴时空有限元方法发展成果,进行时空后验误差估计及自适应算法.

(6) 多柔体系统动力学降阶算法.研究基于数据驱动的参数化多柔体系统动力学降阶方法^[36],探索局部标架建模方法在构造降阶基矢量以及降阶模型动力学仿真的优势. 例如,消除刚体运动的几何非线性, 多柔体系统的降阶基矢量仅需描述柔体弹性变形部分,这将有望减少降阶模型的基矢量个数. 同时,多柔体系统的切线刚度矩阵更新次数降低,能有效提高降阶模型的动力学仿真计算效率. 此外, 通过对比局部标架运动规律,研究单元局部标架自适应合并与分离方法, 有望实现柔体的自适应降阶算法.当柔性体合并至唯一局部标架时, 与传统浮动坐标降阶方法实现兼容.

4 数值算例

本节通过若干数值算例, 阐述上述基于SE(3)群局部标架的多柔体系统动力学建模和计算方法的可行性.

4.1 多柔体系统保能量(角)动量算法

对于线性系统, 当时间积分算法谱半径$\rho ≤ $1时, 该算法无条件稳定.但上述性质无法直接推广至非线性系统,评价总能量是否守恒是评价非线性系统时间积分算法稳定性的重要指标之一.而若采用传统高频耗散的a类算法, 例如HHT或广义a方法,可能会引起非线性系统能量增加, 引起时间积分算法的发散. 因此,对于柔性多体系统动力学仿真, 开发保能量-(角)动量算法具有极其重要意义.

Simo等^[37]首次针对 $\mathbb{R}^{3}$ $\times$ SO(3)群几何精确梁模型, 提出了保能量-(角)动量时间积分方法.尽管随后研究表明^[38]该算法并不能够精确保持系统能量与角动量,但上述方法为非线性系统, 尤其是基于SO(3)群的梁板壳有限单元,构造保能量-(角)动量方法提供了基本思路.本小节以SE(3)群局部标架的几何精确梁单元为例,给出一种保能量-(角)动量时间积分算法.并通过含线性/非线性约束、小/大变形的空间双摆动力学模型,分析其处理多体系统动力学问题的数值特性.

(1) 考察含线性约束的双摆模型的小变形动力学问题. 如图3(a)所示,双摆的两根杆分别长0.3 m和0.5 m, 其宽度与高度均为0.005 m, 材料参数为: $\rho=2700$ kg/m$^{3}$, $E=7.0 \times 10^{10}$ Pa.杆I与地面以及两杆间均通过球铰连接, 系统约束方程均为线性约束. 在初始时刻,杆I静止放置于$X$轴上, 杆II绕$B$点存在初始角速度[3 0 $-4$]$^{T}$rad/s.系统不受外力作用, 在初始速度作用下开始运动.通过基于SE(3)群局部标架的几何精确梁单元对双摆系统进行离散,采用保能量-(角)动量时间积分算法求解系统动力学方程, 仿真时间为100 s,时间步长为$h=1.0 \times 10^{-4}$ s, 总仿真步数为10$^{6}$.

图3

新窗口打开|下载原图ZIP|生成PPT
图3柔性双摆系统的初始构型: (a) 线性约束; (b) 非线性约束

Fig.3Initial configuration of a flexible double pendulum: (a) Linear constrains; (b) nonlinear constrains



图4给出了双摆模型在前30 s内末端$C$点的位置随时间变化规律,其余70 s内位置变化趋势大体类似. 由图可见,小变形双摆模型长时间历程动力学响应非常复杂, 运动非线性程度较高.由于强非线性系统对初值异常敏感, 难以评价算法精度,而系统守恒量能够作为评价算法精度的一个重要指标.由于双摆系统仅在与地面连接处受到约束力作用,该系统对惯性坐标系原点的角动量以及系统总能量守恒.图5与图6分别给出了系统能量、系统沿三个方向角动量的演化规律,其中曲线上下方的数字表示该物理量波动的最大值与最小值. 本质上来说,含线性约束多体系统在计算过程中可等效为一个无约束系统,半离散平衡方程为一常微分方程组. 对于此类情况,时间积分算法可精确地保持系统能量及角动量, 系统能量的相对误差仅为$5.0 \times 10^{-12}$, 系统角动量的相对误差为$7.0 \times 10^{-10}$.

图4

新窗口打开|下载原图ZIP|生成PPT
图4含线性约束双摆模型在小变形运动中末端$C$点位置

Fig.4Displacement of end point $C$ of the double pendulum with linear constrains subject to small deformation and overall motion

图5

新窗口打开|下载原图ZIP|生成PPT
图5含线性约束双摆模型在小变形运动中的总能量变化

Fig.5Total energy variation of the double pendulum with linear constrains subject to small deformation and overall motion

图6

新窗口打开|下载原图ZIP|生成PPT
图6含线性约束双摆模型在小变形运动中的角动量变化(2) 为了考察含线性约束的双摆模型在大变形运动中的角动量及能量守恒特性,将(1)模型中两根杆的截面宽度与高度均改为0.000 2 m, 仿真时间为10 s,时间步长为$h=1.0 \times 10^{-5}$ s, 总仿真步仍取为10$^{6}$步. 此时,计算得到的最大变形率达到15%. 图7与图8

Fig.6Angular momentum variation of the double pendulum with linear constrains subject to small deformation and overall motion

图7

新窗口打开|下载原图ZIP|生成PPT
图7含线性约束双摆模型在大变形运动中的总能量变化

Fig.7Total energy variation of the double pendulum with linear constrains subject to large deformation and overall motion

图8

新窗口打开|下载原图ZIP|生成PPT
图8含线性约束双摆模型在大变形运动中的角动量变化分别给出系统能量与角动量的时间历程.由图可见, 该算法同样能精确地保持系统的守恒量,其中能量相对误差为$1.7 \times 10^{-5}$, 角动量相对误差为$1.0 \times 10^{-5}$. 从理论上来看, 该系统的能量严格守恒.此处的误差主要来自求解非线性代数方程,且该误差随着非线性方程收敛标准的降低而减少. 此外, 从能量的时间历程可见,系统的低频振动激发起高频振动响应, 而有限元模型无法精确描述高频振动模态.同时, 系统保留的大量高频响应将导致计算过程中需要非常小的时间步长,导致计算效率低下. 因此, 对于多柔体系统动力学仿真, 有必要研究高频能量耗散可控的时间积分算法.

Fig.8Angular momentum variation of the double pendulum with linear constrains subject to large deformation and overall motion



(3) 本小节考察含非线性约束的保能量-(角)动量算法的数值特性. 此时,半离散系统动力学方程为一指标三的微分代数方程组,非线性约束的引入将对多体系统动力学数值特性带来显著的改变.如图3(b)所示, 两杆间由旋转铰相连接, 为一典型的含非线性约束多体系统.双摆模型两杆的几何及材料属性与模型(1)一致, 杆II绕$B$点存在初始角速度[300]$^{T}$ rad/s. 仿真时间为100 s, 时间步长为$h=1.0 \times 10^{-4}$ s,总仿真步仍取为10$^{6}$步. 图9与图10分别给出系统总能量与角动量的时间历程.由图可见, 系统能量相对误差为$4.9 \times 10^{-9}$,角动量相对误差为$9.1 \times 10^{-3}$.非线性约束对系统总能量的守恒特性几乎无影响. 然而,这对系统角动量守恒精度带来了较大的影响.这是由于直接求解指标三的含非线性约束微分代数方程组时,需要引入$h^{2}$量级的缩放系数, 导致平衡方程中广义惯性力与广义约束力数值精度差$h^{2}$量级.而系统角动量误差与广义惯性力、广义弹性力、广义约束力的数值精度以及当前构型位置相关.随着时间积分步长的累计, 角动量计算过程中相对误差仅为$9.1 \times 10^{-3}$, 但数值精度仍在可接受范围内. 理论上, 可通过降微分代数方程指标技术, 提高角动量守恒精度. 但在实际工程问题中, 需统筹考虑计算效率与计算精度. 此外, 对于其他类型的非线性约束, 都能够构造相应的广义约束力, 使得系统能量守恒.由于保能量-(角)动量算法的构造依赖于运动参数化方法、系统平衡、几何、本构以及约束方程等一些细节因素, 针对其他类型局部标架单元的保能量-(角)动量算法值得进一步研究. 需要指出的是, 对于非保守系统, 可将耗散力所做的功作为系统一部分能量, 从数值角度亦可其视为一个守恒系统, 系统能量依然是一个守恒量. 例如, 4.4小节所涉及的摩擦碰撞动力学问题, 保能量-(角)动量算法能够精确表征法向碰撞势能与切向摩檫力耗散能量,整体考虑总能量依然守恒.

图9

新窗口打开|下载原图ZIP|生成PPT
图9含非线性约束双摆模型在大变形运动中的总能量变化

Fig.9Total energy variation of the double pendulum with nonlinear constrains subject to large deformation and overall motion

图10

新窗口打开|下载原图ZIP|生成PPT
图10含非线性双摆模型在大变形运动中的角动量变化此外, 在计算数学领域, Marsden等^[39]提出了变分积分概念,该算法被证明能在长时间动力学仿真过程中保持系统的辛几何结构. 对于保守系统,离散系统的动量与角动量守恒, 离散系统的能量在守恒量附近波动. 在欧氏空间内,变分积分算法与上述保能量-(角)动量算法之间的对比可参考文献[40].从工程应用角度来看, 李群SE(3)的变分积分算法仍需要进一步发展.

Fig.10Angular momentum variation of the double pendulum with nonlinear constrains subject to large deformation and overall motion

4.2 广义a方法后验误差估计方法

本小节通过一个自由刚体运动问题和刚性双摆运动问题, 验证基于SE(3)群描述的广义a方法误差估计方法的有效性.

首先, 考察自由刚体的运动. 设初始时刻, 刚体质心位于惯性坐标系原点,刚体的连体坐标系与惯性坐标系重合. 刚体在连体坐标系下的主转动惯量分别为3783.42 kgm$^{2}$, 3776.16 kgm$^{2}$, 2630.94 kgm$^{2}$, 惯性积分别为352.9 kgm$^{2}$, 484.74 kgm$^{2}$, 469.74kgm$^{2}$. 初始时刻,刚体质心速度为[10\ 10\ 10]$^{T}$ m/s, 刚体角速度为[10\ 10\ 10]$^{T}$ rad/s. 此外,刚体在质心处受到[50$|$sin(10$t)$$|$ 52$|$sin (10$t+\pi $/3)$|$ 50$|$sin (10$t+\pi $/2)$|$]$^{T}$ kN的集中力作用,同时受到[50$|$sin(10$t)$$|$ 52$|$sin (10$t+\pi $/3)$|$ 50$|$sin (10$t+\pi $/2)$|$]$^{T}$ kN$\cdot$m的集中力矩作用.现采用基于SE(3)群描述的广义a方法进行仿真计算, 总仿真时间为1 s,时间步长为10$^{-3}$ s, 算法谱半径为0.8.

本算例将每一时间步长细化10倍的数值仿真结果视为精确解, 由此得到时间离散截断后验误差的参考值. 将参考误差与后验误差结果进行对比,以验证基于SE(3)群描述的广义a方法误差估计方法的正确性. 图11给出了刚体位移与转动的时间离散误差估计值与参考值的二范数. 由于本算例的动力学方程具有高度光滑性, 误差估计值与参考解极其吻合.

图11

新窗口打开|下载原图ZIP|生成PPT
图11基于SE(3)群描述的广义a方法计算单刚体运动的后验误差估计其次, 考察刚性双摆系统, 其几何参数与材料属性与4.1算例(3)中模型完全一致. 设该系统在沿着$Z$轴负方向的重力作用下开始运动.采用基于SE(3)群描述的广义a方法进行仿真计算, 总仿真时间为1 s, 时间步长为$1.0 \times 10^{-3}$ s, 算法参数谱半径为0.8.图12与图13分别给出了杆OB的估计误差与参考误差的二范数. 由于后验误差计算方法是基于常微分方程推导, 而多柔体系统约束的引入, 导致估计误差与参考误差之间存在一定的差异, 但两者的误差量级与整体趋势一致, 通过估计误差进行变步长时间积分算法设计是可行的.

Fig.11Posterior error estimation of the generalized alpha method of SE(3) group for a single rigid body

图12

新窗口打开|下载原图ZIP|生成PPT
图12基于SE(3)群的广义a方法对刚性双摆模型平动位移的后验误差估计

Fig.12Posterior error estimation of the generalized alpha method of SE(3) group for the translational displacement of a double pendulum

图13

新窗口打开|下载原图ZIP|生成PPT
图13基于SE(3)群的广义a方法对刚性双摆模型转动位移的后验误差估计上述基于SE(3)群描述的广义a方法的误差估计方法可直接推广至多柔体系统的动力学仿真中. 但随着系统刚度矩阵与阻尼矩阵的引入,时间离散误差估计值与参考误差将进一步扩大, 但两者的变化趋势依然是一致的.

Fig.13Posterior error estimation of the generalized alpha Lie group SE(3) method for the rotational displacement of a double pendulum

4.3 局部标架壳单元静力学分析

本小节采用文献[41]中圆柱壳和球壳静力学模型, 验证本文所提出的SE(3)群局部标架板壳单元的正确性.

如图14所示, 一个无约束的圆柱壳在$A$, $E$两点受到一对大小相等、方向相反的集中力$F=40$ kN的作用,其中$A$点为圆柱体上母线的中点, $E$点为下母线的中点. 圆柱壳的长度$L$为10.35 mm,内径$R$与厚度分别为4.953 mm与0.094 mm, 材料杨氏模量$E=10.5 \times 10^{6}$ N/mm$^{2}$, 泊松比$\upsilon =0.312 5$. 考虑到圆柱壳模型的对称性, 只需要对该模型上半部分的四分之一($ABCD$)进行建模.现将在圆柱壳上施加外力F的过程平均分成100个增量步, 并设静力学迭代容许误差为$1.0 \times 10^{-6}$.

图14

新窗口打开|下载原图ZIP|生成PPT
图14受拉载荷的圆柱壳初始构形本小节分别采用局部标架GESF壳单元与局部标架缩减ANCF壳单元离散$ABCD$圆柱壳, 并分析了24 $\times$ 16与36 $\times$ 24的两种网格数目下圆柱壳的静力学变形. 图15中给出了壳体上$A$, $B$和$C$三点的位移随外力增加的变化曲线. 当载荷接近20 kN时, 圆柱壳发生了屈曲变形, 屈曲前圆柱壳以弯曲变形为主, 屈曲后圆柱壳以拉伸变形为主. 因此, 本算例中GEF壳单元不能对忽略薄膜应变对弯曲变形的影响. 数值结果表明, 上述两类壳单元数值结果与参考文献均吻合较好, 均能够精确地描述圆柱壳的大变形特性. 关于一般性弹性体的屈曲、后屈曲及分叉模拟算法可参见作者前期工作^[42].

Fig.14Initial configuration of a pinched cylinder

图15

新窗口打开|下载原图ZIP|生成PPT
图15圆柱壳上$A$, $B$, $C$点的位移大小其次, 本节对图16所示的顶部存在18$^{\circ}$开口薄球壳进行静力学计算,以说明SE(3)局部标架几何精确壳单元同样能够适用于描述薄壳结构的大变形.该球壳的半径为 10.0 mm, 杨氏模量为 $E=6.825 \times 10^{7}$ N/mm$^{2}$,泊松比为 $\upsilon =0.3$, 厚度为 $t=0.01$ mm($R/t=1000$). 该球壳在底部的$A$, $B$, $C$, $D$四点分别受到集中力 $F=2\lambda F_{0}$ 的作用, 其中$F_{0}=1.0$ N, $\lambda $取为 2.5. 考虑到球壳的对称性,现对其四分之一的模型进行建模.本小节分别采用局部标架GESF壳单元与局部标架缩减ANCF壳单元离散四分之一开口球壳,并分析了24 $\times$ 24与32 $\times$ 32的两种网格数目下球壳的静力学变形.图17给出了不同单元数目下球壳上点$A$与点$B$的位移随外载荷变化曲线. 数值结果表明,上述两类壳单元数值结果与参考文献均吻合较好,本文提出的两类单元能够精确描述薄壳结构的大变形特性.

Fig.15Magnitudes of displacements at nodes $A$, $B$ and $C$ of a pinched cylinder

图16

新窗口打开|下载原图ZIP|生成PPT
图16顶部18$^{\circ}$ 开口球壳初始构形

Fig.16Initial configuration of a hemispherical shell with 18$^{\circ}$ hole

图17

新窗口打开|下载原图ZIP|生成PPT
图17开口球壳上$A$、$B$点的位移大小考虑到板壳单元惯性力与梁单元惯性力计算方法类似,本节中不再加以分析板壳结构的动力学问题. 由于局部标架的引入,板壳结构的单元质量矩阵与切线刚度矩阵满足刚体运动的不变性. 在数值计算中,六自由度的板壳单元质量矩阵与切线刚度矩阵的更新次数将急剧减少,具体数值性质可参考4.5小节中的大型桁架结构展开动力学模拟.

Fig.17Magnitudes of displacements at nodes $A$ and $B$ of a hemispherical shell

4.4 多柔体系统碰撞动力学仿真

本小节通过空间飞网抓捕收口动力学算例以及刚性杆盒碰撞动力学算例,展示基于SE(3)群局部标架描述的多柔体系统与多刚体系统碰撞动力学算法.

首先, 将基于SE(3)群局部标架的柔性梁建模方法与罚函数算法相结合, 处理碰撞问题.

图18为正六边形的空间飞网,其在完全展开状态下以一定的速度抓捕正前方的固定星体. 当抓捕至适当位置时,通过收口装置与收口绳进行缩紧网口.现通过基于SE(3)群局部标架的DER单元对飞网系统进行有限元离散,对飞网与刚性星体、收口绳与刚性环间的碰撞,采用点-线接触策略与保能量罚函数算法进行模拟,利用4.1节提出的保能量-(角)动量时间积分算法求解系统动力学方程.

图19中给出了飞网抓捕星体以及收口阶段的示意图. 在整个抓捕与收口阶段,该系统总能量、动量以及角动量均守恒. 计算结果表明,罚函数方法能够较好地处理大变形柔体间的碰撞问题.若进一步考虑线-线接触、线-面接触问题, 可采用Mortar方法细化接触区域,从而得到高精度的碰撞模型, 其具体算法细节可见文献[43]. 需要指出的是,对于此类强非线性问题, 若采用传统算法(如Newmark-$\beta $、广义a方法等),即使谱半径为1, 系统能量也无法守恒, 甚至会出现总能量增大.

图18

新窗口打开|下载原图ZIP|生成PPT
图18大型空间飞网展开状态及被捕获物体示意图

Fig.18Deployed configuration of the large tethered-net and the target

图19

新窗口打开|下载原图ZIP|生成PPT
图19大型空间飞网抓捕与收口动力学模拟其次, 考察基于SE(3)群局部标架的建模方法与Lagrange乘子法相结合,解决多体系统接碰撞动力学问题. 此时,通过将碰撞时的单边约束表示为一类互补条件, 可通过互补算法求解接触冲量.以下通过刚性杆与空心盒间的两点碰撞算例,说明Lagrange乘子法可直接与上述基于SE(3)群局部标架的方法相兼容,并与罚函数法进行比较.

Fig.19Capturing and closing dynamics of the large tethered-net



如图20所示, 正方形截面杆的边长为0.2 m.刚性空心盒的外表面长、宽、高分别为1.8 m, 1.6 m, 1.4 m, 厚度为0.1 m.杆与空心盒的密度均为7850 kg/m$^{3}$. 在初始时刻,杆两端点在惯性坐标系下的位置坐标分别为$[0 -0.2 1]^{T}$ m与$[1 0.2 0]^{T}$ m. 空心盒不受重力, 杆在沿Z轴负方向的重力作用下开始运动.现采用SE(3)群刚体单元描述刚性杆与空心盒, 并通过非光滑广义$\alpha$锥互补方法对摩擦系数为0, 0.2, 0.4以及0.6的4种情况进行动力学仿真,算法参数谱半径设为0.8, 其中非光滑广义a锥互补方法可参考文献[44].为了验证该算法的正确性, 同时采用罚函数方法对上述几种情况进行数值仿真.

图20

新窗口打开|下载原图ZIP|生成PPT
图20刚性杆与空心盒碰撞模型图21分别给出了摩擦系数为0, 0.2, 0.4以及0.6时,罚函数方法与非光滑广义$\alpha $锥互补方法得到的杆上端点$z$方向位移变化曲线.计算结果表明, 上述两类算法在无摩擦时的位移完全一致. 当存在摩擦时,两中方法存在微小差异. 这是由于锥互补条件推导过程中引入了放松条件,导致该方法处理滑动摩擦时存在小的误差. 上述研究表明,在SE(3)群内构造的两种算法均能处理多柔体系统的多点碰撞问题. 但两者相比较,罚函数方法较为简单, 能方便处理点-面接触或面-面接触问题;而非光滑互补算法相对复杂, 应用于大变形柔性间接触问题的难度较大. 然而,在处理刚性碰撞问题时, 由于罚参数的引入将在系统中引起高频振动,而采用互补条件描述的非光滑动力学能较为精确地处理刚性碰撞问题.

Fig.20Contact model between a rigid rod and a hollow box



在本算例中, 上述两种算法均与隐式时间积分算法相结合.

罚函数方法与锥互补算法同样能应用于显式时间积分算法, 将其与异步变分积分算法结合^[45], 有望能处理复杂多柔体系统的多尺度碰撞问题, 但此类算法的计算精度与计算效率仍需进一步研究.

图21

新窗口打开|下载原图ZIP|生成PPT
图21不同摩擦系数时锥互补方法与罚函数方法对比

Fig.21Comparison between the CCP method and penalty method for the different friction coefficients

4.5 大型桁架结构展开动力学模拟

本小节通过一种空间桁架结构展开动力学模拟, 展示基于SE(3)群局部标架的建模方法在计算效率方面的优势. 图22所示为一单胞元可展开空间桁架结构, 可用于构建大型空间结构系统.

图22

新窗口打开|下载原图ZIP|生成PPT
图22单胞元空间桁架结构展开示意图^[46]

Fig.22The deployed configuration of a unit of the space truss structure^[46]



该胞元结构由28根杆通过旋转铰与滑移铰连接而成, 其具体描述见文献[46].胞元的截面为边长2 m的正方形, 横杆长度为4 m, 杆为空心圆截面,材料的拉伸模量为$2.1 \times 10^{11}$ Pa, 泊松比为0.3,密度为7850 kg/m$^{3}$. 设桁架结构在端部与地面固连, 各胞元通过控制角度$\beta$同时驱动展开, 展开时间设为20 s, 展开角的驱动为余弦控制函数. 在初始时刻,展开角度与展开锁定时刻角度分别为$\beta_{0}=5^{\circ}$与$\beta_{f}=90^{\circ}$. 对每个桁架胞元, 采用10个基于SE(3)群局部标架的几何精确梁建模,用基于SE(3)群描述的广义a方法进行动力学仿真, 算法谱半径为0.8.

首先, 进行两个桁架胞元的展开动力学仿真,分析展开动力学模拟中系统迭代矩阵的复用情况,以说明基于SE(3)群局部标架的建模方法能有效消除刚体运动的非线性,极大地降低整个系统的几何非线性. 当梁截面内外直径分别选为0.055 m与0.08 m时,桁架胞元在展开过程中几乎不发生弹性变形. 此时,系统广义质量矩阵与广义刚度矩阵均无需更新,仅根据初始构型下的初值就能完成整个展开过程的动力学仿真. 这意味着,此时无须计算梁单元的几何刚度矩阵, 即系统迭代矩阵为对称矩阵.

为了考察具有大变形的多柔体系统动力学仿真过程具有上述数值特性,将梁截面内外直径设为0.023 m与0.024 m, 进行同样工况的动力学仿真.此时桁架胞元的最大变形率达到3.3%, 在其展开过程中能够观察到明显的弹性振动.

表1中给出了两种不同时间步长所对应的迭代矩阵$K$与 $\varPhi_{q}$的更新次数以及总的牛顿迭代次数. 其中, 矩阵$K$与系统广义质量矩阵与切线刚度矩阵相关, $\varPhi_{q}$ 为约束方程的雅可比矩阵, "$NNI$ $\geqslant i$"表示一个时间步内, 单个时间步内牛顿迭代次数大于等于$i$步时, 迭代矩阵强制更新; "Frozen"表示迭代矩阵永不更新. 当牛顿迭代步数超过预设步数时, 优先更新约束方程的雅可比矩阵.


Table 1
Table 1The number of times that the iteration matrix and Jacobian of constrain equations are updated as well as the number of total Newton iterations

新窗口打开|下载CSV

首先考察时间步长为$2.0 \times 10^{-3}$ s时的工况. 当$NNI$ $\geqslant $ 1时,即每个时间步中牛顿迭代均更新迭代矩阵, 共需进行32 410次牛顿迭代,平均每个时间步需要3.2次牛顿迭代; 当$NNI$ $\geqslant $ 4时,即若每个时间步中牛顿迭代达到四步时更新一次迭代矩阵, 则迭代矩阵仅需更新19次,约束方程雅可比矩阵需更新3083次, 且总迭代次数较之前增加32%.而当迭代矩阵不更新时, 总迭代次数比$NNI$ $\geqslant $ 4时略有增大. 由此可说明,基于SE(3)群局部标架的建模方法能极大地降低刚体运动导致的几何非线性,迭代矩阵可视为一个常数矩阵.

再考察时间步长进一步减少的情况. 此时, 系统迭代矩阵中的广义质量矩阵占比提升.一般来说, 柔体的广义质量矩阵与弹性变形相关, 而切向刚度矩阵与变形梯度相关.因此, 相比于切向刚度矩阵, 广义质量矩阵往往变化较小. 因此,时间步长为$1.0 \times 10^{-3}$ s, $NNI$ $\geqslant $ 4时, 迭代矩阵几乎不用更新.同时, 约束方程雅可比矩阵次数也大幅降低. 这是由于迭代矩阵趋于常数矩阵,从而减少了对约束方程雅可比矩阵的更新次数.

在大规模的多柔体系统动力学隐式算法仿真中,最耗时的两部分是计算系统迭代矩阵和求解高维线性代数方程组.基于SE(3)群局部标架的建模方法能最大限度地降低系统迭代矩阵更新次数,由此可将迭代矩阵提前进行计算并进行LU分解, 或用于设计迭代算法的预条件算子.相比于其他建模方法, 基于SE(3)群局部标架的建模方法的计算效率显著提高.

最后, 基于上述建模方法, 对具有8个胞元的桁架结构进行展开动力学模拟.图23给出了不同时刻该桁架结构的展开构型.

图23

新窗口打开|下载原图ZIP|生成PPT
图23具有8个胞元的空间桁架展开动力学模拟

Fig.23Deployment simulation of a truss structure with eight units

4.6 基于区域分解的多柔体系统动力学通用并行算法

本小节通过环形桁架-索网天线与大型空间飞网的展开动力学问题,展示一种基于有限元网格撕裂对接(finite element tearing and interconnecting dual-primal, FETI-DP)的多柔体系统动力学通用并行算法. 其中,环形桁架-索网天线展开动力学算例是为了说明并行算法能够适用于含多种不同类型单元与多种不同运动副的复杂多体系统动力学计算;大型空间飞网展开动力学算例是为了并行算法能够处理百万广义坐标量级的柔性多体系统动力学仿真.在该方法中, 采用Dirichlet预条件算子迭代来并行求解界面问题.该方法已成功应用于求解大规模结构静/动力学问题, 具体算法见文献[47,48].

图24是由中国空间技术研究西安分院与北京理工大学联合研制的4 m直径环形桁架-索网天线模型,其几何参数和材料参数详见文献[49].采用基于SE(3)群局部标架的几何精确梁单元和DER单元分别对环形桁架和索网进行建模,将桁架划分为18个单元, 将索网划分为20个单元, 共5.3万个自由度.在该多柔体系统动力学模型中, 存在众多不同类型的运动副,如球铰、旋转铰、滑移铰以及齿轮副. 如图25所示, 通过多层图分解技术,将其有限元网格分解为16区域, 对不同区域通过不同颜色进行表征,实线表示索网系统, 虚线表示环形桁架. 经区域分解后, 每个区约含3800个自由度,在区域界面上约有1000个自由度.

图24

新窗口打开|下载原图ZIP|生成PPT
图244 m直径的环形桁架-索网天线^[49]: (a) 地面实验系统;(b)动力学模拟过程

Fig.24Deployable hoop truss antenna of 4 m in diameter^[49]:(a) Experiment system; (b) dynamic simulation

图25

新窗口打开|下载原图ZIP|生成PPT
图25环形桁架-索网天线的区域分解示意图

Fig.25Domain decomposition of the deployable hoop truss antenna



由于在作者的前期研究中已系统分析过该天线的展开动力学特性^[50], 此处仅给出索网天线在无重力环境下的收回动力学过程,但仿真过程并没有考虑索网间以及与桁架间的碰撞问题.图26给出了不同时刻该天线的构型图,验证了上述并行计算方法可适用于含复杂约束的多柔体系统动力学仿真.

图26

新窗口打开|下载原图ZIP|生成PPT
图26基于区域分解的环形桁架-索网天线收回动力学模拟

Fig.26Folding simulation of the hoop truss antenna based on the domain decomposition



最后, 通过对大型空间飞网系统进行展开动力学分析,说明该并行算法能够处理自由度达百万量级规模的多柔体系统大变形动力学仿真问题.本算例的空间飞网拓扑结构与4.4节中飞网一致, 但是增加了几何尺寸以及网目个数.对于此类超柔性的多体系统,只有采用较为稠密的网格才能较精确地描述其动力学特性.本研究通过基于SE(3)群局部标架的DER单元对该空间飞网系统进行离散,将每一段绳索离散为80个DER单元, 系统自由度达到202万.由于此处仅关心并行算法的可行性,故不再详细介绍绳索的几何参数、材料参数以及具体展开工况. 通过多层图分解技术,将该系统的有限元网格分解为72个区域, 每个区域具有2.8万个自由度,区域界面上具有0.5万个自由度. 图27给出了空间飞网在不同时刻的展开构型.

图27

新窗口打开|下载原图ZIP|生成PPT
图27大型空间飞网展开动力学的并行模拟需要指出的是, 本算法在处理界面问题时的效率仍有较大提升空间.在粗网格构造以及界面问题高效预条件算子等方面, 还需要进入深入研究. 目前,仅实现了空间梁结构的展开动力学并行计算. 对于大型薄膜结构的展开动力学问题,还值得进一步研究.

Fig.27Parallel simulation of the deployment of a large tethered-net based on the domain decomposition

5 结束语

本文基于SE(3)局部标架,讨论了如何发展一套新的多柔体系统动力学建模与计算方法体系.该方法体系在SE(3)局部标架下描述单元应变、应力、运动参数、运动参数变分及增量等物理量,可有效避免大范围刚体运动导致的几何非线性. 对于数值计算而言,该方法体系具有如下优势: 通过最少转动参数(三参数)可完全避免转动奇异性问题;系统广义质量矩阵的主项为常数, 满足刚体转动的不变性,而且对大多数问题广义质量矩阵无须更新;单元切线刚度矩阵能满足刚体转动的不变性,可最大限度地减少切线刚度矩阵的更新次数, 有效提高计算效率. 在计算精度方面,该方法体系能继承欧氏空间中有限单元的性质,即已有提高收敛性的单元技术均能直接应用于由SE(3)局部标架描述的有限单元;由于在该局部标架下描述柔体运动, 可自然消去部分自由度来满足约束方程,例如构造五自由度的板壳单元.

上述新的多柔体系统动力学建模和计算方法体系包含如下核心方法和算法:基于SE(3)群局部标架的复杂多柔体系统建模方法;多柔体系统动力学的长时间历程时间积分器; 多柔体系统的多尺度碰撞算法;多柔体系统的分布式通用并行异步算法; 基于时空后验误差估计的自适应算法;多柔体系统的动力学降阶算法等. 本文通过若干简单算例和工程案例,简要介绍了部分算法及其可行性.

目前, 对上述方法和算法的研究仍处于初步阶段, 尚需进行完善并使其相互融合,进而发展成为一套新的多柔体系统动力学建模和计算方法体系,并形成相应的柔性多体系统动力学软件.

参考文献 View Option 原文顺序
文献年度倒序
文中引用次数倒序
被引期刊影响因子

[1] Simo JC, Vu-Quoc L . On the Dynamics in space of rods undergoing large motions-a geometrically exact approach
Computer Methods in Applied Mechanics and Engineering, 1988, 66 (2): 125-161
[本文引用: 3]

Shabana AA . An absolute nodal coordinates formulation for the large rotation and deformation analysis of flexible bodies.
Chicago: University of Illinois at Chicago, 1996
[本文引用: 3]

[3] Meier C, Popp A, Wall WA . Geometrically exact finite element formulations for slender beams: Kirchhoff-Love theory versus Simo-Reissner theory
Archives of Computational Methods in Engineering, 2019,26(1):163-243
[本文引用: 2]

[4] Gerstmayr J, Sugiyama H, Mikkola A . Review on the absolute nodal coordinate formulation for large deformation analysis of multibody systems
Journal of Computational and Nonlinear Dynamics, 2013,8(3):031016
[本文引用: 1]

[5] Romero I . A comparison of finite elements for nonlinear beams: The absolute nodal coordinate and geometrically exact formulations
Multibody System Dynamics, 2008,20(1):51-68
[本文引用: 1]

[6] Bauchau OA, Han SL, Mikkola A , et al. Comparison of the absolute nodal coordinate and geometrically exact formulations for beams
Multibody System Dynamics, 2014,32(1):67-85
[本文引用: 1]

[7] 王珍, 刘铖 . 几何精确方法与绝对节点坐标方法对比研究//第十届全国多体动力学与控制暨第五届全国航天动力学与控制学术会议, 青岛 , 2017 年9月22-25日
( Wang Zhen, Liu Cheng. Comparison of the geometrically exact formulation and the absolute nodal coordinate formulation//The 10th National Conference on Multibody dynamics and Control & the 5th National Conference on Astrodynamics and Control, Qingdao, 2017-9-22-25 (in Chinese))
[本文引用: 1]

[8] Reissner E . On one-dimensional finite-strain beam theory: The plane problem
Zeitschrift F$\ddot{u}$r Angewandte Mathematik und Physik Zamp, 1972,23(5):795-804
[本文引用: 1]

[9] Ren H, Fan W, Zhu WD . An accurate and robust geometrically exact curved beam formulation for multibody dynamic analysis
ASME Journal of Vibration and Acoustics, 2017,140(1):011012
[本文引用: 1]

[10] Liu JP, Cheng ZB, Ren GX . An arbitrary Lagrangian--Eulerian formulation of a geometrically exact Timoshenko beam running through a tube
Acta Mechanica, 2018,229(8):3161-3188
[本文引用: 1]

[11] Cardona A , G$\acute{e}$radin M. A beam finite element non-linear theory with finite rotations
International Journal for Numerical Methods in Engineering, 1988,26:2403-2438
[本文引用: 2]

[12] Sonneville V, Cardona A , Br$\ddot{u}$ls O. Geometrically exact beam finite element formulated on the special Euclidean group SE(3)
Computer Methods in Applied Mechanics and Engineering, 2014,268:451-474
[本文引用: 1]

[13] 刘铖, 胡海岩. 浅谈多体系统动力学几何非线性与转动非线性// 第十一届全国多体动力学与控制暨第六届全国航天动力学与控制和第十四届全国分析力学联合学术会议, 长沙
2019 -09-20-23日.
[本文引用: 1]

Liu Cheng, Hu Haiyan. A brief discussion on the geometrically nonlinear and rotational nonlinear of the multibody systems//The 11th National Conference on Multibody dynamics and Control & the 6th National Conference on Astrodynamics and Control & the 14th National Conference on analytical dynamics, Changsha
2019-09-20-23 (in Chinese))
[本文引用: 1]

[14] Hughes TJR, Cottrell JA, Bazilevs Y . Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement
Computer Methods in Applied Mechanics and Engineering, 2005,194(39-41):4135-4195
[本文引用: 2]

[15] Simo JC, Fox DD . On a stress resultant geometrically exact shell model Part I: Formulation and optimal parametrization
Computer Methods in Applied Mechanics and Engineering, 1989,72(3):267-304
[本文引用: 2]

[16] Simo JC, Rifai MS, Fox DD . On a stress resultant geometrically exact shell model Part IV: Variable thickness shells with through-the-thickness stretching
Computer Methods in Applied Mechanics and Engineering, 1990,81(1):91-126
[本文引用: 1]

[17] Simo JC, Kennedy JG . On a stress resultant geometrically exact shell model Part V: Nonlinear plasticity formulation and integration algorithms
Computer Methods in Applied Mechanics and Engineering, 1992,96:133-171
[本文引用: 2]

[18] Simo JC, Fox DD, Hughes TJR . Formulations of finite elasticity with independent rotations
Computer Methods in Applied Mechanics and Engineering, 1992,95(2):277-288
[本文引用: 1]

[19] Fox DD, Simo JC . A drill rotation formulation for geometrically exact shells
Computer Methods in Applied Mechanics and Engineering, 1992,98(3):329-343
[本文引用: 1]

[20] 潘亦甦, 陈大鹏 . 具有Drilling自由度的膜单元的杂交／混合有限元法
西南交通大学学报, 1995,30(2):128-134.
[本文引用: 1]

( Pan Yisu, Chen Dapeng . A hybrid/mixed model for membrane elements with drilling degrees of freedom
Journal of Southwest Jiaotong University, 1995,30(2):128-134 (in Chinese))
[本文引用: 1]

[21] Sonneville V . A geometric local frame approach for flexible multibody systems.
Liege, Belgique: University of Liege, 2015
[本文引用: 4]

[22] Shabana AA, Yakoub RY . Three-dimensional absolute nodal coordinate formulation for beam elements: Theory
ASME Journal of Mechanical Design, 2001,123(4):606-613
[本文引用: 1]

[23] Yakoub RY, Shabana AA . Three-dimensional absolute nodal coordinate formulation for beam elements implementation and applications
ASME Journal of Mechanical Design, 2001,123:614-621
[本文引用: 1]

[24] Garc$\acute{i}$a De Jal\'{o}n J, Bayo E . Kinematic and Dynamic Simulation of Multibody Systems the Real-Time Challenge
New York: Springer, 1994
[本文引用: 1]

[25] Bauchau OA . Flexible Multibody Dynamics
Springer Science & Business Media, 2010
[本文引用: 1]

[26] Jeleni$\acute{c}$ G, Cristfield MA . Geometrically exact 3D beam theory: implementation of a strain-invariant finite element for statics and dynamics
Computer Methods in Applied Mechanics and Engineering, 1999,171(1):141-171
[本文引用: 1]

[27] Ibrahimbegovic A, Taylor RL . On the role of frame-invariance in structural mechanics models at finite rotations
Computer Methods in Applied Mechanics and Engineering, 2002,191(45):5159-5176
[本文引用: 1]

[28] Bauchau OA, Han SL . Interpolation of rotation and motion
Multibody System Dynamics, 2014,31(3):339-370
[本文引用: 1]

[29] Patel M, Shabana AA . Locking alleviation in the large displacement analysis of beam elements: the strain split method
Acta Mechanica, 2018,229(7):2923-2946
[本文引用: 1]

[30] 刘铖, 田强, 胡海岩 . 基于绝对节点坐标的多柔体系统动力学高效计算方法
力学学报, 2010,42(6):1197-1205.
[本文引用: 1]

( Liu Cheng, Tian Qiang, Hu Haiyan . Efficient computational method for dynamics of flexible multibody systems
[本文引用: 1]

based on absolute nodal coordinate. Chinese Journal of Theoretical and Applied Mechanics, 2010,42(6):1197-1205 (in Chinese))
[本文引用: 1]

[31] Liu C, Tian Q, Hu HY . New spatial curved beam and cylindrical shell elements of gradient-deficient absolute nodal coordinate formulation
Nonlinear Dynamics, 2012,70(3):1903-1918
[本文引用: 1]

[32] Zhao ZH, Ren GX . A quaternion-based formulation of Euler--Bernoulli beam without singularity
Nonlinear Dynamics, 2012,67(3):1825-1835
[本文引用: 1]

[33] Bergou M, Wardetzky M, Robinson S , et al. Discrete elastic rods
ACM Transactions on Graphics, 2008,27(3):459-470
[本文引用: 1]

[34] Liu C, Tian Q, Hu HY . Dynamics of a large scale rigid--flexible multibody system composed of composite laminated plates
Multibody System Dynamics, 2011,26(3):283-305
[本文引用: 1]

[35] Liu C, Tian Q, Yan D , et al. Dynamic analysis of membrane systems undergoing overall motions, large deformations and wrinkles via thin shell elements of ANCF
Computer Methods in Applied Mechanics and Engineering, 2013,258:81-95
[本文引用: 1]

[36] Hou YS, Liu C, Hu HY . Component-level proper orthogonal decomposition for flexible multibody systems
Computer Methods in Applied Mechanics and Engineering, 2020,361:112691
[本文引用: 1]

[37] Simo JC, Tarnow N, Doblare M . Nonlinear dynamic of three-dimensional rods: Exact energy and momentum conserving algorithm
International Journal for Numerical Methods in Engineering, 1995,38(9):1431-1473
[本文引用: 1]

[38] Romero I, Armero F . An objective finite element approximation of the kinematics of geometrically exact rods and its use in the formulation of an energy--momentum conserving scheme in dynamics
International Journal for Numerical Methods in Engineering, 2002,54:1683-1716
[本文引用: 1]

[39] Marsden JE, West M . Discrete mechanics and variational integrators
Acta Numerica, 2001: 357-514
[本文引用: 1]

[40] Betsch P, Hesch C, Sanger N , et al. Variational integrators and energy-momentum schemes for flexible multibody dynamics
Journal of Computational and Nonlinear Dynamics, 2010,5(3):031001
[本文引用: 1]

[41] Schwarze M, Reese S . A reduced integration solid-shell finite element based on the EAS and the ANS concept-Large deformation problems
International Journal for Numerical Methods in Engineering, 2011,85:289-329
[本文引用: 1]

[42] Luo K, Liu C, Tian Q , et al. An efficient model reduction method for buckling analyses of thin shells based on IGA
Computer Methods in Applied Mechanics and Engineering, 2016,309:243-268
[本文引用: 1]

[43] Sun DW, Liu C, Hu HY . Dynamic computation of 2D segment-to-segment frictionless contact for a flexible multibody system subject to large deformation
Mechanism and Machine Theory, 2019,140:350-376
[本文引用: 1]

[44] Heyn T . On the modeling, simulation, and visualization of many-body dynamics problems with friction and contact. [PhD Thesis]. University of
Wisconsin-Madison, 2013
[本文引用: 1]

[45] Lew AJ, Marsden JE, Ortiz M , et al. Asynchronous variational integrators
Archive for Rational Mechanics and Analysis, 2003,167(2):85-146
[本文引用: 1]

[46] 许焕宾, 李伟杰, 史文华 等. 一种可重复展收的桁架结构及其胞元, 中国专利, 公开号: CN105923170B
2016 -04-26
[本文引用: 3]

[47] Farhat C, Lesoinne M, Letallec P , et al. FETI-DP: a dual--primal unified FETI method-part I: A faster alternative to the two-level FETI method
International Journal for Numerical Methods in Engineering, 2001,50(7):1523-1544
[本文引用: 1]

[48] 刘铖, 叶子晟, 胡海岩 . 基于区域分解的柔性多体系统高效并行算法
中国科学: 物理学力学天文学, 2017,47(10):12-22.
[本文引用: 1]

( Liu Cheng, Ye Zisheng, Hu Haiyan . An efficient parallel algorithm for flexible multibody systems based on domain decomposition method
Scientia Sinica Physica, Mechanica and Astronomica, 2017,47(10):12-22 (in Chinese))
[本文引用: 1]

[49] 李培, 马沁巍, 宋燕平 等. 大型空间环形桁架天线反射器展开动力学模拟与实验研究
中国科学: 物理学力学天文学, 2017,47(10):3-11.
[本文引用: 3]

( Li Pei, Ma Qinwei, Song Yanping , et al. Deployment dynamics simulation and ground test of a space hoop trus antenna reflector
Scientia Sinica Physica, Mechanica and Astronomica, 2017,47(10):3-11 (in Chinese))
[本文引用: 3]

[50] Li P, Liu C, Tian Q , et al. Dynamics of a deployable mesh reflector of satellite antenna: Parallel computation and deployment simulation
Journal of Computational and Nonlinear Dynamics, 2016,11(6):061005
[本文引用: 1]