王 悦 崔雅琦 於祖庆 兰 朋,2) 陆念力

* (哈尔滨工业大学机电工程学院,哈尔滨 150001)

† (河海大学机电工程学院,江苏常州 213002)

引言

随着柔性轻质结构在航天、汽车以及船舶工程等领域日益广泛的应用,柔性多体系统的动力学分析也愈加重要.然而,以线弹性小变形假设为基础的,将节点位移和转动作为坐标的传统有限单元不适合解决存在大位移、大变形的柔性多体系统的动力学问题.此外,对于传统结构有限元,将几何数据转换为有限元网格数据是困难且耗时的[1].据估计,在航空航天、船舶制造和汽车工业中,大约80%的分析时间用于网格生成[2].在计算机辅助工程(CAE)领域,有限单元只是计算机辅助设计(CAD)几何模型的近似表示而非精确描述.这种CAE 网格与CAD 几何之间的区别将会导致精确性问题,例如在壳分析中,随着几何缺陷的增加,屈曲荷载有相当大的降低[2].

为了准确刻画柔性体大位移、大变形的动力学行为,研究者将弹性力学与有限元方法和多体动力学理论相结合,形成了绝对节点坐标方法(ANCF).ANCF 利用梯度代替转角对柔性体位移场进行插值并与多体系统动力学理论相结合,使得其具备描述柔性体大变形、大转动的能力,从而成为多柔体系统动力学领域的研究热点[3].此外,研究者又发现其与计算机辅助分析软件中的非均匀有理B 样条(NURBS)几何体家族之间存在线性转化关系,从而可以实现从几何模型到ANCF 单元网格之间的快速转换,避免了网格划分的步骤.相关的研究始于缆索单元与B 样条曲线之间的转化关系[4],并延伸至有理的ANCF 缆索单元与NURBS 曲线间关系[5].相关研究深刻揭示了二者在几何体形状描述、连续性控制等方面的高度相似性[6-7].此后,学者们将研究拓展至ANCF 双参数单元[8-9].研究结果表明,当单元维数升高,NURBS 几何体张量积形式的参数方程会引入冗余自由度,从而需要对ANCF 单元进行改造从而使其与之匹配.例如增加混合梯度向量,或在板单元中间添加额外节点等.这一问题在三维单元中更加严重.Yu 和Cui[10]提出了一种48 自由度的八节点实体梁单元,可以精确离散使用相同的基函数阶次Bézier 体所代表的几何模型.Ma 等[11]研究了一种基于三次有理Bézier 体积的三维有理绝对节点坐标公式(RANCF)流体单元,准确描述初始曲线形状的液柱.

对计算机辅助工程和计算机辅助设计进行整合的另外一个方向就是由Hughes 等[2]于2005年提出的等几何分析(IGA),其概念是将CAD 和有限元分析(FEA)两个领域统一起来,使用相同的基础进行几何描述和分析[12].IGA 方法因为具有实现无缝整合设计和分析的潜力[13]而受到了力学领域的关注[2].学者们进行了很多针对壳体和板的等几何分析[14-17],这些研究多是基于当前CAD 系统的行业标准NURBS.然而,由于NURBS 控制点呈矩形网格分布,从而导致其拓扑结构中存在大量冗余控制点[18],进而导致分析效率降低.此外,修剪后的NURBS 表面不可避免地存在间隙和重叠是影响CAD,CAM和CAA 系统互操作性的最严重的障碍之一[19].为了克服NURBS 的这些缺点,Sederberg 等[20]提出了T 样条.与NURBS 曲面相比,T 样条曲面的一个优点是允许局部细化.由于T 交叉点(T-junction)的存在使得T 样条允许控制点局部插入到控制网格中,而非整行或整列地增加控制点[20],从而大大减少了多余控制点的数量[18].T 样条的另一个优点是T 样条模型是水密的.多个NURBS 补丁可以合并成一个单一的水密T 样条,任何修剪过的NURBS 对象都可以转换为未修剪的T 样条[21].因此,T 样条被认为是未来CAD 行业的新标准,将在CAD 与CAE 的集成中发挥重要作用.Casquero 等[22]利用适合分析的任意度T 样条曲面进行了完全非线性薄壳的结构分析.Casquero 等[23]还利用适合分析的T 样条来解决Kirchhoff−Love 壳问题.Dimitri 等[24]提出了一种基于T 样条的等几何分析方法,应用于大变形条件下变形体间的无摩擦接触问题.

