两类基于局部标架梁单元的闭锁缓解方法1)

汤惠颖 张志娟 刘 铖,2) 刘绍奎

∗(北京理工大学宇航学院,北京 100081)

†(北京空间飞行器总体设计部,北京 100094)

引言

几何精确方法[1](geometrically exact formulation,GEF)与绝对节点坐标方法[2](absolute nodal coordinate formulation,ANCF)是描述大转动、大变形柔性多体系统的两类常用的建模方法.20 世纪80 年中期,李群/李代数基本理论逐渐与非线性有限元方法融合,以解决梁、板/壳结构的大转动、大变形动力学问题.几何精确建模方法在有限元理论框架下利用节点位置矢量与转动伪矢量作为节点参数描述三维梁大转动、大变形运动,对位置矢量与转动伪矢量独立插值[3].其中涉及的两个核心算法在于:第一、刚体旋转矩阵在李群SO(3)内采用乘法更新,对应的转动伪矢量由对数映射确定,满足刚体转动合成的几何意义[4].第二、通过转动参数定义梁截面局部标架,并在局部标架下研究梁微元的动能与弹性势能,即拉伸、剪切、扭转、弯曲应变以及截面角速度均定义在局部标架内[5-6].在多体系统动力学领域,Shabana[7]于1996 年提出了绝对节点坐标方法,其本质上也属于一种非线性有限元方法[8].该方法在处理物体转动问题上,摒弃复杂的转动参数,采用节点位置矢量以及沿物质坐标的斜率矢量来描述物体转动,可方便地建立大转动、大变形柔性多体系统动力学模型[9-10].

已有研究表明,在局部标架思想下,可构造出一类SE(3)几何精确梁单元[11]以及绝对节点坐标梁单元.局部标架梁单元能够规避刚体运动带来的几何非线性问题,离散数值模型中广义质量矩阵与切线刚度矩阵满足刚体变换的不变性,可明显地提高柔性多体系统动力学问题的计算效率.在有限元方法中,梁、板/壳单元普遍存在剪切、薄膜以及泊松闭锁问题.多年来,有限元领域学者提出了许多单元技术来有效地缓解闭锁.本文主要关注如何提高局部标架梁单元的收敛性,关于李群局部标架建模与计算方法涉及的其它细节本文不再赘述,可参见作者的相关研究工作.

缩减积分是解决闭锁问题最广泛使用的技术之一.Zienkiewicz 等[12]最早将缩减积分应用于壳单元、等参四边形单元和梁单元[13],提高了单元精度.Hughes 等[14]将选择缩减积分应用于板单元,缓解了薄板的剪切闭锁.Simo 等[15]采用缩减积分技术缓解了几何精确梁中的剪切闭锁.Malkus 和Hughes[16]证明了某些混合列式与使用缩减积分的位移公式的等价性,并且表明缩减积分单元可以达到多变量有限元的性能.Noor 和Peters[17]在曲梁中应用选择缩减积分,并通过比较刚度矩阵讨论了混合列式和位移列式的等价性和“近等价性”.

混合法不仅假设位移场,同时假设应力场或应变场,是一种高性能的多变量有限元方法,已被广泛用于解决体积闭锁、剪切闭锁和薄膜闭锁等问题.1985 年,Pian[18]提出了一种基于Hellinger-Reissner两场变分原理的混合单元,采用假设应力和位移场改善了梁和板的弯曲性能.1988 年,Liu 等[19]提出了基于Hu-Washizu 三场变分原理的弯曲超收敛单元,该单元在粗网格中显示出良好的精度.1994 年,Dorfi和Busby[20]提出了一种基于Hellinger-Reissner 两场变分原理的混合曲梁单元,该单元的位移和应力具有良好的收敛性.1990 年,Simo 和Rifai[21]提出了增强假设应变法,该方法是一种具有变分基础的闭锁缓解技术.随后,Simo 和Armero[22]将增强应变法推广到了几何非线性问题中.1993 年,Andelfinger 和Ramm[23]利用增强应变法开发了二维和三维板和壳单元,并证明了该方法与基于Hellinger-Reissner 两场变分原理的假设应力法[24]是等价的.

