DYNAMIC ANALYSIS OF VARIABLE MESH ISOGEOMETRIC THIN PLATE BASED ON T-SPLINE
-
摘要: 具有大位移、大变形的薄板在接触碰撞等工况下, 其局部应变会产生剧烈变化. 为了保证对其进行动力学分析的精度和计算效率, 本文整合计算机辅助设计(CAD)与计算机辅助工程(CAE)系统, 提出了一种基于T样条曲面的变网格柔性系统等几何分析方法. 首先, 建立基于T样条曲面单元的基尔霍夫薄板运动学模型, 并根据非线性格林−拉格朗日应变建立由T样条曲面单元离散的薄板弹性模型. 其次, 通过在T网格中的局部区域插入节点的方式, 达到T样条曲面网格局部更新的目的. 利用T样条混合函数细化算法得到计算新广义变量的转换矩阵, 并结合广义α法创建了变自由度系统动力学方程的求解算法, 形成了系统的T样条单元局部细化算法. 最后, 静力学算例与柔性单摆模型分别验证了T样条薄板弹性模型的正确性, 以及T样条薄板单元在动力学分析上的精度和收敛性. 通过对受冲击柔性薄板的动力学分析表明, 本文所提出T样条单元及局部细化算法可以只在接触碰撞等应变剧烈变化的区域实现局部网格细化, 从而控制系统自由度数, 提高计算效率.Abstract: The local strain of a thin plate with large displacement and large deformation will change dramatically under contact and collision working conditions. In order to ensure the accuracy and computational efficiency of flexible thin plate system’s dynamic analysis, this investigation integrates computer aided design (CAD) and computer aided engineering (CAE) systems, and proposes an isogeometric analysis (IGA) method for variable mesh flexible multibody system based on T-spline surface elements. Firstly, the kinematic description of a Kirchhoff thin plate based on T-spline surface elements is modeled, and the elastic model of a thin plate discretized by T-spline surface elements is established according to the nonlinear Green−Lagrange strain. Secondly, the goal of updating mesh of T-spline surface locally is achieved by inserting knots into the local region of the corresponding T-mesh. The transformation matrix, which is used to calculate the new generalized coordinates, generalized velocities and generalized accelerations of the refined system, is obtained by using the T-spline blending function refinement algorithm. The calculating solution algorithm for the dynamic equation of the system with variable degrees of freedom is created by combining the generalized α method with geometry update routine, and thus the local mesh refinement algorithm for the surface which is modeled by T-spline is formed. Finally, statics examples and flexible pendulum model verify the correctness for the elastic model of Kirchhoff thin plate based on T-spline surface, as well as computation precision and convergence for the proposed method in the dynamics analysis respectively. The dynamic analysis of the impacted flexible thin plate shows that the T-spline element and local refinement algorithm proposed in this paper can realize local mesh update only in the area where the strain changes violently, such as contact and collision, so as to control the degree of freedom of the system and improve the computational efficiency.
-
引言
随着柔性轻质结构在航天、汽车以及船舶工程等领域日益广泛的应用, 柔性多体系统的动力学分析也愈加重要. 然而, 以线弹性小变形假设为基础的, 将节点位移和转动作为坐标的传统有限单元不适合解决存在大位移、大变形的柔性多体系统的动力学问题. 此外, 对于传统结构有限元, 将几何数据转换为有限元网格数据是困难且耗时的[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]{{\boldsymbol{S}}_{\rm{T}}}(u,v) = \sum\limits_{i = 1}^n {{{\boldsymbol{P}}_i}{T_i}(u,v)} (1) 式中,
{{\boldsymbol{S}}_{\rm{T}}} 为T样条曲面上物质坐标为\left( {u,v} \right) 的点的全局坐标,{{\boldsymbol{P}}_i} = \left( {{x_i},{y_i},{z_i}} \right) 为第i 个控制点的全局坐标,n 为控制点的数量,{T_i} 为第i 个控制点所对应的有理T样条基函数, 其表达式为{T_i}(u,v) = \frac{{{\omega _i}{B_i}(u,v)}}{{\displaystyle\sum\limits_{i = 1}^n {{\omega _i}{B_i}(u,v)} }}{\rm{ = }}\frac{{{\omega _i}{N_i}(u){R_i}(v)}}{{W(u,v)}} (2) 式中, 混合函数
{B_i} 为B样条基函数{N_i}\left( u \right) 和{R_i}\left( v \right) 的乘积.{\omega _i} 是第i 个控制点所对应的权重,W 为权值与基函数的乘积之和.与NURBS不同的是, 每一个T样条基函数都是基于局部结点矢量而不是全局结点矢量来计算的. 局部结点矢量的建立依赖于T网格和锚点的建立, 根据这两者可以得到T样条曲面上每一个控制点所对应的局部结点矢量
{{\boldsymbol{u}}_i} = \left[ {{u_1},{u_2}, \cdots ,{u_{p + 2}}} \right] 和{{\boldsymbol{v}}_i} = \left[ {{v_1},{v_2}, \cdots ,{v_{q + 2}}} \right] . T网格和锚点的具体定义参见文献[26-27], 并可依据参考文献[28]得到三次T样条基函数以及其导数的表达式.1.2 T样条薄板的运动学描述
由于基尔霍夫薄板忽略横向剪切变形, 法向量始终垂直于中性层, 其运动学描述可以简化为其中性层的运动学描述. 为此, 首先需要采用T样条单元对薄板中性层进行运动学描述.
T样条单元是由T样条基函数的简化连续线构成的物理区域[29]. 在等几何分析中, 利用IEN数组可以确定每个样条单元的连接关系. IEN中的I表示数组是整数值, EN是“element nodes”的首字母缩写, 表示“单元节点”, 该数组可以将局部基函数号和单元号映射到相应的全局控制点号[29]. 对于第
i 个T样条单元而言, 其形函数以及单元节点向量{{\boldsymbol{e}}_i} 可以表示为{{\boldsymbol{e}}_i} = {\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{P}}_{{\rm{IEN}}\left( 1 \right)}}}&{{{\boldsymbol{P}}_{{\rm{IEN}}\left( 2 \right)}}}& \cdots &{{{\boldsymbol{P}}_{{\rm{IEN}}\left( j \right)}}}& \cdots &{{{\boldsymbol{P}}_{{\rm{IEN}}\left( n \right)}}} \end{array}} \right]^{\rm{T}}} (3) \begin{split} {{\boldsymbol{S}}_i}\left( {u,v} \right) = &\left[ {{T_{{\rm{IEN}}\left( 1 \right)}}\left( {u,v} \right){\boldsymbol{I}}}\qquad{{T_{{\rm{IEN}}\left( 2 \right)}}\left( {u,v} \right){\boldsymbol{I}}}\qquad \cdots \right. \\ &\left. {{T_{{\rm{IEN}}\left( j \right)}}\left( {u,v} \right){\boldsymbol{I}}}\qquad \cdots \qquad{{T_{{\rm{IEN}}\left( n \right)}}\left( {u,v} \right){\boldsymbol{I}}} \right] \end{split}\;\;\;\; (4) 式中,
n 是一个T样条单元中包含的控制点的个数,{\rm{IEN}}\left( j \right) 为该T样条单元的IEN数列中的第j 项,{{\boldsymbol{P}}_{{\rm{IEN}}\left( j \right)}} 为第{\rm{IEN}}\left( j \right) 个控制点的坐标,{\boldsymbol{I}} 为3 \times 3 的单位矩阵,{T_{{\rm{IEN}}\left( j \right)}}\left( {u,v} \right) 为第{\rm{IEN}}\left( j \right) 个控制点所对应的有理T样条基函数.综上所述, 基尔霍夫薄板的运动学描述为
{\boldsymbol{r}}\left( {u,v} \right) = {{\boldsymbol{S}}_i}\left( {u,v} \right) \cdot {{\boldsymbol{e}}_i} (5) 式中,
\left( {u,v} \right) 为当前构型下物质点{\boldsymbol{r}} 的参数坐标,{{\boldsymbol{e}}_i} 为参数\left( {u,v} \right) 所属的第i 个T样条单元的单元节点向量,{{\boldsymbol{S}}_i}\left( {u,v} \right) 是第i 个T样条单元的形函数.1.3 T样条薄板的弹性模型
基于T样条曲面的基尔霍夫薄板的总弹性能可表示为[30]
U = {U_{{\rm{mid}}}} + {U_{\rm{\kappa }}} = \frac{1}{2}\int\limits_V {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}^{\rm{T}}{\boldsymbol{E}}{{\boldsymbol{\varepsilon}} _{{\rm{mid}}}}{\rm{d}}V} + \frac{1}{2}\int\limits_V {{{\boldsymbol{\kappa}} ^{\rm{T}}}{\boldsymbol{E\kappa}} {\rm{d}}V} (6) 式中,
{U_{{\rm{mid}}}} 为只与中性层应变{{\boldsymbol{\varepsilon}} _{{\rm{mid}}}} 有关的薄板中性层的能量,{U_{\rm{\kappa }}} 为只与弯曲应变{\boldsymbol{\kappa}} 有关的弯曲应变能.{\boldsymbol{E}} 为各向同性线弹性材料的广义胡克定律弹性系数矩阵. 根据文献[31],{{\boldsymbol{\varepsilon}} _{{\rm{mid}}}} 和{\boldsymbol{\kappa}} 可以表示为\begin{split} & {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}}{\rm{ = }}\\ &\dfrac{1}{2}{\boldsymbol{T}}_0^{\rm{T}}\left( {\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {\boldsymbol{r}}}}{{\partial u}} \cdot \dfrac{{\partial {\boldsymbol{r}}}}{{\partial u}}}\!\!&\!\!{\dfrac{{\partial {\boldsymbol{r}}}}{{\partial u}} \cdot \dfrac{{\partial {\boldsymbol{r}}}}{{\partial v}}} \\ {\dfrac{{\partial {\boldsymbol{r}}}}{{\partial v}} \cdot \dfrac{{\partial {\boldsymbol{r}}}}{{\partial u}}}\!\!&\!\!{\dfrac{{\partial {\boldsymbol{r}}}}{{\partial v}} \cdot \dfrac{{\partial {\boldsymbol{r}}}}{{\partial v}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial u}} \cdot \dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial u}}}\!\!&\!\!{\dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial u}} \cdot \dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial v}}} \\ {\dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial v}} \cdot \dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial u}}}\!\!&\!\!{\dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial v}} \cdot \dfrac{{\partial {{\boldsymbol{r}}_0}}}{{\partial v}}} \end{array}} \right]} \right){{\boldsymbol{T}}_0}\\ & {\boldsymbol{\kappa}} =\\ & {\rm{ - }}z{\boldsymbol{T}}_0^{\rm{T}}\left( {\left[ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^2}{\boldsymbol{r}}}}{{\partial {u^2}}} \cdot {\boldsymbol{n}}}\!\!&\!\!{\dfrac{{{\partial ^2}{\boldsymbol{r}}}}{{\partial u\partial v}} \cdot {\boldsymbol{n}}} \\ {\dfrac{{{\partial ^2}{\boldsymbol{r}}}}{{\partial u\partial v}} \cdot {\boldsymbol{n}}}\!\!&\!\!{\dfrac{{{\partial ^2}{\boldsymbol{r}}}}{{\partial {v^2}}} \cdot {\boldsymbol{n}}} \end{array}} \right] \!-\! \left[ {\begin{array}{*{20}{c}} {\dfrac{{{\partial ^2}{{\boldsymbol{r}}_0}}}{{\partial {u^2}}} \cdot {{\boldsymbol{n}}_0}}\!\!&\!\!{\dfrac{{{\partial ^2}{{\boldsymbol{r}}_0}}}{{\partial u\partial v}} \cdot {{\boldsymbol{n}}_0}} \\ {\dfrac{{{\partial ^2}{{\boldsymbol{r}}_0}}}{{\partial u\partial v}} \cdot {{\boldsymbol{n}}_0}}\!\!&\!\!{\dfrac{{{\partial ^2}{{\boldsymbol{r}}_0}}}{{\partial u\partial v}} \cdot {{\boldsymbol{n}}_0}} \end{array}} \right]} \right){{\boldsymbol{T}}_0} \end{split} 式中,
{\boldsymbol{r}} 和{{\boldsymbol{r}}_0} 分别表示当前构型和初始构型中, 薄板中性层上任意点{\boldsymbol{P}}\left( {u,v} \right) 的坐标.{\boldsymbol{n}} = \dfrac{{{{\boldsymbol{r}}_u} \times {{\boldsymbol{r}}_v}}}{{\left| {{{\boldsymbol{r}}_u} \times {{\boldsymbol{r}}_v}} \right|}} 表示中性层的单位法向量.{{\boldsymbol{T}}_0} 表示中性层的局部笛卡尔坐标系{\left( {{{\boldsymbol{e}}_0}} \right)_1}{\rm{ - }}{\left( {{{\boldsymbol{e}}_0}} \right)_2}{\rm{ - }}{\left( {{{\boldsymbol{e}}_0}} \right)_3} 与曲面坐标系{\left( {{{\boldsymbol{g}}_0}} \right)_1}{\rm{ - }}{\left( {{{\boldsymbol{g}}_0}} \right)_2}{\rm{ - }}{{\boldsymbol{n}}_0} 的变换矩阵, 其表达式为\left. {\begin{array}{*{20}{l}}{\left({\boldsymbol{g}}_{0}\right)}_{1}=\dfrac{\partial {\boldsymbol{r}}_{0}}{\partial u}\begin{array}{c},\end{array}{\left({\boldsymbol{g}}_{0}\right)}_{{2}}=\dfrac{\partial {\boldsymbol{r}}_{0}}{\partial v}\\ {\boldsymbol{n}}_{0}=\dfrac{{\left({\boldsymbol{g}}_{0}\right)}_{1}\times {\left({\boldsymbol{g}}_{0}\right)}_{2}}{\left|{\left({\boldsymbol{g}}_{0}\right)}_{1}\times {\left({\boldsymbol{g}}_{0}\right)}_{2}\right|}\end{array}} \right\} (7) \left. {\begin{array}{*{20}{l}} {{{\left( {{{\boldsymbol{e}}_0}} \right)}_1} = {{\left( {{{\boldsymbol{g}}_0}} \right)}_1}} \\ {{{\left( {{{\boldsymbol{e}}_0}} \right)}_{\rm{3}}} = {{\boldsymbol{n}}_0}} \\ {{{\left( {{{\boldsymbol{e}}_0}} \right)}_{\rm{2}}} = {{\left( {{{\boldsymbol{e}}_0}} \right)}_{\rm{3}}} \times {{\left( {{{\boldsymbol{e}}_0}} \right)}_{\rm{1}}}} \end{array}} \right\} (8) \;\;\;{{\boldsymbol{T}}_0} = {\left[ {\begin{aligned} {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \cdot {{\left( {{{\boldsymbol{e}}_0}} \right)}_1}} \quad {{{\left( {{{\boldsymbol{g}}_0}} \right)}_1} \cdot {{\left( {{{\boldsymbol{e}}_0}} \right)}_2}} \\ {{{\left( {{{\boldsymbol{g}}_0}} \right)}_2} \cdot {{\left( {{{\boldsymbol{e}}_0}} \right)}_1}} \quad {{{\left( {{{\boldsymbol{g}}_0}} \right)}_2} \cdot {{\left( {{{\boldsymbol{e}}_0}} \right)}_2}} \end{aligned}}\right]^{ - {\rm{T}}}} (9) 弹性力
{{\boldsymbol{Q}}_{\rm{s}}} 是弹性能U 的相对于单元节点向量{\boldsymbol{e}} 的导数. 对于一个基于T样条曲面单元的薄板系统, 其弹性力表达式为{{\boldsymbol{Q}}_{\rm{s}}} = \frac{{\partial U}}{{\partial {\boldsymbol{e}}}}{\rm{ = }}\int_V {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}^{\rm{T}}{\boldsymbol{E}}\frac{{\partial {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}}}}{{\partial {\boldsymbol{e}}}}} {\rm{d}}V{\rm{ + }}\int_V {{{\boldsymbol{\kappa}} ^{\rm{T}}}{\boldsymbol{E}}\frac{{\partial {\boldsymbol{\kappa}} }}{{\partial {\boldsymbol{e}}}}} {\rm{d}}V (10) 切线刚度矩阵
{\boldsymbol{J}}{{\boldsymbol{Q}}_{\rm{s}}} 是弹性力{{\boldsymbol{Q}}_{\rm{s}}} 对{\boldsymbol{e}} 的导数, 其表达式为\begin{split} {\boldsymbol{J}}{{\boldsymbol{Q}}_{\rm{s}}} =& \int_V {\left( {{{\dfrac{{\partial {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}^{\rm{T}}}}}{{\partial {{\boldsymbol{e}}_j}}}}}{\boldsymbol{E}}\dfrac{{\partial {{\boldsymbol{\varepsilon}} _{{\rm{mid}}}}}}{{\partial {{\boldsymbol{e}}_i}}} + {\boldsymbol{\varepsilon}} _{{\rm{mid}}}^{\rm{T}}{\boldsymbol{E}}\dfrac{{{\partial ^2}{{\boldsymbol{\varepsilon}} _{{\rm{mid}}}}}}{{\partial {{\boldsymbol{e}}_j}\partial {{\boldsymbol{e}}_i}}}} \right){\rm{d}}V} +\\ & \int_V {\left( {{{\dfrac{{\partial {\boldsymbol{\kappa}}^{\rm{T}} }}{{\partial {{\boldsymbol{e}}_j}}}}}{\boldsymbol{E}}\dfrac{{\partial {\boldsymbol{\kappa}} }}{{\partial {{\boldsymbol{e}}_i}}} + {{\boldsymbol{\kappa}} ^{\rm{T}}}{\boldsymbol{E}}\dfrac{{{\partial ^2}{\boldsymbol{\kappa}} }}{{\partial {{\boldsymbol{e}}_j}\partial {{\boldsymbol{e}}_i}}}} \right){\rm{d}}V} \end{split} (11) 图1给出了等几何基尔霍夫薄板的转换示意图. 首先, 体积积分简化为厚度
h 与中性层面积分的乘积. 然后, 当前构型下物理空间中的面微元{\rm{d}}A 被映射到参数空间中的面微元{\rm{d}}{A_0} , 其定义为{\rm{d}}A{\rm{ = }}\left\| {{{\boldsymbol{r}}_u} \times {{\boldsymbol{r}}_v}} \right\|{\rm{d}}{A_0} = {J_2}{\rm{d}}{A_0} (12) 最后, 对薄板的积分表达式可以写为
\int_V f {\rm{d}}V = h{J_2}\int_0^v {\int_0^u {f{\rm{d}}u{\rm{d}}v = } } {J_1}{J_2}h\sum\limits_{j = 1}^m {\sum\limits_{i = 1}^n {{\omega _i}{\omega _j}f\left( {\hat u,\hat v} \right)} } (13) 式中,
{\omega _i} 和{\omega _j} 是积分点\hat u 和\hat v 所对应的积分系数;{J_1} 是参数空间到父单元的转换矩阵的行列式值,{J_1} 的表达式为{J_1} = \frac{{{u_i} - {u_{i - 1}}}}{2} \cdot \frac{{{v_i} - {v_{i - 1}}}}{2} (14) ![]()
图 1 等几何基尔霍夫薄板积分过程Figure 1. The process of isogeometric Kirchhoff thin plate integration2. 网格局部细化
与NURBS单元网格相比, T样条允许在局部细化单元网格. 在处理碰撞接触等局部发生应变的剧烈变化的问题时具有优势. 本节将讨论T样条薄板的网格局部细分及新控制点的计算方法.
如图2所示, 一个边界为
\left[ {{u_i},{u_{i + 1}}} \right] \times \left[ {{v_i},{v_{i + 1}}} \right] 的T样条单元被均匀分为4个单元, 将会产生两个新的结点,{u_{{\rm{mid}}}} = {{\left( {{u_i} + {u_{i + 1}}} \right)} / 2} 和{v_{{\rm{mid}}}} = {{\left( {{v_i} + {v_{i + 1}}} \right)} / 2} . 细化后, 初始T样条单元的区域将被4个新的子单元所替换, 其单元边界分别为\left[ {u_i}, {u_{{\rm{mid}}}} \right] \times \left[ {{v_j},{v_{{\rm{mid}}}}} \right] ,\left[ {u_{{\rm{mid}}}},\right. \left.{ u _{i + 1}} \right] \times \left[ {{v_j},{v_{{\rm{mid}}}}} \right] ,\left[ {{u_i},{u_{{\rm{mid}}}}} \right] \times \left[ {{v_{{\rm{mid}}}},{v_{j + 1}}} \right] 以及\left[ {{u_{{\rm{mid}}}},{u_{i{\rm{ + }}1}}} \right] \times \left[ {{v_{{\rm{mid}}}},{v_{j + 1}}} \right] . 每一个T样条单元可以被循环细化直到达到停止细化的指标.曲面
H 初始由一个双三次T样条{S_1} 表示, 它的控制点数量为m , 笛卡尔空间中的控制点矩阵{\boldsymbol{P}} 的表达式在式(15)中给出. 由于T网格中一个交点对应一个控制点, 细化后的几何模型将生成新的控制点, 所以细化后的曲面{S_2} 中控制点个数增至n , 其控制点坐标为\tilde {\boldsymbol{P}} . 初始T样条曲面{S_1} 被称为是细化后的T样条曲面{S_2} 的子空间(表示为{S_1} \subset {S_2} ){\boldsymbol{P}} = {\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{P}}_1^{\rm{T}}}&{{\boldsymbol{P}}_2^{\rm{T}}}& \cdots &{{\boldsymbol{P}}_i^{\rm{T}}}& \cdots &{{\boldsymbol{P}}_m^{\rm{T}}} \end{array}} \right]^{\rm{T}}} (15) 根据参考文献[18], 可以得到T样条混合函数的局部细化算法: 初始T样条曲面上任意控制点
{{\boldsymbol{P}}_i} 的局部结点矢量为{{\boldsymbol{u}}_i} = \left[ {u_i^0,u_i^1,u_i^2,u_i^3,u_i^4} \right] 和{{\boldsymbol{v}}_i} = \left[ v_i^0,v_i^1,v_i^2,\right. \left.v_i^3,v_i^4 \right] , 则由{{\boldsymbol{u}}_i} 和{{\boldsymbol{v}}_i} 得到的B样条基函数{N_i}\left( u \right) 和{R_i}\left( v \right) 分别为{N_i}\left[ {u_i^0,u_i^1,u_i^2,u_i^3,u_i^4} \right]\left( u \right) 和{R_i}\left[ {v_i^0,v_i^1,v_i^2,v_i^3,v_i^4} \right]\left(v \right) , 所以其混合函数{B_i} 可表示为{B_i} = {N_i}\left( u \right) \cdot {R_i}\left( v \right) . 如图2所示, 初始T样条曲面被细化后, 两个新的结点\hat u = {u_{{\rm{mid}}}} 和\hat v = {v_{{\rm{mid}}}} 被插入到原始单元的边界上. 控制点{{\boldsymbol{P}}_i} 所对应的基函数{N_i}\left( u \right) 和{R_i}\left( v \right) 将使用插入\hat u 和\hat v 的新的局部结点矢量表示. 因此, 基函数{N_i}\left( u \right) 可以写作{{c_1}}\tilde N_i^1\left( u \right) 与{c_2}\tilde N_i^2\left( u \right) 之和,{R_i}\left( v \right) 可以写作{d_{ 1}}\tilde R_i^1\left( v \right) 与{d_2}\tilde R_i^2\left( v \right) 之和. 一个初始控制点{{\boldsymbol{P}}_i} 的混合函数{B_i} 可以由\tilde N_i^s\left( u \right) 和\tilde R_i^r\left( v \right) 表示{B_i}\left( {u,v} \right){\rm{ = }}\sum\limits_{s = 1}^2 {\sum\limits_{r = 1}^2 {\left( {{c_s} \cdot {d_r}} \right)\tilde N_i^s\tilde R_i^r} } {\rm{ = }}\sum\limits_{t = 1}^4 {f_i^t\tilde B_i^t\left( {u,v} \right)} (16) 式中, 初始T样条曲面
{S_1} 中的第i 个混合函数{B_i} 被分解为四个新混合函数\tilde B_i^t 与其对应系数f_i^t 乘积之和. 如果\tilde B_i^t 等于细化后T样条曲面第j 个混合函数{\tilde B_j} , 即\tilde B_i^t{\rm{ = }}{\tilde B_j} , 那么{S_1} 中的混合函数{B_i} 可以被写作{S_2} 中混合函数{\tilde B_j} 的线性组合{B_i}\left( {u,v} \right) = \sum\limits_{j = 1}^n {f_i^j{{\tilde B}_j}\left( {u,v} \right)} (17) 因为
{S_1} \subset {S_2} , 由参数\left( {u,v} \right) 定义的{S_1} 上的点{{\boldsymbol{r}}_1} 和{S_2} 上的点{{\boldsymbol{r}}_2} 是相等的, 也就是{{\boldsymbol{r}}_1}\left( {u,v} \right) \equiv {{\boldsymbol{r}}_2}\left( {u,v} \right) . 根据式(1)和式(2),{{\boldsymbol{r}}_1} 和{{\boldsymbol{r}}_2} 可表示为\left. {\begin{array}{*{20}{l}} {{{\boldsymbol{r}}_1}\left( {u,v} \right){\rm{ = }}\displaystyle\sum\limits_{i = 1}^m {\frac{{{{\boldsymbol{P}}_i}{\omega _i}{B_i}\left( {u,v} \right)}}{{W\left( {u,v} \right)}}} } \\ {{{\boldsymbol{r}}_2}\left( {u,v} \right){\rm{ = }}\displaystyle\sum\limits_{j = 1}^n {\frac{{{{\tilde {\boldsymbol{P}}}_j}{{\tilde \omega }_j}{{\tilde B}_j}\left( {u,v} \right)}}{{\tilde W\left( {u,v} \right)}}} } \end{array}} \right\} (18) 根据
{B_i} 与{\tilde B_j} 之间的转换关系, 可得到\tilde {\boldsymbol{P}}_j^4 和{\boldsymbol{P}}_i^4 之间的转换关系, 转换方程如式(19)和式(20)所示\sum\limits_{i = 1}^m {c_i^j{\boldsymbol{P}}_i^4 = \tilde {\boldsymbol{P}}_j^4} (19) 式中, 带有上标4的
{\boldsymbol{P}}_i^4 和\tilde {\boldsymbol{P}}_j^4 表示包括四维空间中的齐次坐标\left.\begin{split} & {\boldsymbol{P}}_i^4 = {{\boldsymbol{P}}_i}{\omega _i} = \left( {{\omega _i}x,{\omega _i}y,{\omega _i}z,{\omega _i}} \right)\\ & {{\tilde {\boldsymbol{P}}}^4} = {{\boldsymbol{M}}_{\rm{m}}}{{\boldsymbol{P}}^4} \\ & {{\boldsymbol{P}}^4}{\rm{ = }}{\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{P}}_1^4}&{{\boldsymbol{P}}_2^4}& \cdots &{{\boldsymbol{P}}_i^4}& \cdots &{{\boldsymbol{P}}_m^4} \end{array}} \right]^{\rm{T}}} \\ & {{\tilde {\boldsymbol{P}}}^4}{\rm{ = }}{\left[ {\begin{array}{*{20}{c}} {\tilde {\boldsymbol{P}}_1^4}&{\tilde {\boldsymbol{P}}_2^4}& \cdots &{\tilde {\boldsymbol{P}}_j^4}& \cdots &{\tilde {\boldsymbol{P}}_n^4} \end{array}} \right]^{\rm{T}}} \end{split}\right\} (20) 式中,
{{\boldsymbol{P}}^4} 和{\tilde {\boldsymbol{P}}^4} 是分别由控制点{\boldsymbol{P}}_i^4 和\tilde {\boldsymbol{P}}_j^4 组成的矩阵. 对于转换矩阵{{\boldsymbol{M}}_{\rm{m}}} , 其第j 行第i 列为式(17)中的系数f_i^j .根据式(20), 可以得到细化后T样条的齐次坐标矩阵
{\tilde {\boldsymbol{P}}^4} . 然而, 本文中T样条单元的节点向量是由控制点笛卡尔坐标构成的而非齐次坐标. 因此, 控制点{\tilde {\boldsymbol{P}}_j} 的笛卡尔坐标和权重的计算表达示为\left. {\begin{array}{*{20}{l}} {\tilde {\boldsymbol{\omega}} = {{\boldsymbol{M}}_{\rm{m}}}{\boldsymbol{\omega}} } \\ {{{\tilde {\boldsymbol{P}}}_j} = {{\tilde {\boldsymbol{P}}_j^4} / {{\omega _j}}}} \end{array}} \right\} (21) 式中
{\boldsymbol{\omega}} 和\tilde {\boldsymbol{\omega}} 分别为初始T样条和细化后T样条中全部控制点的权重向量,{\tilde \omega _j} 为向量\tilde {\boldsymbol{\omega}} 的第j 项,{\tilde {\boldsymbol{P}}_j} 是\tilde {\boldsymbol{P}} 的第j 项. 至此, 计算细化后T样条的新控制点的目标得以完成.综上所述, 根据本节给出的局部细化策略和算法, 可以根据原T样条曲面和插入的结点得到细化后的T样条曲面, 实现几何模型拓扑结构的变化.
3. 变网格T样条薄板的求解算法
本文中基于T样条的柔性等几何薄板系统的动力学方程可表示为
\tilde {\boldsymbol{M}}\ddot {\boldsymbol{q}} + {\tilde {\boldsymbol{Q}}_{\rm{s}}} - {\tilde {\boldsymbol{Q}}_{\rm{e}}} = {\bf{0}} (22) 式中,
\tilde {\boldsymbol{M}} 表示缩减自由度之后的系统质量阵,\ddot {\boldsymbol{q}} 表示广义加速度,{\tilde {\boldsymbol{Q}}_{\rm{s}}} 和{\tilde {\boldsymbol{Q}}_{\rm{e}}} 表示自由度缩减后的系统广义弹性力与包含接触力的广义外力. 需要指出的是, 引入了基于T样条的网格细化算法后, 式(22)所含有的自由度数量是时变的, 因此, 需要对现有求解方法进行改进, 以考虑拓扑结构的变化. 图3给出了变网格T样条薄板系统动力学分析的求解算法流程图, 图中{{\boldsymbol{q}}_0} ,{\dot {\boldsymbol{q}}_0} 和{\ddot {\boldsymbol{q}}_0} 分别为初始系统的广义坐标、广义速度和广义加速度. 当接触力{\tilde {\boldsymbol{Q}}_{\rm{c}}} 被检测为非零向量时, 执行细化子程序对接触区域施加进行网格局部细化操作.在细化子程序中, 当前时刻
t 的收敛结果{\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} ,\dot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} ,\ddot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} 为细化程序的输入变量. 首先, 模型的几何形状需要更新. 确定被插入结点\left( {\hat u,\hat v} \right) 后, 可得到转换矩阵{{\boldsymbol{M}}_{\rm{m}}} 以计算细化后的T样条控制点, 此时系统自由度由m 增加为n . 为了计算新系统的广义节点坐标、速度和加速度, 首先按照式(23)将初始系统的广义变量{\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} ,\dot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} 和\ddot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} 还原为初始系统的节点坐标{\boldsymbol{e}} , 节点速度\dot {\boldsymbol{e}} 和节点加速度\ddot {\boldsymbol{e}} \left. {\begin{array}{*{20}{l}} {{{\boldsymbol{e}}_{{\rm{temp}}}} = {{\boldsymbol{e}}_0} + {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} \cdot {\boldsymbol{B}}} \\ {{{\dot {\boldsymbol{e}}}_{{\rm{temp}}}} = {{\dot {\boldsymbol{e}}}_0} + \dot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} \cdot {\boldsymbol{B}}} \\ {{{\ddot {\boldsymbol{e}}}_{{\rm{temp}}}} = {{\ddot {\boldsymbol{e}}}_0} + \ddot {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} \cdot {\boldsymbol{B}}} \end{array}} \right\} (23) 然后根据转换矩阵
{{\boldsymbol{M}}_{\rm{m}}} 得到细化后新系统的节点坐标\tilde {\boldsymbol{e}} , 节点速度{\dot {\tilde {\boldsymbol{e}}}} 和节点加速度{\ddot {\tilde {\boldsymbol{e}}}} , 如下式(24)所示\left. {\begin{array}{*{20}{l}} {{{\tilde {\boldsymbol{e}}}_{{\rm{temp}}}} = \frac{{{{\boldsymbol{M}}_{\rm{m}}} \cdot {{\boldsymbol{e}}_{{\rm{temp}}}}}}{{\tilde {\boldsymbol{\omega}} }}} \\ {{{\dot {\tilde {\boldsymbol{e}}}}_{{\rm{temp}}}} = \frac{{{{\boldsymbol{M}}_{\rm{m}}} \cdot {{\dot {\boldsymbol{e}}}_{{\rm{temp}}}}}}{{\tilde {\boldsymbol{\omega}} }}} \\ {{{\ddot {\tilde {\boldsymbol{e}}}}_{{\rm{temp}}}} = \frac{{{{\boldsymbol{M}}_{\rm{m}}} \cdot {{\ddot {\boldsymbol{e}}}_{{\rm{temp}}}}}}{{\tilde {\boldsymbol{\omega}} }}} \end{array}} \right\} (24) 同理可得新系统在
t = 0 时刻的节点变量{\tilde {\boldsymbol{e}}_0} ,{{\dot {\tilde {\boldsymbol{e}}}}_0} 和{{\ddot {\tilde {\boldsymbol{e}}}}_0} .更新几何模型后, 约束方程也会发生变化, 所以将会得到新的布尔矩阵
\tilde {\boldsymbol{B}} . 根据\tilde {\boldsymbol{B}} 可以得到新系统的广义坐标\tilde {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} 、广义速度{\dot {\tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} 和广义加速度{\ddot{ \tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} , 其计算公式如下\left. {\begin{array}{*{20}{l}} {\tilde {\boldsymbol{B}} \cdot \tilde {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} = {{\tilde {\boldsymbol{e}}}_{{\rm{temp}}}} - {{\tilde {\boldsymbol{e}}}_0}} \\ {\tilde {\boldsymbol{B}} \cdot {\dot {\tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} = {{\dot {\tilde {\boldsymbol{e}}}}_{{\rm{temp}}}} - {{\dot {\tilde {\boldsymbol{e}}}}_0}} \\ {\tilde {\boldsymbol{B}} \cdot {\ddot{ \tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} = {{\ddot {\tilde {\boldsymbol{e}}}}_{{\rm{temp}}}} - {{\ddot {\tilde {\boldsymbol{e}}}}_0}} \end{array}} \right\} (25) 获取了上述变量之后, 系统动力学方程中的变量将会得到更新. 通过与
\tilde {\boldsymbol{B}} 的乘积, 将会得到新系统的广义质量阵\tilde {\boldsymbol{M}} , 广义外力{\tilde {\boldsymbol{Q}}_{\rm{e}}} , 广义弹性力{\tilde {\boldsymbol{Q}}_{\rm{s}}} 和广义接触力{\tilde {\boldsymbol{Q}}_{\rm{c}}} . 将更新后的动力学方程输入到广义α算法中进行求解, 便可得到网格局部细化、自由度增加后的新系统下一时间步的广义变量\tilde {\boldsymbol{q}}_{n + 1}^{\left( {k + 1} \right)} ,{\dot {\tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} 和{\ddot{ \tilde {\boldsymbol{q}}}}_{n + 1}^{\left( {k + 1} \right)} , 实现了对变自由度系统动力学方程的求解.4. 数值算例
4.1 受端弯矩作用的悬臂薄板
图4为一末端受弯矩
M 的悬臂薄板. 该薄板长度L{\rm{ = }}12\;{\rm{m}} , 宽度b = 1\;{\rm{m}} , 厚度h{\rm{ = }}0.1\;{\rm{m}} , 弹性模量E{\rm{ = }}1.2\;{\rm{MPa}} , 泊松比\nu {\rm{ = }}0 , 惯性矩I = \dfrac{{b{h^3}}}{{12}} . 图5给出了不同端弯矩下的薄板构型图.根据文献[32]可知, 取最大弯矩为
{M_0}{\rm{ = }}\dfrac{{2{\rm{{\text{π}} }}EI}}{L} , 当末端所受弯矩M = {M_0} , 此时薄板将弯曲成一个半径为R = \dfrac{{EI}}{M} = \dfrac{L}{{2{\rm{{\text{π}} }}}} 的圆. 当末端所受弯矩为M = 0.75{M_0} ,0.5{M_0} 和0.25{M_0} 时, 薄板的端截面转角分别为1.5{\rm{{\text{π}} }} ,{\rm{{\text{π}} }} 和0.5{\rm{{\text{π}} }} .![]()
图 5 端弯矩作用下, 变形后悬臂梁的构形图Figure 5. Configuration of deformed cantilever under external moments4.2 受内压圆环
一平面圆环构件, 在其中心孔边缘受沿半径向外方向的均布载荷
P = 10\;{\rm{MPa}} . 根据对称性, 取圆环构件的四分之一并在两侧施加滑动边界约束, 如图6(a)所示. 圆环构件的内与外径分别为{R_1}{\rm{ = }}1\;{\rm{m}} 和{R_2} = 2 m, 厚度h = 0.1 m. 材料密度\;\rho\; {\rm{ = }}\;7800\; kg/m3, 杨氏模量E = 210\; GPa, 泊松比\nu = 0.3.![]()
图 6 四分之一圆环构件及其局部细化示意图Figure 6. Diagram of a quarter of circular thin plate and its local refinement diagram采用T样条曲面单元对构件离散, 单元网格分布图如图6(b)所示. 因为最大米塞斯应力出现在AB边界附近, 内部区域被局域细化. 模型共含有28个单元和72个控制点. 图7给出了局部细化后构件的米塞斯应力云图.
因为该算例属于平面问题, 所以轴向应力
{\sigma _{{z}}} = 0 . 根据拉美方程可以得到圆环在任意半径r处米塞斯应力{\sigma _{\rm{s}}} 的理论解\left. {\begin{array}{*{20}{l}} {{\sigma _{{\theta }}}{{ = }}\dfrac{{{{{P}}_0}{K^2}}}{{{K^2} - 1}}\left( {1 - \dfrac{{R_2^2}}{{{r^2}}}} \right)} \\ {{\sigma _{{r}}}{{ = }}\dfrac{{{{{P}}_0}{K^2}}}{{{K^2} - 1}}\left( {1 + \dfrac{{R_2^2}}{{{r^2}}}} \right)} \\ {{\sigma _{{s}}} = \sqrt {\dfrac{1}{2}\left[ {{{\left( {{\sigma _{{r}}} - {\sigma _{{\theta }}}} \right)}^2} + {{\left( {{\sigma _{{r}}} - {\sigma _{{z}}}} \right)}^2} + {{\left( {{\sigma _{{\theta }}} - {\sigma _{{z}}}} \right)}^2}} \right]} } \end{array}} \right\} (26) 式中,
{\sigma _{{r}}} 为径向应力,{\sigma _{\rm{\theta }}} 为环向应力, K为圆环外径与内径之比,K = {R_2}:{R_1} .在有限元软件ANSYS中采用SHELL63单元运行相同算例并对比构件中最大和最小应力结果如表1所示. 可以看出, 相较于ANSYS结果, IGA结果更加接近于理论解而且与理论解的误差小于1‰.
表 1 圆环构件最小和最大米塞斯应力Table 1. Minimum and maximum von-Mises stress on circular thin plate{\sigma _{\min }}/MPa {\sigma _{\max }}/MPa theoretical solution 14.53 26.67 ANSYS solution (200 elements) 14.45 26.37 IGA solution (28 elements) 14.52 26.68 4.3 柔性单摆
该算例测试了基于T样条的等几何薄板的动力学特性. 如图8所示, 一个柔性单摆铰接固定在A处. 柔性薄板单摆在重力加速度
g = 9.81\;{\rm{m}}/{{\rm{s}}^2} 的作用下坠落. 薄板的参数在表2中给出.表 2 材料参数Table 2. Material parametersParameter Value length a/m 2 width b/m 2 height h/m 0.01 density \rho /(kg∙m−3) 7850 young’s modulus E/MPa 2 poisson ratio \mu 0.3 图9给出了含有2 × 2个T样条单元的薄板在1 s内动能
{E_{\rm{k}}} 、重力势能{E_{\rm{g}}} 、弹性势能{E_{\rm{e}}} 及总能量{E_{\rm{t}}} 的变化情况. 可以看出, 基于T样条的等几何薄板满足能量守恒定律. 为了检验本文提出的方法的收敛性, 分别使用2 × 2, 4 × 4和8 × 8的T样条曲面单元对图8所示的单摆进行离散. 图10给出了单摆自由端C的z向位移曲线. 可以看出, 基于T样条的等几何薄板具有较好的收敛性. 图11给出了含有2 × 2个T样条单元的单摆在1 s内的连续构型.![]()
图 11 含有2 × 2个T样条单元的单摆的连续构型Figure 11. Configuration of the pendulum with 2 × 2 T-spline elements4.4 刚性球与柔性板的接触
本算例将T样条局部细化算法应用于柔体动力学分析中. 如图12所示, 一个刚性球从四边固定的薄板中心正上方自由下落. 接触发生后, 对薄板中心受冲击区域进行局部细化. 刚性球与板的参数表3中给出.
![]()
图 12 自由坠落刚性球与柔性薄板的碰撞Figure 12. The collision between a rigid ball and a flexible thin plate表 3 接触算例的参数Table 3. Parameters of the contact exampleParameter Name Value the flexible
thin platelength a/m 2 width b/m 2 height h/m 0.04 density \;{\rho _1}/(kg∙m−3) 7850 Young’s modulus E/MPa 20 Poisson ratio \mu 0.3 the rigid ball radius R/m 0.1 density \;{\rho _2}/(kg∙m−3) 1500 initial position {{\boldsymbol{P}}_{\rm{b}}} (1,1,0.2) 在本算例中, 采用罚值法来完成刚体球与柔性板的接触实现. 在每个T样条单元中, 刚体球与均匀检测点的距离用式(27)表示
\left. {\begin{array}{*{20}{l}} {{{\boldsymbol{r}}_{\rm{s}}} = {\boldsymbol{r}}\left( {u,v} \right) + {\boldsymbol{n}} \cdot \dfrac{h}{2}} \\ {p = \left| {{{\boldsymbol{r}}_{\rm{s}}} - {{\boldsymbol{r}}_{\rm{c}}}} \right| - R} \\ {{\boldsymbol{n}} = \dfrac{{{{\boldsymbol{r}}_u} \times {{\boldsymbol{r}}_v}}}{{\left| {{{\boldsymbol{r}}_u} \times {{\boldsymbol{r}}_v}} \right|}}} \end{array}} \right\} (27) 式中,
p 为嵌入深度,{\boldsymbol{r}}\left( {u,v} \right) 为薄板中性层上的接触检测点,{{\boldsymbol{r}}_{\rm{s}}} 为任意点检测点在薄板上表面所对应的物质点,{{\boldsymbol{r}}_{\rm{c}}} 为球心位置,h 为薄板厚度,R 为球心半径,{\boldsymbol{n}} 表示接触检测点的单位法向量.一旦小球与薄板发生接触, 开始计算薄板所受的接触力. 计算切向接触力时, 采用文献[33]中计及临界滑动速度
{v_{\rm{0}}} 的平滑化库伦摩擦模型. 接触探测点i 的法向接触力{\boldsymbol{F}}_{\rm{n}}^i 和切向接触力{\boldsymbol{F}}_{\rm{t}}^i 表达式分别为{\boldsymbol{F}}_{\rm{n}}^i = -\left( {k{p^{1.5}} + c{v_{\rm{p}}}\left| p \right|} \right) \cdot {\boldsymbol{n}} (28) \begin{split} & {\boldsymbol{F}}_{\rm{t}}^i = \left\{ {\begin{array}{*{20}{l}} {{\rm{ - }}\mu \left| {{\boldsymbol{F}}_{\rm{n}}^i} \right|{{\boldsymbol{v}}_{{\rm{et}}}}},&{\left| {{{\boldsymbol{v}}_{\rm{t}}}} \right| < {{{v}}_0}} \\ {{\rm{ - }}\dfrac{{\left| {{{\boldsymbol{v}}_{\rm{t}}}} \right|}}{{{{{v}}_0}}}\mu \left| {{\boldsymbol{F}}_{\rm{n}}^i} \right|{{\boldsymbol{v}}_{{\rm{et}}}}},&{\left| {{{\boldsymbol{v}}_{\rm{t}}}} \right| \geqslant {{{v}}_0}} \end{array}} \right.\\ & {{\boldsymbol{v}}_{{\rm{et}}}} = \dfrac{{{{\boldsymbol{v}}_{\rm{t}} }}}{{\left| {{{\boldsymbol{v}}_{\rm{t}}}} \right|}} \end{split} (29) 式中,
{{{v}}_{\rm{p}}} 为嵌入速度,k 和c 表示刚度系数与阻尼系数. 在式中,\mu 为薄板与球之间的摩擦系数,{{\boldsymbol{v}}_{\rm{t}}} 为切向接触速度,{{{v}}_0} 为假定的临界滑动速度,{{\boldsymbol{v}}_{{\rm{et}}}} 为单元切向接触速度.图13介绍了三种网格构型图. 网格1和网格2都对表面进行了均匀的网格细化, 其单元数量分别为14 × 14和10 × 10. 网格3的初始网格较为粗糙, 含有10 × 10个单元. 当接触发生后, 薄板的接触区域被局部细化, 单元数量由100变为112个. 三组薄板的详细信息见表4.
表 4 三组薄板的详细信息Table 4. The detailed information of three groupsGroup Number of
control pointNumber of
elementElement
size/m2mesh 1 289 196 1/7 × 1/7 (all region) mesh 2 169 100 1/5 × 1/5 (all region) mesh 3 185 112 1/10 × 1/10 (contact region) 图14给出了三组小球球心在1 s内的z向位移曲线. 图15和图16则为三组薄板的质心z向位移曲线与弹性能曲线及其局部放大图. 在
t = 0.216\;{\rm{s}} 时, 刚性球与柔性板首次接触, 对薄板进行如图13所示的3种不同网格细化策略. 三组球的质心垂直位移具有较好的一致性, 薄板质心的垂直位移曲线在1 s内也呈现相似的趋势. 从薄板的质心位移与弹性能计算结果来看, 允许局部网格细化的T样条等几何薄板可以达到与退化为NURBS曲面的均匀网格板相同的计算精度, 且相较于网格2, 局部细化的网格3所对应曲线的变化趋势更接近网格1.为了比较不同网格构型的薄板对计算效率的影响, 图17给出了三组薄板的仿真时间柱状图.
三组算例均是在配备Intel Core i5 CPU和8GB RAM的笔记本电脑上执行的. 由图17可知, 网格2与网格3用时接近, 约为网格1计算用时的一半. 值得注意的是, 由于在接触区域附近执行了网格细化, 使得求解收敛更快, 网格3虽然自由度数略多于网格2, 用时反而更少. 这也进一步体现了T样条单元网格局部细化的优势.
图18展示了带有112个T样条单元局部细化薄板在不同时刻的米塞斯应力云图, 分别为
t = 0.18\;{\rm{s}} 接触前的薄板发生最大变形,t = 0.218\;{\rm{s}} 发生接触后首次使用细化后模型的时刻,t = 0.548\;{\rm{s}} 接触后的发生最大变形的时刻以及仿真结束时刻t = 1\;{\rm{s}} .![]()
图 18 局部细化薄板的米塞斯应力云图Figure 18. von-Mises stress distribution of the locally refined thin plate5. 结论
本文提出了一种基于T样条曲面的等几何分析方法, 建立了使用T样条曲面单元离散的基尔霍夫薄板模型. 单元基函数和节点坐标分别为T样条基函数和控制点坐标, 无需网格划分, 既保证模型的几何精确, 又能在没有约束方程的情况下保证期望的连续性条件. 给出了基于T样条的薄板运动学模型、弹性力模型及其雅克比矩阵的计算方法. 为了体现T样条局部细化特性在减少单元数目、提高计算效率方面的优势, 提出了一种基于T样条单元的局部网格细化算法, 包括单元网格拓扑的改变和相应的新控制点坐标计算方法. 将该细化算法与广义α法相结合, 建立了带有网格局部细化的变自由度系统动力学方程的求解算法. 通过受端弯矩的悬臂薄板以及环形构件受内压的静力学算例证明了T样条单元弹性力模型的正确性. 柔性摆算例验证了T样条单元在动力学问题中的收敛性和机械能守恒特性. 最后, 为验证T样条局部细化特性在模拟接触碰撞等柔性体局部发生应变剧烈变化等问题中的优势, 建立了刚性球落在柔性板上的动力学实例并进行了仿真. 首先, 采用10 × 10单元网格对柔性板进行离散. 接触发生后, 对冲击区域进行局部细化, 得到112个单元的网格. 对10 × 10和14 × 14两种网格进行了仿真, 并将仿真结果作为基准. 可以看到, 112单元网格的计算结果与14 × 14网格具有良好的一致性, 但所消耗的时间只有其1/2. 以上算例证明了基于T样条的局部网格细化等几何薄板在柔性多体系统动力学分析中的应用价值.
-
![]()
图 1 等几何基尔霍夫薄板积分过程
Figure 1. The process of isogeometric Kirchhoff thin plate integration
![]()
图 5 端弯矩作用下, 变形后悬臂梁的构形图
Figure 5. Configuration of deformed cantilever under external moments
![]()
图 6 四分之一圆环构件及其局部细化示意图
Figure 6. Diagram of a quarter of circular thin plate and its local refinement diagram
![]()
图 11 含有2 × 2个T样条单元的单摆的连续构型
Figure 11. Configuration of the pendulum with 2 × 2 T-spline elements
![]()
图 12 自由坠落刚性球与柔性薄板的碰撞
Figure 12. The collision between a rigid ball and a flexible thin plate
![]()
图 18 局部细化薄板的米塞斯应力云图
Figure 18. von-Mises stress distribution of the locally refined thin plate
表 1 圆环构件最小和最大米塞斯应力
Table 1 Minimum and maximum von-Mises stress on circular thin plate
{\sigma _{\min }}/MPa {\sigma _{\max }}/MPa theoretical solution 14.53 26.67 ANSYS solution (200 elements) 14.45 26.37 IGA solution (28 elements) 14.52 26.68 
下载: 导出CSV
表 2 材料参数
Table 2 Material parameters
Parameter Value length a/m 2 width b/m 2 height h/m 0.01 density \rho /(kg∙m−3) 7850 young’s modulus E/MPa 2 poisson ratio \mu 0.3
下载: 导出CSV
表 3 接触算例的参数
Table 3 Parameters of the contact example
Parameter Name Value the flexible
thin platelength a/m 2 width b/m 2 height h/m 0.04 density \;{\rho _1}/(kg∙m−3) 7850 Young’s modulus E/MPa 20 Poisson ratio \mu 0.3 the rigid ball radius R/m 0.1 density \;{\rho _2}/(kg∙m−3) 1500 initial position {{\boldsymbol{P}}_{\rm{b}}} (1,1,0.2)
下载: 导出CSV
表 4 三组薄板的详细信息
Table 4 The detailed information of three groups
Group Number of
control pointNumber of
elementElement
size/m2mesh 1 289 196 1/7 × 1/7 (all region) mesh 2 169 100 1/5 × 1/5 (all region) mesh 3 185 112 1/10 × 1/10 (contact region)
下载: 导出CSV
-
[1] Zienkiewicz OC, Taylor RL. The Finite Element Method for Solid and Structural Mechanics, Seventh Edition. Stoneham/London: Butterworth-Heinemann, 2009: 1-20.
[2] 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 doi: 10.1016/j.cma.2004.10.008
[3] Shabana AA, Xu L. Rotation-based finite elements: reference-configuration geometry and motion description. Acta Mechanica Sinica, 2021, 37(1): 105-126 doi: 10.1007/s10409-020-01030-6
[4] Lan P, Shabana AA. Integration of B-spline geometry and ANCF finite element analysis. Nonlinear Dynamics, 2010, 61(1): 193-206
[5] Lan P, Shabana AA. Rational finite elements and flexible body dynamics. Journal of Vibration and Acoustics of the ASME, 2010, 132(4): 041007 doi: 10.1115/1.4000970
[6] Lan P, Yu Z, Du L, et al. Integration of non-uniform Rational B-splines geometry and rational absolute nodal coordinates formulation finite element analysis. Acta Mechanica Solida Sinica, 2014, 27(5): 486-495 doi: 10.1016/S0894-9166(14)60057-4
[7] Yu Z, Lan P, Lu N. A piecewise beam element based on absolute nodal coordinate formulation. Nonlinear Dynamics, 2014, 77(1-2): 1-15 doi: 10.1007/s11071-014-1248-x
[8] Mikkola A, Shabana AA, Sanchez-Rebollo C, et al. Comparison between ANCF and B-spline surfaces. Multibody System Dynamics, 2013, 30(2): 119-138 doi: 10.1007/s11044-013-9353-z
[9] Pappalardo CM, Yu Z, Zhang X, et al. Rational ANCF thin plate finite element. Journal of Computational and Nonlinear Dynamics, 2016, 11(5): 051009 doi: 10.1115/1.4032385
[10] Yu Z, Cui Y. A new ANCF solid-beam element: relationship with bézier volume and application on leaf spring modeling. Acta Mechanica Sinica, 2021 (in press)
[11] Ma L, Wei C, Ma C, et al. Modeling and verification of a RANCF fluid element based on cubic rational bezier volume. Journal of Computational and Nonlinear Dynamics, 2020, 15(4): 041005
[12] Hamed AM, Jayakumar P, Letherwood MD, et al. Ideal compliant joints and integration of computer aided design and analysis. Journal of Computational and Nonlinear Dynamics, 2015, 10(2): 021015 doi: 10.1115/1.4027999
[13] Mizuno Y, Sugiyama H. Sliding and nonsliding joint constraints of B-spline plate elements for integration with flexible multibody dynamics simulation. ASME Journal of Computational and Nonlinear Dynamics, 2014, 9(1): 011001 doi: 10.1115/1.4025277
[14] 李新康, 张继发, 郑耀. Kirchhoff-Love壳的等几何分析. 固体力学学报, 2014, 35(S1): 129-133 (Li Xinkang, Zhang Jifa, Zheng Yao. Isogeometric analysis of Kirchhoff-Love shells. Chinese Journal of Solid Mechanics, 2014, 35(S1): 129-133 (in Chinese) [15] 李新康, 张继发, 郑耀. Mindlin板的等几何分析. 固体力学学报, 2013, 33(S1): 198-203 (Li Xinkang, Zhang Jifa, Zheng Yao. Isogeometric analysis of mindlin plates. Chinese Journal of Solid Mechanics, 2013, 33(S1): 198-203 (in Chinese) [16] Thai TQ, Rabczuk T, Zhuang X. Isogeometric cohesive zone model for thin shell delamination analysis based on Kirchhoff-Love shell model. Frontiers of Structural and Civil Engineering, 2020, 14(2): 267-279 doi: 10.1007/s11709-019-0567-x
[17] 刘涛, 汪超, 刘庆运等. 基于等几何方法的压电功能梯度板动力学及主动振动控制分析. 工程力学, 2020, 37(12): 228-242 (Liu Tao, Wang Chao, Liu Qingyun, et al. Analysis for dynamic and active vibration control of piezoelectric functionally graded plates based on isogeometric method. Engineering Mechanics, 2020, 37(12): 228-242 (in Chinese) doi: 10.6052/j.issn.1000-4750.2020.04.0266 [18] Sederberg TW, Cardon DL, Finnigan GT, et al. T-spline simplification and local refinement. Acm Transactions on Graphics, 2004, 23(3): 276-283 doi: 10.1145/1015706.1015715
[19] Kasik DJ, Buxton W, Ferguson DR. Ten CAD challenges. IEEE Computer Graphics Applications, 2005, 25(2): 81-92 doi: 10.1109/MCG.2005.48
[20] Sederberg TW, Zheng J, Bakenov A, et al. T-splines and T-NURCCs. Acm Transactions on Graphics, 2003, 22(3): 477-484 doi: 10.1145/882262.882295
[21] Sederberg TW, Finnigan GT, Li X, et al. Watertight trimmed NURBS. Acm Transactions on Graphics, 2008, 27(3): 1-8
[22] Casquero H, Liu L, Zhang Y, et al. Arbitrary-degree T-splines for isogeometric analysis of fully nonlinear Kirchhoff–Love shells. Computer-Aided Design, 2017, 82: 140-153 doi: 10.1016/j.cad.2016.08.009
[23] Casquero H, Wei X, Toshniwal D, et al. Seamless integration of design and Kirchhoff–Love shell analysis using analysis-suitable unstructured T-splines. Computer Methods in Applied Mechanics and Engineering, 2020, 360: 112765
[24] Dimitri R, De Lorenzis L, Scott MA, et al. Isogeometric large deformation frictionless contact using T-splines. Computer Methods in Applied Mechanics and Engineering, 2014, 269: 394-414 doi: 10.1016/j.cma.2013.11.002
[25] Uhm TK, Youn SK. T-spline finite element method for the analysis of shell structures. International Journal for Numerical Methods in Engineering, 2010, 80(4): 507-536
[26] Bazilevs Y, Calo VM, Cottrell JA, et al. Isogeometric analysis using T-splines. Computer Methods in Applied Mechanics and Engineering, 2015, 199(5-8): 229-263
[27] 刘登学, 张友良, 刘高敏. 基于适合分析T样条的高阶数值流形方法. 力学学报, 2017, 49(1): 212-222 (Liu Dengxue, Zhang Youliang, Liu Gaomin. Higher-order numerical manifold method based on analysis-suitable T-spline. Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(1): 212-222 (in Chinese) doi: 10.6052/0459-1879-16-217 [28] Li X, Zheng J, Sederberg TW, et al. On linear independence of T-spline blending functions. Computer Aided Geometric Design, 2012, 29(1): 63-76 doi: 10.1016/j.cagd.2011.08.005
[29] Scott MA, Borden MJ, Verhoosel CV, et al. Isogeometric finite element data structures based on Bézier extraction of T-splines. International Journal for Numerical Methods in Engineering, 2011, 88(2): 126-156 doi: 10.1002/nme.3167
[30] 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 doi: 10.1115/1.4023487
[31] 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
[32] Sze KY, Liu XH, Lo SH. Popular benchmark problems for geometric nonlinear analysis of shells. Finite Elements in Analysis and Design, 2004, 40(11): 1551-1569 doi: 10.1016/j.finel.2003.11.001
[33] 兰朋, 崔雅琦, 於祖庆. 完备化ANCF薄板单元及在钢板弹簧动力学建模中的应用. 力学学报, 2018, 50(5): 1156-1167 (Lan Peng, Cui Yaqi, Yu Zuqing. The completed form of elastic model for ANCF thin plate element and its application on dynamic modeling of the leaf spring. Chinese Journal of Theoretical and Applied Mechanics, 2018, 50(5): 1156-1167 (in Chinese) doi: 10.6052/0459-1879-18-133
计量
- 文章访问数: 1312
- HTML全文浏览量: 434
- PDF下载量: 84