在实际应用中,由于ANCF 单元使用梯度作为节点向量,其建模过程较为复杂.另外,如果将一个薄板离散为几个ANCF 薄板单元,其单元边界上的梯度不连续,因此格林−拉格朗日应变也是不连续的.而使用IGA 方法可以直接使用CAD 软件对分析对象进行几何建模,而且几何模型的连续性条件可以自然得通过结点重复性保证.为此,本文将在等几何分析的框架下,开展基于T 样条曲面单元的基尔霍夫薄板动力学分析方法研究.本文的主要贡献在于:

(1)基于可局部细化的T 样条曲面,提出局部细化策略解决动力学分析中局部应变变化较大的问题,并根据T 样条曲面几何模型的局部细化算法创建变自由度系统动力学方程的求解算法;

(2)通过变网格柔性薄板与刚性球的碰撞问题,展示所提出局部细化T 样条曲面单元在接触碰撞问题中的应用价值.

1 基于T 样条的等几何基尔霍夫薄板

1.1 T 样条曲面

T 样条可视为控制网格上带有T 交叉点的NURBS 曲面[25],可以通过控制点与T 样条基函数Ti(u,v)相乘得到.T 样条曲面的定义式如下[26]

式中,ST为T 样条曲面上物质坐标为 (u,v)的点的全局坐标,Pi=(xi,yi,zi) 为第i个控制点的全局坐标,n为控制点的数量,Ti为第i个控制点所对应的有理T 样条基函数,其表达式为

式 中,混合函数Bi为B 样条基函数Ni(u)和Ri(v)的乘积. ωi是第i个控制点所对应的权重,W为权值与基函数的乘积之和.

与NURBS 不同的是,每一个T 样条基函数都是基于局部结点矢量而不是全局结点矢量来计算的.局部结点矢量的建立依赖于T 网格和锚点的建立,根据这两者可以得到T 样条曲面上每一个控制点所对应的局部结点矢量和T 网格和锚点的具体定义参见文献[26-27],并可依据参考文献[28]得到三次T 样条基函数以及其导数的表达式.

1.2 T 样条薄板的运动学描述

由于基尔霍夫薄板忽略横向剪切变形,法向量始终垂直于中性层,其运动学描述可以简化为其中性层的运动学描述.为此,首先需要采用T 样条单元对薄板中性层进行运动学描述.

T 样条单元是由T 样条基函数的简化连续线构成的物理区域[29].在等几何分析中,利用IEN 数组可以确定每个样条单元的连接关系.IEN 中的I 表示数组是整数值,EN 是“element nodes”的首字母缩写,表示“单元节点”,该数组可以将局部基函数号和单元号映射到相应的全局控制点号[29].对于第i个T 样条单元而言,其形函数以及单元节点向量ei可以表示为

式中,n是一个T 样条单元中包含的控制点的个数,IEN(j) 为该T 样条单元的IEN 数列中的第j项,PIEN(j)为第 I EN(j)个控制点的坐标,I为 3 ×3的单位矩阵,TIEN(j)(u,v) 为第IEN(j)个控制点所对应的有理T 样条基函数.

综上所述,基尔霍夫薄板的运动学描述为

式中, (u,v)为当前构型下物质点r的参数坐标,ei为参数 (u,v)所属的第i个T 样条单元的单元节点向量,Si(u,v)是第i个T 样条单元的形函数.

1.3 T 样条薄板的弹性模型

基于T 样条曲面的基尔霍夫薄板的总弹性能可表示为[30]

式中,Umid为只与中性层应变 εmid有关的薄板中性层的能量,Uκ为只与弯曲应变 κ 有关的弯曲应变能.E为各向同性线弹性材料的广义胡克定律弹性系数矩阵.根据文献[31], εmid和 κ 可以表示为

式中,r和r0分别表示当前构型和初始构型中,薄板中性层上任意点P(u,v)的坐标.表示中性层的单位法向量.T0表示中性层的局部笛卡尔坐标系(e0)1−(e0)2−(e0)3与曲 面坐标系(g0)1−(g0)2−n0的变换矩阵,其表达式为

弹性力Qs是弹性能U的相对于单元节点向量e的导数.对于一个基于T 样条曲面单元的薄板系统,其弹性力表达式为

切线刚度矩阵JQs是弹性力Qs对e的导数,其表达式为

图1 给出了等几何基尔霍夫薄板的转换示意图.首先,体积积分简化为厚度h与中性层面积分的乘积.然后,当前构型下物理空间中的面微元 dA被映射到参数空间中的面微元 dA0,其定义为