ANCF 全参数(fully parameterized)梁单元基于三维连续介质力学方法计算广义弹性力,当泊松比不为零时,其计算结果不能收敛到正确解.该问题即为泊松闭锁[25].起初,学者们令泊松比为零[26]或通过简化应力张量[27]除去泊松效应,缓解了泊松闭锁.随后,Gerstmayr 等[28]采用选择缩减积分技术缓解了泊松闭锁,同时在单元的变形模式中保留了泊松效应.2010 年,Matikainen 等[29]提出了一种高阶三维ANCF 梁单元,该单元对截面进行二次插值,有效地缓解了泊松闭锁,但是单元自由度过多.2018 年,Patel 和Shabana[30]提出了一种新型闭锁缓解技术——应变分解法(strain split method,SSM).该方法通过分解Green-Lagrange 应变张量并修改本构模型,假设仅与梁中线变形相关的低阶应变具有泊松效应,忽略正应变的高阶项的泊松效应,显著改善了梁的弯曲性能.

本文研究局部标架下几类梁单元的闭锁缓解方法.采用Hu-Washizu 三场变分原理缓解局部标架下李群SE(3)几何精确梁的剪切闭锁,采用应变分解法缓解局部标架下ANCF 全参数梁的泊松闭锁.通过对比几何精确与绝对节点坐标两类建模方法,分析基于局部标架的李群SE(3)几何精确梁在计算精度与效率方面的优势.本文中的李群SE(3)几何精确梁与ANCF 梁单元均基于局部标架,为了方便起见,下文省去“基于局部标架”几个字.

1 基于李群SE(3)局部标架的几何精确梁理论

1.1 梁的运动学描述及本构模型

基于经典的Timoshenko 梁假设,李群几何精确梁理论融合非线性有限元方法与几何力学思想,在李群SE(3)内建立系统的平衡方程,可精确地描述柔性梁的大转动、大变形刚柔耦合动力学特性,同时可规避刚体运动带来的几何非线性[31].

梁的初始构型和当前构型如图1 所示.引入参考坐标系{O;e1,e2,e3}.在当前构型中,在梁中线的每一点定义随体基矢量(t,n,b),引入梁中线位置矢量x和截面旋转矩阵R,初始构型对应的量为(·)0.

图1 几何精确梁初始构型与当前构型Fig.1 Initial and current configurations of the geometrically exact beam

在梁的当前构型中,截面上任意一点p的位置矢量为

式中,弹性系数矩阵D=diag(EA,ks2GA,ks3GA,GJs,EI2,EI3),其中E和G分别为杨氏模量和剪切模量,A,I2和I3分别表示梁初始构型的截面面积和截面惯性矩,Js表示圣维南扭转常数,ks2和ks3为截面的剪切修正系数.

1.2 SE3 几何精确梁控制方程及线性化

系统的动力学平衡方程为

式中,δWint为内力所做的虚功,δWext为外载荷所做的虚功,δWine为惯性力所做的虚功.

为获得一个对称的切线刚度矩阵,本文在前一个收敛步所在切平面进行物理量的一次与二次变分.因此,应变的变分可以表示为

式中,T 为SE(3)群的切空间算符,δh=[δdTδηT]T,δd和δη分别为局部标架下的虚位移和虚转角,B 为应变算符.内部虚功可以表示为

外部虚功可以表示为

式中,pext为作用于p点的外载荷.惯性虚功可以表示为

式中,v=[UTωT]T,U和ω分别为局部标架下的线速度与角速度,顶标·表示变量对时间t的全导数,ρA为单位长度的材料密度,I3表示3×3 的单位矩阵,J为截面形心惯性矩阵.

对内部虚功做线性化可得

应变算符B 的线性化的转置为

为了进一步简化∆(δWint)的几何部分,引入两个矩阵ΞT和ΞT′,它们的定义如下

式中,A为一个任意的6×1 的矢量,记号和分别是矩阵ΞT′的前半部分和后半部分,其具体形式可参见文献[11].因此,几何部分为

最后,内部虚功的线性化可以写为

需要指出的是,由于δh为局部标架下的平动与转动虚位移,∆h为局部标架下的平动与转动无限小位移,它们均与该点的任意刚体位移无关.因此,单元的弹性力及其切线刚度矩阵均满足刚体变换的不变性.对于广义惯性力及其广义质量矩阵也同时具有类似的形式.

1.3 SE3 几何精确梁有限元离散格式

对梁中线位置矢量x和截面相对旋转伪矢量Θr[11]进行空间离散

式中,NI是节点I的形函数,xI是梁中线上节点I的位置矢量,Rr是参考点所在截面的旋转矩阵,RI是节点I所在截面的旋转矩阵,旋转伪矢量为转轴单位矢量与转角的乘积.

离散应变算符B 得到应变−位移矩阵BI,它的表达式为

弹性力可以表示为

材料刚度和几何刚度矩阵可以表示为

对于非常细长的梁结构,采用一阶插值描述几何精确梁弯曲时将会产生剪切闭锁.由于梁非常细长,截面内横向剪切应变处处接近于0.而由于平动参数x和转动参数θ采用了相同的插值函数,使式(5)中横向剪切应变γs中的R(θ)和x′两项具有不同阶次,γs=0 不可能处处满足,导致能量泛函中剪切变形能项的量级不正确,由此带来了剪切闭锁.此时,单元表现得十分刚硬,即计算出来的位移值远小于参考解,总体刚度矩阵很大.

2 基于局部标架的绝对节点坐标方法的三维全参数梁理论

本节简要介绍基于局部标架的绝对节点坐标全参数梁单元建模方法.某质点从初始位置X运动到当前位置r,变形梯度可以表示为

Green-Lagrange 应变张量可以表示为

应变能可以表示为

式中,D为弹性矩阵,V为梁的体积.应变能对广义坐标e求偏导得弹性力为

弹性力对广义坐标求偏导得刚度阵为

相比于几何精确梁单元,ANCF 单元没有转动参数,不能直接定义局部标架,需要通过平动位移场定义转动场.此处采用ANCF 斜率矢量定义该点的局部标架基矢量

式中,rX和rY为当前位置矢量r对物质坐标的偏导数.由于平动位移场在局部标架中描述,单元的刚体运动可自然被消除.因此,系统质量矩阵以及切线刚度矩阵满足刚体运动的不变性,从而可极大地减少系统质量矩阵与刚度矩阵在仿真过程中的更新次数.

ANCF 全参数梁单元会遇到泊松闭锁问题,闭锁产生的原因如下.由于ANCF 全参数梁是连续介质单元,梁的中线为三阶Hermite 插值,在泊松效应的作用下,截面应该变形为曲面.然而横截面为一阶插值,不能描述截面的曲面变形模式,因此导致泊松闭锁.

3 两类局部标架梁单元的闭锁缓解方法

3.1 基于Hu-Washizu 三场变分原理的剪切闭锁缓解方法

Hu-Washizu 三场变分原理不仅假设位移场,同时假设应力和应变场,且三类场变量的插值相互独立.假设剪切应变和假设剪切应力的分布为

式中,Ss=[0,,0,0,0]T.剪切部分的弹性力可以表示为

剪切部分的刚度矩阵可以表示为

式中,

离散方程组可以表示为

式中,Km和Kb分别为薄膜和弯曲部分的刚度矩阵,∆hN为节点运动矢量的增量,Res为残差.将上式静力缩聚化为单场形式

一阶插值的几何精确梁会产生严重的剪切闭锁,二阶及二阶以上插值的几何精确梁剪切闭锁较轻.因此,本文关于几何精确梁剪切闭锁处理方法着重于一阶梁单元.Hu-Washizu 三场变分原理中,位移场、应变场和应力场均独立插值.式(28)假设横向剪切应变为常数,可避免产生过大的剪切应变能,有效缓解了剪切闭锁.

此外,对于梁单元,缩减积分技术也是缓解剪切闭锁的有效手段.数值算例展示了采用Hu-Washizu三场变分原理与缩减积分技术缓解梁单元剪切闭锁的效果.相比较而言,缩减积分技术由于其简便性在梁单元中应用较为广泛,但用在板壳单元中会产生沙漏变形模式.进一步,上述Hu-Washizu 三场变分原理的研究将为低阶几何精确板壳单元的闭锁缓解起到一定的借鉴作用.

3.2 基于应变分解法的泊松闭锁缓解方法

ANCF 全参数梁单元的位置场和变形梯度可以写为

εc中只有低阶应变,且εc仅与梁中线参数X相关;εk中包含高阶应变,这些高阶应变会导致截面变形和弯曲变形.为了缓解闭锁,泊松耦合只用在正应变的低阶项中,即仅在εc中考虑泊松效应,忽略正应变的高阶项的泊松效应.因此,梁截面不会变形为曲面,缓解了泊松闭锁.第二类Piola-Kirchhoff应力的Voigt形式可以写为

式中,应变也是Voigt 形式.基于平面应变假设,弹性系数矩阵可以被定义为

式中,λ=Eν/[(1+ν)(1 −2ν)]和µ=E/[2(1+ν)]为拉梅常数,E和ν 分别为杨氏模量和泊松比.应变能可以表示为

弹性力可以表示为

材料和几何刚度可以表示为

4 三类常用的有限元插值方法