图1 等几何基尔霍夫薄板积分过程Fig.1 The process of isogeometric Kirchhoff thin plate integration

最后,对薄板的积分表达式可以写为

式中, ωi和 ωj是积分点u和v所对应的积分系数;J1是参数空间到父单元的转换矩阵的行列式值,J1的表达式为

2 网格局部细化

与NURBS 单元网格相比,T 样条允许在局部细化单元网格.在处理碰撞接触等局部发生应变的剧烈变化的问题时具有优势.本节将讨论T 样条薄板的网格局部细分及新控制点的计算方法.

如图2 所示,一个边界为 [ui,ui+1]×[vi,vi+1]的T 样条单元被均匀分为4 个单元,将会产生两个新的结点,umid=(ui+ui+1)/2 和vmid=(vi+vi+1)/2.细化后,初始T 样条单元的区域将被4 个新的子单元所替换,其单元边界分别为以及[umid,ui+1]每一个T 样条单元可以被循环细化直到达到停止细化的指标.

图2 局部细化策略的示意图Fig.2 Schematic diagram of local refinement strategy

曲面H初始由一个双三次T 样条S1表示,它的控制点数量为m,笛卡尔空间中的控制点矩阵P的表达式在式(15)中给出.由于T 网格中一个交点对应一个控制点,细化后的几何模型将生成新的控制点,所以细化后的曲面S2中控制点个数增至n,其控制点坐标为P˜.初始T 样条曲面S1被称为是细化后的T 样条曲面S2的子空间(表示为S1⊂S2)

根据参考文献[18],可以得到T 样条混合函数的局部细化算法:初始T 样条曲面上任意控制点Pi的局部结点矢量为则由ui和vi得到的B 样条基函数Ni(u) 和Ri(v)分别为所以其混合函数Bi可表示为Bi=Ni(u)·Ri(v).如图2所示,初始T 样条曲面被细化后,两个新的结点=umid和=vmid被插入到原始单元的边界上.控制点Pi所对应的基函数Ni(u) 和Ri(v)将使用插入u和v的新的局部结点矢量表示.因此,基函数Ni(u)可以写 作之和,Ri(v) 可以写作与之和.一个初始控制点Pi的混合函数Bi可以由表示

式中,初始T 样条曲面S1中的第i个混合函数Bi被分解为四个新混合函数与其对应系数乘积之和.如果等于细化后T 样条曲面第j个混合函数那么S1中的混合函数Bi可以被写作S2中混合函数的线性组合

因为S1⊂S2,由参数 (u,v)定义的S1上的点r1和S2上的点r2是相等的,也就是r1(u,v)≡r2(u,v).根据式(1)和式(2),r1和r2可表示为

根据Bi与B˜j之间的转换关系,可得到之间的转换关系,转换方程如式(19)和式(20)所示

式中,P4和P˜4是分别由控制点组成的矩阵.对于转换矩阵Mm,其第j行第i列为式(17)中的系数.

根据式(20),可以得到细化后T 样条的齐次坐标矩阵然而,本文中T 样条单元的节点向量是由控制点笛卡尔坐标构成的而非齐次坐标.因此,控制点的笛卡尔坐标和权重的计算表达示为

综上所述,根据本节给出的局部细化策略和算法,可以根据原T 样条曲面和插入的结点得到细化后的T 样条曲面,实现几何模型拓扑结构的变化.

3 变网格T 样条薄板的求解算法

本文中基于T 样条的柔性等几何薄板系统的动力学方程可表示为

图3 变网格系统动力学分析流程图Fig.3 Dynamic analysis flow chart of variable mesh system

然后根据转换矩阵Mm得到细化后新系统的节点坐标,节点速度和节点加速度,如下式(24)所示

更新几何模型后,约束方程也会发生变化,所以将会得到新的布尔矩阵.根据可以得到新系统的广义坐标广 义 速 度和广义加速度其计算公式如下量将会得到更新.通过与的乘积,将会得到新系统的广义质量阵,广义外力,广义弹性力和广

获取了上述变量之后,系统动力学方程中的变义接触力.将更新后的动力学方程输入到广义α 算法中进行求解,便可得到网格局部细化、自由度增加后的新系统下一时间步的广义变量实现了对变自由度系统动力学方程的求解.

4 数值算例

4.1 受端弯矩作用的悬臂薄板