本节简要介绍三类有限元插值方法,包括Lagrange 插值、Hermite 插值以及非均匀有理B-样条(NURBS).其中,Lagrange 插值在有限元中应用最为广泛;绝对节点坐标方法则采用Hermite 插值,以保证位移场的C1连续;NURBS 的基函数可以构造任意阶连续的近似函数,建立了计算机辅助设计(CAD)与计算机辅助工程(CAE)间的桥梁,该方法被称为等几何分析.

为对比分析不同插值方法的优劣势,本文涉及一阶与三阶Lagrange 插值、三阶Hermite 插值以及三阶NURBS 插值.

4.1 Lagrange 插值

给定函数f(x)在n+1 个互不相同的点xi(i=0,1,2,···,n)上的函数值f(xi)=yi,n次Lagrange 插值多项式可以表示为

4.2 Hermite 插值

给定函数f(x)在n+1 个互不相同的点xi(i=0,1,2,···,n)上的函数值yi及一阶导数值,Hermite插值多项式可以表示为

Hermite 插值基函数Ai(x)和Bi(x)为

4.3 非均匀有理B-样条

在区间[a,b]内,设U={u0,u1,u2,···,um}是一不减的实数序列.以U作为样条结点,p次的第i个B-样条基函数Ni,p(u)定义为

式中,Pi是控制点,ωi是权系数.

5 直/曲梁结构大变形静力学分析

本节考察XOY平面内的直梁与曲梁模型,验证上述闭锁处理方法的有效性.其中,悬臂直梁长L=1 m;悬臂曲梁半径R=1 m,初始构型为1/4圆.梁的截面均为矩形,宽和高为0.01 m;材料的杨氏模量为E=2.1×1011Pa,泊松比为ν=0.3.首先,在直梁和曲梁末端分别施加集中力(30,20,400)N (直梁)、(30,20,100)N (曲梁),用ABAQUS 计算得到梁端点沿z方向的位移分别为0.5143 m (直梁)、0.6138 m(曲梁).计算的过程中,上述两类载荷均分为10 个加载步均匀加载.此外,对直梁与曲梁进行纯弯曲测试.在直梁和曲梁末端均施加z方向的力矩Mz=−πEI/2L,由材料力学解答知,梁端沿z方向转角位移的解析解为θz=−π/2.

5.1 悬臂直梁和曲梁低阶单元闭锁缓解方法数值测试

本节考查SE(3)几何精确一阶直梁与曲梁单元在三维复杂应力状态下,缩减积分技术和Hu-Washizu三场变分原理缓解剪切闭锁的能力.

悬臂直梁和曲梁末端受集中力的作用,变形如图2 所示.单元位移场采用一阶Lagrange 插值.对比分析完全积分(exact integration,EI)、缩减积分(reduced integration,RI)和Hu-Washizu 三场变分原理(Hu-Washizu variational principle,HWVP)三种情况,在不同自由度(degrees of freedom,Dofs)下悬臂梁末端点z方向的位移uz,计算相对误差,将结果列于表1和表2 中.同时,采用ABAQUS 软件中B31 梁单元(该单元为三维二节点一阶梁单元,每个节点6 个广义坐标,分别是3 个平移和3 个转动参数)对本算例进行仿真计算.

图2 (a)直梁和(b)曲梁变形前后示意图Fig.2 Initial and deformed configurations of(a)the straight beam and(b)the curved beam

表1 悬臂直梁末端z 方向位移和相对误差Table 1 End displacements of z direction and relative error of the straight cantilever beam

表2 悬臂曲梁末端z 方向位移和相对误差Table 2 End displacements of z direction and relative error of the curved cantilever beam

表1 和表2 的数值结果表明:采用缩减积分技术的SE(3)几何精确梁与ABAQUS 中B31 单元收敛结果基本一致.采用Hu-Washizu 三场变分原理缓解闭锁的SE(3)几何精确梁与其他两类方法的收敛结果不同,其中,直梁模型收敛值偏大,曲梁模型收敛值偏小.在三维复杂应力状态下,一阶直梁和曲梁均有剪切闭锁问题.缓解闭锁后,直梁的计算精度提高了98%以上,曲梁的计算精度提高了52%以上,效果十分显著.对于一阶直梁,以上两种缓解闭锁的方法计算精度无明显差别.对于一阶曲梁,缩减积分技术的计算精度略高于Hu-Washizu 三场变分原理.

5.2 悬臂直梁和曲梁纯弯曲测试

本节考查SE(3)几何精确一阶直梁与曲梁单元在纯弯曲应力状态下,缩减积分技术和Hu-Washizu三场变分原理缓解剪切闭锁的能力,同时评价单元描述纯弯曲的能力.

在悬臂直梁和曲梁末端施加集中弯矩,此时,悬臂梁发生纯弯曲,处于简单应力状态.单元位移场采用一阶Lagrange 插值.对比分析完全积分、缩减积分和Hu-Washizu 三场变分原理三种情况下,悬臂梁末端点z方向的转角位移θz,计算相对误差,将结果列于表3 和表4 中.同时,采用ABAQUS 软件中B31 梁单元对本算例进行仿真计算.

表3 悬臂直梁末端z 方向角位移和相对误差Table 3 Angular end displacements of z direction and relative error of the straight cantilever beam

表4 悬臂曲梁末端z 方向角位移和相对误差Table 4 Angular displacements of z direction and relative error of the curved cantilever beam

表3 和表4 的数值结果表明:在纯弯曲应力状态下,一阶直梁和曲梁均有剪切闭锁问题.缓解闭锁后,直梁和曲梁模型的计算精度均提高了98%以上,效果十分显著.而且,以上两种缓解闭锁的方法和ABAQUS 中B31 单元的计算精度完全一致.表3中,一阶直梁在任意单元数目下,转角与解析解的相对误差均为零.这是由于SE(3)几何精确直梁模型能够精确地构造一个常弯曲单元,因此仅用1 个单元就能够精确模拟纯弯曲状态下的转动场.表4 中,一阶曲梁模型的转角渐近收敛,这是由于一阶插值函数无法精确地描述一段圆弧.

同时,采用三阶Lagrange 插值构造了高阶几何精确直梁与曲梁单元.通过上述悬臂直梁与曲梁算例,测试了单元收敛精度,其中同样采用Hu-Washizu三场变分原理处理剪切闭锁.数值结果表明:处理闭锁不会提升单元收敛性能,进一步说明了高阶直梁与曲梁单元不存在剪切闭锁.在5.3 节中,将进一步对比处理闭锁后的一阶单元与高阶单元的收敛性.

5.3 几类局部标架梁单元收敛性对比分析

本节对比一阶和三阶SE(3)几何精确梁单元、一阶R3×R3几何精确梁、ANCF 缩减(slope deficiency)梁、ANCF 全参数梁与ABAQUS 软件中B31 梁单元的收敛性.其中,一阶几何精确梁与ANCF 全参数梁单元均已处理闭锁.

在悬臂直梁和曲梁末端施加集中力.计算以上几类建模方法在不同自由度下悬臂梁末端点z方向的位移uz和相对误差,将结果列于表5 和表6 中.ANCF 缩减曲梁单元无算进行本算例的计算。

表5 悬臂直梁末端z 方向位移和相对误差Table 5 End displacements of z direction and relative error of straight cantilever beams

表6 悬臂曲梁末端z 方向位移和相对误差Table 6 End displacements of z direction and relative error of curved cantilever beams

表5 和表6 的数值结果表明:直梁模型的五种建模方法收敛结果基本一致.对于曲梁模型,SE(3)几何精确梁、R3×R3几何精确梁和ABAQUS 中B31 单元收敛结果基本一致,ANCF 全参数梁单元收敛结果略大于其他三类方法的收敛值.需要指出的是,表5 中三阶NURBS 插值的SE(3)几何精确梁单元在24 个自由度时,虽然计算结果的相对误差为6.6×10−4,但此时单元还没有收敛,60 个自由度时,单元才收敛.对于高阶单元,三阶Lagrange 比三阶NURBS 插值的SE(3)几何精确梁的计算精度高,几类几何精确梁均比ANCF 梁单元的计算精度高.ANCF 缩减梁比ANCF 全参数梁单元的计算精度高.

对于一阶直梁单元,两类几何精确梁具有完全相同的计算精度.对于一阶曲梁单元,R3×R3几何精确梁的计算精度略高于SE(3)几何精确梁.对于一阶直梁和曲梁单元,两类几何精确方法均比ABAQUS中B31 单元的计算精度高.显然高阶单元比低阶单元精度高,但是高阶单元计算效率低,即自由度相同时,高阶单元刚度矩阵内非零元素的个数多.

6 高转速/大变形柔性多体系统动力学验证算例

本节通过高转速曲柄滑块机构以及空间双摆的动力学仿真,分析局部标架下的梁单元在动力学仿真中的计算精度差异,以及消除刚体运动带来的几何非线性后梁单元的数值特性.由于商业软件无法直接进行大变形柔性多体系统动力学仿真,因此本节不再采用商业软件进行仿真对比.为了设计特定的小变形或大变形的柔性多体系统动力学模型,本节算例中部分材料参数未与真实材料对应.