图4 为一末端受弯矩M的悬臂薄板.该薄板长度L=12 m, 宽度b=1 m , 厚度h=0.1 m, 弹性模量E=1.2 MPa,泊松比 ν =0,惯性矩图5 给出了不同端弯矩下的薄板构型图.根据文献[32] 可知,取最大弯矩为当末端所受弯矩M=M0,此时薄板将弯曲成一个半径为的圆.当末端所受弯矩为M=0.75M0,0.5M0和 0 .25M0时,薄板的端截面转角分别为 1 .5π , π和 0 .5π.

图4 受端弯矩作用的悬臂薄板Fig.4 Cantilever subjected to end bending moment

图5 端弯矩作用下,变形后悬臂梁的构形图Fig.5 Configuration of deformed cantilever under external moments

4.2 受内压圆环

一平面圆环构件,在其中心孔边缘受沿半径向外方向的均布载荷P=10 MPa.根据对称性,取圆环构件的四分之一并在两侧施加滑动边界约束,如图6(a)所示.圆环构件的内与外径分别为R1=1 m和R2= 2 m,厚度h= 0.1 m.材料密度ρ=7800kg/m3,杨氏模量E=210 GPa,泊松比 ν =0.3.

图6 四分之一圆环构件及其局部细化示意图Fig.6 Diagram of a quarter of circular thin plate and its local refinement diagram

采用T 样条曲面单元对构件离散,单元网格分布图如图6(b) 所示.因为最大米塞斯应力出现在AB边界附近,内部区域被局域细化.模型共含有28 个单元和72 个控制点.图7 给出了局部细化后构件的米塞斯应力云图.

图7 局部细化后构件米塞斯应力云图Fig.7 von-Mises stress nephogram after local refinement

因为该算例属于平面问题,所以轴向应力 σz=0.根据拉美方程可以得到圆环在任意半径r处米塞斯应力 σs的理论解

式中, σr为径向应力, σθ为环向应力,K为圆环外径与内径之比,K=R2:R1.

在有限元软件ANSYS 中采用SHELL63 单元运行相同算例并对比构件中最大和最小应力结果如表1 所示.可以看出,相较于ANSYS 结果,IGA 结果更加接近于理论解而且与理论解的误差小于1‰.

表1 圆环构件最小和最大米塞斯应力Table 1 Minimum and maximum von-Mises stress on circular thin plate

4.3 柔性单摆

该算例测试了基于T 样条的等几何薄板的动力学特性.如图8 所示,一个柔性单摆铰接固定在A 处.柔性薄板单摆在重力加速度g=9.81 m/s2的作用下坠落.薄板的参数在表2 中给出.

表2 材料参数Table 2 Material parameters

图8 一端铰接的柔性薄板单摆Fig.8 Flexible thin plate pendulum

图9 给出了含有2 × 2 个T 样条单元的薄板在1 s 内动能Ek、重力势能Eg、弹性势能Ee及总能量Et的变化情况.可以看出,基于T 样条的等几何薄板满足能量守恒定律.为了检验本文提出的方法的收敛性,分别使用2 × 2,4 × 4 和8 × 8 的T 样条曲面单元对图8 所示的单摆进行离散.图10 给出了单摆自由端C 的z向位移曲线.可以看出,基于T 样条的等几何薄板具有较好的收敛性.图11 给出了含有2 ×2 个T 样条单元的单摆在1 s 内的连续构型.

图9 能量曲线Fig.9 The energy balance curve

图10 等几何薄板中T 样条单元的收敛性Fig.10 Convergence of T-spline element in thin plate

图11 含有2 × 2 个T 样条单元的单摆的连续构型Fig.11 Configuration of the pendulum with 2 × 2 T-spline elements

4.4 刚性球与柔性板的接触

本算例将T 样条局部细化算法应用于柔体动力学分析中.如图12 所示,一个刚性球从四边固定的薄板中心正上方自由下落.接触发生后,对薄板中心受冲击区域进行局部细化.刚性球与板的参数表3中给出.

表3 接触算例的参数Table 3 Parameters of the contact example

图12 自由坠落刚性球与柔性薄板的碰撞Fig.12 The collision between a rigid ball and a flexible thin plate

在本算例中,采用罚值法来完成刚体球与柔性板的接触实现.在每个T 样条单元中,刚体球与均匀检测点的距离用式(27)表示

式中,p为嵌入深度,r(u,v)为薄板中性层上的接触检测点,rs为任意点检测点在薄板上表面所对应的物质点,rc为球心位置,h为薄板厚度,R为球心半径,n表示接触检测点的单位法向量.