本节数值算例共涉及如下3 类建模方法:(1)一阶和三阶Lagrange 插值的SE(3)几何精确梁单元;(2)一阶Lagrange 插值的R3×R3几何精确梁单元;(3)ANCF 全参数梁单元.需要指出的是,除了三阶SE(3)几何精确梁单元,其他单元均已处理闭锁.

6.1 发生小变形和中等变形的高速旋转的平面曲柄滑块机构

本节对一个典型的多体系统——曲柄滑块机构进行动力学仿真.如图3 所示,平面曲柄滑块机构由两根杆与一个滑块组成,杆I 一端与地面采用球铰相连接,另一端与杆II 一端也通过球铰连接,杆II 末端与一个仅能在X轴上运动的刚体滑块铰接.

图3 曲柄滑块机构Fig.3 A crank-slider mechanism

本算例中,通过杆II 中点C到其首尾连线的距离d度量杆II 的弹性变形,变形率为d与杆原长之间的比值.同时,以1944 个自由度的R3空间ANCF全参数梁单元所得到的数值结果作为参考值,建模方法数值结果收敛标准设为变形场最大值与参考解的相对误差在1%以内.

6.1.1 模型I

杆I 与杆II 的长度分别为0.3 m 和0.5 m,两杆的截面均为矩形,宽与高为0.05 m,材料参数设为:ρ=2000 kg/m3,E=8.2×1012Pa,ν=0.3.杆II 末端滑块的质量为0.1 kg.系统初始静止放置于XOY水平面内,杆I 与X轴夹角为0◦.在曲柄上施加Z方向力矩900 N·m,持续时间为0.2 s,杆I 的最大转速约为1033 rad/s.时间积分算法采用广义α 方法,仿真时间为0.2 s,时间步长设为1.0×10−4s,算法参数谱半径设为0.8.

通过试算,三阶SE(3)几何精确梁84 个自由度收敛,采用缩减积分技术和Hu-Washizu 三场变分原理的一阶SE(3)几何精确梁与一阶R3×R3几何精确梁均在156 个自由度收敛,ANCF 全参数梁312 个自由度收敛.图4 给出了上述三类建模方法收敛时杆II 变形率随时间的变化曲线.随着转速增大,杆II 的变形逐渐增大,最大变形能够达到杆长度的0.04%.

6.1.2 模型II

杆I 与杆II 的长度分别为0.3 m 和0.5 m,两杆的截面均为矩形,宽与高为0.03 m,材料参数设为:ρ=5600 kg/m3,E=4.1×1011Pa,ν=0.3.杆II 末端滑块质量为0.1 kg.系统初始静止放置于XOY水平面内,杆I 与X轴夹角为0◦.在曲柄上施加Z方向力矩500 N·m,持续时间为0.2 s,杆I 的最大转速约为518 rad/s.时间积分算法采用广义α 方法,仿真时间为0.2 s,时间步长设为3.0×10−5s,算法参数谱半径设为0.8.

通过试算,三阶SE(3)几何精确梁84 个自由度收敛,采用缩减积分技术和Hu-Washizu 三场变分原理的一阶SE(3)几何精确梁与一阶R3×R3几何精确梁均在300 个自由度收敛,ANCF 全参数梁504 个自由度收敛.图5 给出了上述三类建模方法收敛时杆II 变形率随时间的变化曲线.随着转速增大,杆II 的变形逐渐增大,最大变形能够达到杆长度的4%.

图5 杆II 变形率随时间变化曲线Fig.5 Deformation rate of the rod II

由图4 和图5 可知,SE(3)几何精确梁、R3×R3几何精确梁和ANCF 全参数梁均能达到收敛的数值结果,验证了这三类建模方法在描述高转速、小变形以及中等变形多体系统动力学问题时的正确性.采用缩减积分技术与Hu-Washizu 三场变分原理可以缓解动力学问题中的剪切闭锁,两种方法的收敛速度无明显差别.一阶R3×R3几何精确梁与一阶SE(3)几何精确梁收敛速度无明显差别.

在一个牛顿迭代步中,最耗时的两个部分包括计算系统刚度矩阵与迭代矩阵的LU 分解/回代.由于本文中的梁单元均采用局部标架建模方法,刚度矩阵更新次数大幅降低.因此,可通过计算迭代矩阵中非零元素的个数来评价一个牛顿迭代步的计算速度.对于模型I,几类一阶几何精确梁迭代矩阵非零元素的个数均为2808,三阶SE(3)几何精确梁和ANCF全参数梁迭代矩阵非零元素的个数分别为3528 和11 232.对于模型II,几类一阶几何精确梁迭代矩阵非零元素的个数均为5400,三阶SE(3)几何精确梁和ANCF 全参数梁迭代矩阵非零元素的个数分别为3528 和18 144.由以上数据可知,在达到相同的计算精度时,几类一阶几何精确梁计算效率相同,SE(3)和R3×R3几何精确梁的计算效率明显高于ANCF 全参数梁.

其次,通过分析系统质量矩阵和刚度矩阵更新情况,验证基于局部标架的SE(3)几何精确方法能够有效地减轻刚体运动带来的非线性问题.若忽略外力对应的Jacobian 矩阵,系统迭代矩阵的表达式为

式中,h为积分步长,β 为广义α 方法算法参数,M为广义质量矩阵,K为切线刚度矩阵,Φq为约束方程的Jacobian 矩阵.本算例中,牛顿迭代收敛标准设为广义坐标修正量||∆||2<1.0×10−6.表7 给出了一阶SE(3)几何精确梁迭代矩阵的主要部分Jq=h2β(M+K)的更新次数(number of update,NOU)以及仿真过程中牛顿迭代总次数(number of iterations,NOI).表7 中“NNI ≥i”表示一个时间步内,牛顿迭代次数大于等于i步时,Jq强制更新.

表7 不同模型Jq 的更新次数及仿真过程牛顿迭代总次数Table 7 Number of updating Jq and number of Newton iterations about different models

模型I 的数值结果表明:对于基于局部标架的SE(3)几何精确梁模型,若每个牛顿迭代步中Jq均更新,则整个仿真过程共需要进行4713 次牛顿迭代步,同时Jq也需要更新4712 次.当NNI ≥3 时(即每个时间步中牛顿迭代达到三步时更新一次Jq),Jq需要更新357 次,总迭代次数增加了30%左右.当NNI≥4 时(即每个时间步中牛顿迭代达到四步时更新一次Jq),Jq需要更新10 次,总迭代次数增加了37%左右.模型II 与模型I 情况类似,这里不再赘述.由此可得,对于此类高速转动的多体动力学问题,基于局部标架的SE(3)几何精确梁能够极大地降低刚体运动带来的几何非线性.

6.2 发生中等变形和超大变形的空间双摆

本节对空间双摆模型进行动力学仿真.如图6 所示,空间双摆模型初始放置于XOY水平面上,并在重力作用下开始运动.与算例6.1 一致,本算例通过杆II 中点C到其首尾连线的距离d度量杆II 的弹性变形,变形率为d与杆原长之间的比值.同时,以1944个自由度的R3空间ANCF 全参数梁单元所得到的数值结果作为参考值,建模方法数值结果收敛标准设为变形场最大值与参考解的相对误差在1%以内.

图6 空间双摆Fig.6 A spatial double pendulum

6.2.1 模型III

双摆两个摆臂的长度分别为0.3 m 和0.5 m,两个摆臂的截面均为矩形,宽与高为0.5 mm,材料参数设为:ρ=2700 kg/m3,E=2.1×1010Pa,ν=0.3,重力加速度为g=9.81 m/s2.时间积分算法采用广义α 方法,仿真时间1 s,时间步长设为5.0×10−4s,算法参数谱半径设为0.8.

通过试算,三阶SE(3)几何精确梁84 个自由度收敛,采用缩减积分技术和Hu-Washizu 三场变分原理的一阶SE(3)几何精确梁均在300 个自由度收敛,ANCF 全参数梁1464 个自由度收敛.300 个自由度的一阶R3×R3几何精确梁在t=0.457 5 s 发散.图7给出了上述三类建模方法收敛时杆II 变形率随时间的变化曲线,最大变形能够达到杆长度的5.5%.

图7 杆II 变形率随时间变化曲线Fig.7 Deformation rate of the rod II

6.2.2 模型IV

双摆两个摆臂的长度分别为0.3 m 和0.5 m,两个摆臂的截面均为矩形,宽与高为0.3 mm,材料参数设为:ρ=2700 kg/m3,E=1.5×1010Pa,ν=0.3,重力加速度为g=9.81 m/s2.时间积分算法采用广义α 方法,仿真时间0.5 s,时间步长设为5.0×10−4s,算法参数谱半径设为0.8.