一旦小球与薄板发生接触,开始计算薄板所受的接触力.计算切向接触力时,采用文献[33]中计及临界滑动速度v0的平滑化库伦摩擦模型.接触探测点i的法向接触力Fin和切向接触力Fit表达式分别为

式中,vp为嵌入速度,k和c表示刚度系数与阻尼系数.在式中, µ 为薄板与球之间的摩擦系数,vt为切向接触速度,v0为假定的临界滑动速度,vet为单元切向接触速度.

图13 介绍了三种网格构型图.网格1 和网格2 都对表面进行了均匀的网格细化,其单元数量分别为14 × 14 和10 × 10.网格3 的初始网格较为粗糙,含有10 × 10 个单元.当接触发生后,薄板的接触区域被局部细化,单元数量由100 变为112 个.三组薄板的详细信息见表4.

图13 三组不同的网格构型Fig.13 The mesh refinement for three groups

表4 三组薄板的详细信息Table 4 The detailed information of three groups

图14 给出了三组小球球心在1 s 内的z向位移曲线.图15 和图16 则为三组薄板的质心z向位移曲线与弹性能曲线及其局部放大图.在t=0.216 s时,刚性球与柔性板首次接触,对薄板进行如图13所示的3 种不同网格细化策略.三组球的质心垂直位移具有较好的一致性,薄板质心的垂直位移曲线在1 s 内也呈现相似的趋势.从薄板的质心位移与弹性能计算结果来看,允许局部网格细化的T 样条等几何薄板可以达到与退化为NURBS 曲面的均匀网格板相同的计算精度,且相较于网格2,局部细化的网格3 所对应曲线的变化趋势更接近网格1.

图14 球心z 向位移Fig.14 Vertical displacement of the center of the ball

图15 薄板质心的z 向位移Fig.15 Vertical displacement of the plate’s centroid

图16 薄板弹性能曲线Fig.16 The elastic energy curve of thin plate

为了比较不同网格构型的薄板对计算效率的影响,图17 给出了三组薄板的仿真时间柱状图.

图17 计算消耗时间Fig.17 Time consumption for three groups

三组算例均是在配备Intel Core i5 CPU 和8GB RAM 的笔记本电脑上执行的.由图17 可知,网格2 与网格3 用时接近,约为网格1 计算用时的一半.值得注意的是,由于在接触区域附近执行了网格细化,使得求解收敛更快,网格3 虽然自由度数略多于网格2,用时反而更少.这也进一步体现了T 样条单元网格局部细化的优势.

图18 展示了带有112 个T 样条单元局部细化薄板在不同时刻的米塞斯应力云图,分别为t=0.18 s接触前的薄板发生最大变形,t=0.218 s发生接触后首次使用细化后模型的时刻,t=0.548 s接触后的发生最大变形的时刻以及仿真结束时刻t=1 s.

图18 局部细化薄板的米塞斯应力云图Fig.18 von-Mises stress distribution of the locally refined thin plate

5 结论

本文提出了一种基于T 样条曲面的等几何分析方法,建立了使用T 样条曲面单元离散的基尔霍夫薄板模型.单元基函数和节点坐标分别为T 样条基函数和控制点坐标,无需网格划分,既保证模型的几何精确,又能在没有约束方程的情况下保证期望的连续性条件.给出了基于T 样条的薄板运动学模型、弹性力模型及其雅克比矩阵的计算方法.为了体现T 样条局部细化特性在减少单元数目、提高计算效率方面的优势,提出了一种基于T 样条单元的局部网格细化算法,包括单元网格拓扑的改变和相应的新控制点坐标计算方法.将该细化算法与广义α 法相结合,建立了带有网格局部细化的变自由度系统动力学方程的求解算法.通过受端弯矩的悬臂薄板以及环形构件受内压的静力学算例证明了T 样条单元弹性力模型的正确性.柔性摆算例验证了T 样条单元在动力学问题中的收敛性和机械能守恒特性.最后,为验证T 样条局部细化特性在模拟接触碰撞等柔性体局部发生应变剧烈变化等问题中的优势,建立了刚性球落在柔性板上的动力学实例并进行了仿真.首先,采用10 × 10 单元网格对柔性板进行离散.接触发生后,对冲击区域进行局部细化,得到112 个单元的网格.对10 × 10 和14 × 14 两种网格进行了仿真,并将仿真结果作为基准.可以看到,112 单元网格的计算结果与14 × 14 网格具有良好的一致性,但所消耗的时间只有其1/2.以上算例证明了基于T 样条的局部网格细化等几何薄板在柔性多体系统动力学分析中的应用价值.