通过试算,三阶SE(3)几何精确梁84 个自由度收敛,采用缩减积分技术和Hu-Washizu 三场变分原理的一阶SE(3)几何精确梁156 个自由度收敛,ANCF 全参数梁1752 个自由度收敛.156 个自由度的一阶R3×R3几何精确梁在t=0.444 5 s 发散.图8给出了上述三类建模方法收敛时杆II 变形率随时间的变化曲线,最大变形能够达到杆长度的15%.

由图7 和图8 可知,SE(3)几何精确梁、R3×R3几何精确梁和ANCF 全参数梁均能达到收敛的数值结果,验证了这三类建模方法在描述中等变形和超大变形多体系统动力学问题时的正确性.采用缩减积分技术与Hu-Washizu 三场变分原理均可有效缓解动力学中的剪切闭锁问题,两种方法的收敛速度无明显差别.

图8 杆II 变形率随时间变化曲线Fig.8 Deformation rate of the rod II

下面通过计算迭代矩阵中非零元素的个数来评价一个牛顿迭代步的计算速度.对于模型III,几类一阶几何精确梁迭代矩阵非零元素的个数均为5400,三阶SE(3)几何精确梁和ANCF 全参数梁迭代矩阵非零元素的个数分别为3528 和52 704.对于模型IV,几类一阶几何精确梁迭代矩阵非零元素的个数均为2808,三阶SE(3)几何精确梁和ANCF 全参数梁迭代矩阵非零元素的个数分别为3528 和63 072.由以上数据可知,在达到相同的计算精度时,几类一阶几何精确梁计算效率相同,SE(3)几何精确梁的计算效率明显高于ANCF 全参数梁.

通过与6.1 节中等变形的平面多体系统对比,可以发现,随着柔性体变形的增大,SE(3)几何精确梁比ANCF 全参数梁单元优势更明显.一阶R3×R3几何精确梁在仿真计算的过程中发散,是由于本文仅采用Rescaling 技术[32]处理转动参数,没有进一步采用其他方法避免转动参数的奇异性.相比较而言,本文提出的局部标架建模方法则可自然地避免转动参数的奇异性问题.

与6.1 节统计方法类似,通过分析系统质量矩阵和刚度矩阵更新情况,以说明基于局部标架的SE(3)几何精确建模方法描述大转动、大变形非线性运动的能力.本算例中,牛顿迭代收敛标准设为广义坐标修正量||∆||2<1.0×10−6.表8 中给出了一阶SE(3)几何精确梁Jq的更新次数以及仿真过程中的牛顿迭代总次数.

模型IV 的数值结果表明:对于基于局部标架的SE(3)几何精确梁模型,若每个牛顿迭代步中Jq均更新,则整个仿真过程共需要进行2479 次牛顿迭代步,同时Jq也需要更新2478 次.当NNI ≥3 时(即每个时间步中牛顿迭代达到三步时更新一次Jq),Jq需更新440 次,总迭代次数增加了40%左右.当NNI ≥4 时(即每个时间步中牛顿迭代达到四步时更新一次Jq),Jq需要更新283 次,总迭代次数增加了54%左右.模型III 与模型IV 情况类似,这里不再赘述.由此可得,对于此类大变形多体动力学问题,基于局部标架的SE(3)几何精确梁仍然能够明显地降低刚体运动带来的几何非线性.

表8 不同模型Jq 的更新次数及仿真过程牛顿迭代总次数Table 8 Number of updating Jq and number of Newton iterations about different models

7 结论

本文基于李群SE(3)局部标架,研究几类梁单元的闭锁缓解方法.采用Hu-Washizu 三场变分原理缓解了SE(3)几何精确梁单元中的剪切闭锁,采用应变分解法缓解了ANCF 全参数梁单元中的泊松闭锁,并对SE(3)几何精确梁、R3×R3几何精确梁和ANCF全参数梁的单元性能进行了算例对比.静力学与动力学算例结果均表明,缓解闭锁后的SE(3)几何精确梁计算精度高.在进行大转动、大变形动力学计算时,R3×R3几何精确梁单元会出现转动参数奇异问题,若处理不当则会不收敛,而SE(3)几何精确梁单元可避免转动参数奇异性问题.SE(3)几何精确梁的广义质量矩阵主项为常数,与转动参数相关的切线刚度矩阵满足刚体转动的不变性,在计算中无须更新,可极大地减少迭代矩阵的更新次数,有效地提高计算效率.基于SE(3)群局部标架的几何精确建模方法更适合用于描述梁结构的大转动、大变形运动.