汤惠颖
* , 张志娟
? , 刘铖
, * , 2) , 刘绍奎
? * 北京理工大学 宇航学院, 北京 100081? 北京空间飞行器总体设计部, 北京 100094LOCKING ALLEVIATION TECHNIQUES OF TWO TYPES OF BEAM ELEMENTS BASED ON THE LOCAL FRAME FORMULATION1) Tang Huiying
* , Zhang Zhijuan
? , Liu Cheng
, * , 2) , Liu Shaokui
? * School of Aerospace Engineering, Beijing Institute of Technology, Beijing 100081, China? Beijing Institute of Spacecraft System Engineering, Beijing 100094, China通讯作者: 2) 刘铖, 副教授, 主要研究方向: 柔性多体系统动力学建模与计算; 非线性有限元方法等. E-mail:
liucheng_bit@aliyun.com 收稿日期: 2020-08-7
接受日期: 2020-11-16
网络出版日期: 2021-02-07
基金资助: 1) 国家自然科学基金 .12072026 国家自然科学基金 .11832005 国家自然科学基金 .11672034 十三五民用航天项目资助.
Received: 2020-08-7
Accepted: 2020-11-16
Online: 2021-02-07
作者简介 About authors
摘要 对于大转动、大变形柔性体的刚柔耦合动力学问题,基于李群SE(3)局部标架(local frame formulation, LFF)的建模方法能够规避刚体运动带来的几何非线性问题,离散数值模型中广义质量矩阵与切线刚度矩阵满足刚体变换的不变性,可明显地提高柔性多体系统动力学问题的计算效率. 有限元方法中,闭锁问题是导致单元收敛性能低下的主要原因, 例如梁单元的剪切以及泊松闭锁.多变量变分原理是缓解梁、板/壳单元闭锁的有效手段. 该方法不仅离散位移场,同时离散应力场或应变场, 可提高应力与应变的计算精度. 本文基于上述局部标架,研究几类梁单元的闭锁处理方法, 包括几何精确梁(geometrically exact beam formulation, GEBF)与绝对节点坐标(absolute nodal coordinate formulation, ANCF)梁单元. 其中, 采用Hu-Washizu三场变分原理缓解几何精确梁单元中的剪切闭锁,采用应变分解法缓解基于局部标架的ANCF全参数梁单元中的泊松闭锁. 数值算例表明,局部标架的梁单元在描述高转速或大变形柔性多体系统时,可消除刚体运动带来的几何非线性, 极大地减少系统质量矩阵和刚度矩阵的更新次数.缓解闭锁后的几类局部标架梁单元收敛性均得到了明显提升.
关键词： 局部标架方法 ;
几何精确梁 ;
绝对节点坐标梁单元 ;
闭锁 ;
Hu-Washizu三场变分原理 Abstract For rigid-flexible coupling dynamic problems with large rotation and large deformation, the modeling method based on the local frame formulation (LFF) of SE(3) group can avoid geometrically nonlinear problem caused by the rigid-body motion. In discretized flexible multibody systems, the generalized mass matrix and the tangent stiffness matrix are invariant under the arbitrary rigid-body motion, which can improve computational efficiency significantly. In the finite element method, locking is the main reason for low convergence rate of elements, such as shear and Poisson locking in beam elements. Mixed methods are effective strategies to alleviate locking in beam and plate/shell elements. In these methods, not only the displacement field but also the stress field and the strain field are discretized, which can increase the accuracy of stress and strain. Based on the local frame formulation, the paper studies locking alleviation techniques of several beam elements, including geometrically exact beam formulation (GEBF) and absolute nodal coordinate formulation (ANCF) beam elements. The Hu-Washizu variational principle is used to alleviate shear locking in the geometrically exact beam, while the strain split method is used to eliminate Poisson locking in the fully parameterized ANCF beam. Numerical examples show that the proposed beam elements based on the local frame formulation can eliminate geometrically nonlinearity caused by the rigid-body motion and can minimize the updating times of mass matrices and tangent stiffness matrices when modeling flexible multibody systems with high rotational speed or large deformation. After locking alleviation, the convergence rate of the above beam elements improves significantly.
Keywords： local frame formulation ;
geometrically exact beam formulation ;
absolute nodal coordinate beam formulation ;
locking ;
Hu-Washizu variational principle PDF (1001KB) 元数据 多维度评价 相关文章 导出 EndNote |
Ris |
Bibtex 收藏本文 本文引用格式 汤惠颖, 张志娟, 刘铖, 刘绍奎. 两类基于局部标架梁单元的闭锁缓解方法
1) .
力学学报 [J], 2021, 53(2): 482-495 DOI:
10.6052/0459-1879-20-274 Tang Huiying, Zhang Zhijuan, Liu Cheng, Liu Shaokui.
LOCKING ALLEVIATION TECHNIQUES OF TWO TYPES OF BEAM ELEMENTS BASED ON THE LOCAL FRAME FORMULATION1) .
Chinese Journal of Theoretical and Applied Mechanics [J], 2021, 53(2): 482-495 DOI:
10.6052/0459-1879-20-274 引言 几何精确方法
[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$; $ e_{1}$, $ e_{2}$, $ e_{3}$}. 在当前构型中, 在梁中线的每一点定义随体基矢量($ t$, $ n$, $ b$), 引入梁中线位置矢量$ x$和截面旋转矩阵$ R$, 初始构型对应的量为($\cdot$)$^{0}$.
图1 新窗口打开 |
下载原图ZIP |
生成PPT 图1几何精确梁初始构型与当前构型 Fig.1Initial and current configurations of the geometrically exact beam 在梁的当前构型中, 截面上任意一点$p$的位置矢量为
(1) $\begin{eqnarray} \label{eq1} { x}_{p} (\xi_{1} ,\xi_{2} ,\xi_{3} )={ x}(\xi_{1} )+{ R}(\xi_{1} ){ y}_{p} (\xi_{2} ,\xi_{3} ) \end{eqnarray}$ 式中, $\xi_{1}\in [0, L]$表示初始构型中$p$点所在截面对应的弧长坐标, $\xi _{2}$和$\xi_{3}$为$p$点在局部坐标系中的坐标分量, $L$为梁单元的长度, $ y_{p} = [0 \ \ \xi_{2} \ \ \xi_{3}]^{\rm T}$. 基于小应变假设, 当前构型的应变度量可以表示为
(2) $\begin{eqnarray} \label{eq2} { f}=\left[ {{\begin{array}{*{20}c} {{ f}_{u} }\\ {{ f}_{\omega } }\\ \end{array} }} \right]=\left[ {{\begin{array}{*{20}c} {{ R}^{\rm T}{ {x}'}}\\ {\mbox{axial}({ R}^{\rm T}{ {R}'})}\\ \end{array} }} \right] \end{eqnarray}$ 式中, $ f_{u}$和$ f_{\omega}$分别为应变度量$ f$的平动和转动部分, 记号$(\cdot )' = \partial (\cdot )/\partial \xi_{1} $, axial$(\cdot )$表示反对称矩阵$\cdot$的轴矢量. 初始构型的应变度量$ f^{0}$可以表示为
(3) $\begin{eqnarray} \label{eq3} { f}^{0}=\left[ {{\begin{array}{*{20}c} {{ R}^{\rm 0T}{ x}^{0'}}\\ {\mbox{axial}({ R}^{\rm 0T}{ R}^{0'})}\\ \end{array} }} \right] \end{eqnarray}$ 引入记号$\widehat{(\cdot )}$和$\check{(\cdot )}$, 它们的定义如下
(4) $\begin{eqnarray} \label{eq4} {\hat{{ f}}}=\left[ {{\begin{array}{c@{\quad }c} {{\tilde{{ f}}}_{\omega } }& {{\tilde{{ f}}}_{u} } \\ {{\bf 0}_{3} }& {{\tilde{{ f}}}_{\omega } }\\ \end{array} }} \right], \ \ \check{ f}=\left[ {{\begin{array}{c@{\quad }c} {{\bf 0}_{3} }& {{\tilde{{ f}}}_{u} }\\ {{\tilde{{ f}}}_{u} }& {{\tilde{{ f}}}_{\omega } }\\ \end{array} }} \right] \end{eqnarray}$ 式中, 记号$\tilde{{\cdot }}$表示矢量$\cdot $的反对称矩阵, ${\bf0}_{3}$表示3 $\times$ 3的零矩阵. 应变可以表示为
(5) $\begin{eqnarray} \label{eq5} { \varepsilon }={ f}-{ f}^{0}=\left[ {{\begin{array}{*{20}c} {{ \gamma }}\\ {{ \kappa }}\\ \end{array} }} \right]=\left[ {{\begin{array}{*{20}c} {{ R}^{\rm T}{ x'}-{ R}^{\rm 0T}{ x}^{0'}}\\ {\mbox{axial}({ R}^{\rm T}{ R'}-{ R}^{\rm 0T}{ R}^{0'})}\\ \end{array} }} \right] \end{eqnarray}$ 式中, $ \gamma = [\gamma_{1} \ \ \gamma_{2} \ \ \gamma_{3}]^{\rm T}$, $\gamma_{1}$表示梁中线的拉伸应变, $\gamma_{2}$和$\gamma_{3}$表示横向剪切应变, $ \kappa = [\kappa _{1} \ \ \kappa_{2} \ \ \kappa_{3}]^{\rm T}$, $\kappa_{1}$表示截面的扭转应变, $\kappa_{2}$和$\kappa_{3}$表示截面绕$ n$轴和$ b$轴的弯曲应变. 应力可以表示为
(6) $\begin{eqnarray} \label{eq6} { S}={ D}{ \varepsilon } \end{eqnarray}$ 式中, 弹性系数矩阵$ D = {\rm diag} (EA$, $k_{\rm s2}GA$, $k_{\rm s3}GA$, $GJ_{\rm s}$, $EI_{2}$, $EI_{3})$, 其中$E$和$G$分别为杨氏模量和剪切模量, $A$, $I_{2}$和$I_{3}$分别表示梁初始构型的截面面积和截面惯性矩, $J_{\rm s}$表示圣维南扭转常数, $k_{\rm s2}$和$k_{\rm s3}$为截面的剪切修正系数.
1.2 SE3几何精确梁控制方程及线性化 系统的动力学平衡方程为
(7) $\begin{eqnarray} \label{eq7} \delta {\cal W}_{\rm int} -\delta {\cal W}_{\rm ext} +\delta {\cal W}_{\rm ine} = {0} \end{eqnarray}$ 式中, $\delta {\cal W}_{\rm int} $为内力所做的虚功, $\delta {\cal W}_{\rm ext} $为外载荷所做的虚功, $\delta {\cal W}_{\rm ine} $为惯性力所做的虚功.
为获得一个对称的切线刚度矩阵, 本文在前一个收敛步所在切平面进行物理量的一次与二次变分. 因此, 应变的变分可以表示为
(8) $\begin{eqnarray} \label{eq8} \delta { \varepsilon }=\lt({\rm T}'+{\rm T}\frac{\rm d}{{\rm d}\xi_{1} }+{\hat{{ f}}{\rm T}})\delta { h} \triangleq {\rm B}\delta { h} \end{eqnarray}$ 式中, T为SE(3)群的切空间算符, $\delta h= [\delta d^{\rm T}\ \ \delta \eta^{\rm T}]^{\rm T}$, $\delta d$和$\delta \eta $分别为局部标架下的虚位移和虚转角, B为应变算符. 内部虚功可以表示为
(9) $\begin{eqnarray} \label{eq9} \delta {\cal W}_{\rm int} =\int_L {(\delta { \varepsilon })^{\rm T} D{ \varepsilon} {\rm d}\xi_{1} } =\int_L \delta h^{\rm T}{\rm B}^{\rm T} D{ \varepsilon }{\rm d}\xi _{1} \end{eqnarray}$ 外部虚功可以表示为
(10) $\begin{eqnarray} \label{eq10} \delta {\cal W}_{\rm ext} =\int_L {\delta { x}_{p}^{\rm T} { p}^{\rm ext}{\rm d}\xi_{1}} \end{eqnarray}$ 式中, $ p^{\rm ext}$为作用于$p$点的外载荷. 惯性虚功可以表示为
(11) $\begin{eqnarray} \label{eq11} \delta {\cal W}_{\rm ine} =\int {\delta{ h}^{\rm T}\left( {{ C\dot{{ v}}}-{\hat{{ v}}}^{\rm T}{ C v}} \right){\rm d}\xi_{1} }, \ \ { C}=\left[ {{\begin{array}{c@{\ \ }c} {\rho_{A} { I}_{3} }& {{\bf 0}_{3} }\\ {{\bf 0}_{3} }& {{ J}}\\ \end{array} }} \right]\quad \end{eqnarray}$ 式中, $ v= [ U^{\rm T}\ \ \omega^{\rm T}]^{\rm T}$, $ U$和 $ \omega $分别为局部标架下的线速度与角速度, 顶标 $\cdot$ 表示变量对时间$t$的全导数, $\rho_{A}$为单位长度的材料密度, $ I_{3}$表示3 $\times$ 3的单位矩阵, $ J$为截面形心惯性矩阵.
对内部虚功做线性化可得
(12) $\begin{eqnarray} \label{eq12} &&\Delta (\delta {\cal W}_{\rm int} )=\int_L {\Big(\delta { h}^{\rm T}{\rm B}^{\rm T}{ D}\Delta { \varepsilon }+\delta { h}^{\rm T}\Delta {\rm B}^{\rm T}{ S}\Big){\rm d}\xi_{1} } =\\&&\qquad \int_L {\Big(\delta { h}^{\rm T}{\rm B}^{\rm T}{ D{\rm B}}\Delta { h}+\delta { h}^{\rm T}\Delta {\rm B}^{\rm T}{ S}\Big){\rm d}\xi_{1} } \end{eqnarray}$ 应变算符B的线性化的转置为
(13) $\begin{eqnarray} \label{eq13} \Delta {\rm B}^{\rm T}=\Delta {\rm T'}^{\rm T}+\frac{\rm d}{{\rm d}\xi_{1} }\Delta {\rm T}^{\rm T}+\Delta {\rm T}^{\rm T}{\hat{{ f}}}^{\rm T}+{\rm T}^{\rm T}\Delta {\hat{{ f}}}^{\rm T} \end{eqnarray}$ 为了进一步简化$\Delta (\delta {\cal W}_{\rm int} )$的几何部分, 引入两个矩阵$\varXi_{{\rm T}}$和$\varXi_{\rm T'}$, 它们的定义如下
(14) $\left.\begin{array}{l}\Delta \boldsymbol{T}^{\mathrm{T}} \boldsymbol{A}=\boldsymbol{\Xi}_{\mathrm{T}}(\boldsymbol{A}) \Delta \boldsymbol{h} \\\Delta \boldsymbol{T}^{\prime \mathrm{T}} \boldsymbol{A}=\boldsymbol{\Xi}_{\mathrm{T}^{\prime}}(\boldsymbol{A})\left[\begin{array}{c}\Delta \boldsymbol{h} \\\Delta \boldsymbol{h}^{\prime}\end{array}\right]=\left[\begin{array}{lll}\boldsymbol{\Xi}_{\mathrm{T}^{\prime}}^{(1)}(\boldsymbol{A}) & \boldsymbol{\Xi}_{\mathrm{T}^{\prime}}^{(2)}(\boldsymbol{A})\end{array}\right]\left[\begin{array}{l}\Delta \boldsymbol{h} \\\Delta \boldsymbol{h}^{\prime}\end{array}\right]\end{array}\right\}$ 式中, $ A$为一个任意的6 $\times$ 1的矢量, 记号${\Xi }_{\rm T'}^{(1)} $和${\Xi }_{\rm T'}^{(2)} $分别是矩阵${\varXi }_{\rm T'} $的前半部分和后半部分, 其具体形式可参见文献[
11 ]. 因此, 几何部分为
(15) $\begin{eqnarray} \label{eq15} \Delta {\rm B}^{\rm T}{ S}=\varUpsilon ({ S})\Delta { h} \end{eqnarray}$ 式中
(16) $\begin{eqnarray} \label{eq16} &&\varUpsilon ({ S})={\varXi }_{\rm T'}^{(1)} ({ S})+{\varXi }_{{\rm T}} \Big( {\hat{{ f}}}^{\rm T}{ S}\Big)+{\rm T}^{\rm T}\check{ S}\Big({ {T}'}+{\hat{{ f}}T}\Big) +\\&&\qquad\Big[{\varXi }_{\rm T'}^{(2)} ({ S})+{\rm T}^{\rm T}\check{ S} {\rm T}\Big]\frac{\rm d}{{\rm d}\xi_{1} }+\frac{\rm d}{{\rm d}\xi _{1} }{\varXi }_{{\rm T}} ({ S}) \end{eqnarray}$ 最后, 内部虚功的线性化可以写为
(17) $\begin{eqnarray} \label{eq17} \Delta (\delta {\cal W}_{\rm int} )=\delta { h}^{\rm T}\int_L {\Big[{\rm B}^{\rm T}{ D{\rm B}}+\varUpsilon ({ S})\Big]{\rm d}\xi_{1} } \Delta { h} \end{eqnarray}$ 需要指出的是, 由于$\delta h$为局部标架下的平动与转动虚位移, $\Delta h$为局部标架下的平动与转动无限小位移, 它们均与该点的任意刚体位移无关. 因此, 单元的弹性力及其切线刚度矩阵均满足刚体变换的不变性. 对于广义惯性力及其广义质量矩阵也同时具有类似的形式.
1.3 SE3几何精确梁有限元离散格式 对梁中线位置矢量$ x$和截面相对旋转伪矢量$\varTheta^{r}$
[11 ] 进行空间离散
(18) $\begin{eqnarray} \label{eq18} { x}=N_{I} { x}_{I},\ \ \varTheta ^{r}=N_{I} \lg ({ R}_{r}^{-1} { R}_{I} ) \end{eqnarray}$ 式中, $N_{I}$是节点$I$的形函数, $ x_{I}$是梁中线上节点$I$的位置矢量, $ R_{r}$是参考点所在截面的旋转矩阵, $ R_{I}$是节点$I$所在截面的旋转矩阵, 旋转伪矢量为转轴单位矢量与转角的乘积.
离散应变算符B得到应变$\!-\!$位移矩阵$ B_{I}$, 它的表达式为
(19) $\begin{eqnarray} \label{eq19} { B}_{I} ={\rm T'}N_{I} +{\rm T}{N}'_{I} +{\hat{{ f}}\rm T}N_{I} \end{eqnarray}$ 弹性力可以表示为
(20) $\begin{eqnarray} \label{eq20} { F}_{\rm int} =\int_L {{ B}_{I}^{\rm T} { S}{\rm d}\xi_{1} } \end{eqnarray}$ 材料刚度和几何刚度矩阵可以表示为
(21) $\begin{eqnarray} \label{eq21} { K}_{\rm mat} =\int_L {{ B}_{I}^{\rm T} { D{\rm B}}_{I} {\rm d}\xi_{1}}, \ \ { K}_{\rm geo} =\int_L {N_{I}^{\rm T} \varUpsilon ({ S})N_{I} {\rm d}\xi_{1} } \end{eqnarray}$ 对于非常细长的梁结构, 采用一阶插值描述几何精确梁弯曲时将会产生剪切闭锁. 由于梁非常细长, 截面内横向剪切应变处处接近于0. 而由于平动参数$ x$和转动参数 $ \theta $采用了相同的插值函数, 使式(5)中横向剪切应变$ \gamma_{s}$中的$ R$ ($\theta)$和$ x'$两项具有不同阶次, $ \gamma_{s} = {\bf 0}$ 不可能处处满足, 导致能量泛函中剪切变形能项的量级不正确, 由此带来了剪切闭锁. 此时, 单元表现得十分刚硬, 即计算出来的位移值远小于参考解, 总体刚度矩阵很大.
2 基于局部标架的绝对节点坐标方法的三维全参数梁理论 本节简要介绍基于局部标架的绝对节点坐标全参数梁单元建模方法. 某质点从初始位置$ X$运动到当前位置$ r$, 变形梯度可以表示为
(22) $\begin{eqnarray} \label{eq22} { J}=\frac{\partial { r}}{\partial { X}} \end{eqnarray}$ Green-Lagrange应变张量可以表示为
(23) $\begin{eqnarray} \label{eq23} { \varepsilon }=\frac{1}{2}({ J}^{\rm T}{ J}-{ I}) \end{eqnarray}$ 应变能可以表示为
(24) $\begin{eqnarray} \label{eq24} U=\frac{1}{2}\int_{{V}} {{ \varepsilon }^{\rm T} D{ \varepsilon}{\rm d}V} \end{eqnarray}$ 式中, $ D$为弹性矩阵, $V$为梁的体积. 应变能对广义坐标$ e$求偏导得弹性力为
(25) $\begin{eqnarray} \label{eq25} { F}^{\rm int}=\frac{\partial U}{\partial { e}}=\int_V { \varepsilon ^{{\rm T}} D\frac{\partial \varepsilon }{\partial { e}}{\rm d}V} \end{eqnarray}$ 弹性力对广义坐标求偏导得刚度阵为
(26) $\begin{eqnarray} \label{eq26} { K}=\frac{\partial { F}^{\rm int}}{\partial { e}}=\int_V \left[ \left( \frac{\partial \varepsilon }{\partial { e}} \right)^{\rm T} D\frac{\partial \varepsilon }{\partial { e}}+ \varepsilon^{\rm T} D\frac{\partial^{2} \varepsilon }{\partial { e}^{2}} \right]{\rm d}V \end{eqnarray}$ 相比于几何精确梁单元, ANCF单元没有转动参数, 不能直接定义局部标架, 需要通过平动位移场定义转动场. 此处采用ANCF斜率矢量定义该点的局部标架基矢量
(27) $\left. \begin{array} \\ { R}=[{ i}_{1}\ \ { i}_{2}\ \ { i}_{3}], \ \ { 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\}$ 式中, $ r_{X}$和$ r_{Y}$为当前位置矢量$ r$对物质坐标的偏导数. 由于平动位移场在局部标架中描述, 单元的刚体运动可自然被消除. 因此, 系统质量矩阵以及切线刚度矩阵满足刚体运动的不变性, 从而可极大地减少系统质量矩阵与刚度矩阵在仿真过程中的更新次数.
ANCF全参数梁单元会遇到泊松闭锁问题, 闭锁产生的原因如下. 由于ANCF全参数梁是连续介质单元, 梁的中线为三阶Hermite插值, 在泊松效应的作用下, 截面应该变形为曲面. 然而横截面为一阶插值, 不能描述截面的曲面变形模式, 因此导致泊松闭锁.
3 两类局部标架梁单元的闭锁缓解方法 3.1 基于Hu-Washizu三场变分原理的剪切闭锁缓解方法 Hu-Washizu三场变分原理不仅假设位移场, 同时假设应力和应变场, 且三类场变量的插值相互独立. 假设剪切应变$\bar{ \gamma}$和假设剪切应力$\bar{ \tau }$的分布为
(28) $\begin{eqnarray} \label{eq28} \bar{ \gamma}={ N}_{\gamma }{ \alpha }_{e} ,\ \ {\bar{{ \tau }}}={ N}_{\tau } {\beta }_{e} \end{eqnarray}$ 式中, $ N_{\gamma } = N_{\tau } =[1, \xi ,\cdots, \xi^{n-2}]$, $ \alpha_{e} = [\alpha_{1}, \alpha_{2}$, $\cdots$, $\alpha_{n-1}]^{\rm T}$, $ \beta_{e}= [\beta_{1}, \beta_{2}, ..., \beta_{n-1}]^{\rm T}$, $n$为单元节点数. 剪切部分的泛函可以表示为
(29) $\begin{eqnarray} \label{eq29} \varPi_{s} =\int_L \left[ \frac{1}{2}\bar{{ \gamma }}_{s}^{\rm T} D_{\rm s} \bar{ \gamma}_{\rm s} +\bar{ \tau }^{\rm T}( \gamma _{\rm s} -\bar{ \gamma}_{\rm s}) \right]{\rm d}\xi_{1} \end{eqnarray}$ 式中, $ \gamma_{\rm s}$代表由等参插值产生的剪切应变, 弹性系数矩阵$ D_{\rm s} = {\rm diag} (k_{\rm s2}GA$, $k_{\rm s3}GA$). 对剪切部分的泛函$\varPi _{\rm s}$求变分得
(30) $\begin{eqnarray} \label{eq30} &&\hspace{-2mm} \delta \varPi_{\rm s}\!=\!\int_L {\left[ {{\bar{{ \gamma }}}_{\rm s}^{\rm T} D_{\rm s} \delta {\bar{{ \gamma }}}_{\rm s} +\delta {\bar{{ \tau }}}^{\rm T}({ \gamma }_{\rm s} -{\bar{{ \gamma }}}_{\rm s} )+\delta ({ \gamma }_{\rm s}^{\rm T} -{\bar{{ \gamma }}}_{\rm s}^{\rm T} ){\bar{{ \tau }}}} \right]{\rm d}\xi_{1} }=\\&&\qquad\hspace{-4mm} \int_L {\left[ {\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} D_{\rm s} {\bar{{ \gamma }}}_{\rm s} +\delta {\bar{{ \tau }}}^{\rm T}({ \gamma }_{\rm s} -{\bar{{ \gamma }}}_{\rm s} )+(\delta h^{\rm T}{\rm B}^{\rm 23T}-\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} ){\bar{{ \tau}}}} \right]{\rm d}\xi_{1} }\\&&\qquad \end{eqnarray}$ 式中, B$^{23}$为应变算符B的第2和第3行. 对$\delta \varPi_{\rm s}$线性化得
(31) $\begin{eqnarray} \label{eq31} &&\Delta (\delta \varPi_{s} )=\int_L \Big[\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} D_{\rm s} \Delta {\bar{{ \gamma }}}_{\rm s} +\delta {\bar{{ \tau }}}^{\rm T}(\Delta { \gamma }_{\rm s} -\Delta {\bar{{ \gamma }}}_{\rm s} )+\\&&\qquad (\delta h^{\rm T}{\rm B}^{\rm 23T}-\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} )\Delta {\bar{{ \tau }}}+\delta h^{\rm T}\Delta {\rm B}^{\rm 23T}{\bar{{ \tau }}}\Big]{\rm d}\xi_{1} =\\&&\qquad\int_L \Big[\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} D_{\rm s} \Delta {\bar{{ \gamma }}}_{\rm s} +\delta {\bar{{ \tau }}}^{\rm T}{\rm B}^{23}\Delta h-\delta {\bar{{ \tau }}}^{\rm T}\Delta {\bar{{ \gamma }}}_{\rm s} +\\&&\qquad\delta h^{\rm T}{\rm B}^{\rm 23T}\Delta {\bar{{ \tau }}}-\delta {\bar{{ \gamma }}}_{\rm s}^{\rm T} \Delta {\bar{{ \tau }}}+\delta h^{\rm T}\varUpsilon ( S_{\rm s} )\Delta h \Big]{\rm d}\xi_{1} \end{eqnarray}$ 式中, $ S_{\rm s}= [{0}, {\bar{{ \tau }}}, {0}, {0}, {0}]^{\rm T}$. 剪切部分的弹性力可以表示为
(32) $\begin{eqnarray} \label{eq32} { F}_{\rm s}^{\rm int} =\int_L {\left[ {\begin{array}{c} {\rm B}^{\rm 23T}\bar{{\tau }} \\ { N}_{\gamma }^{\rm T} ({ D}_{\rm s} \bar{{\gamma }}-\bar{{\tau }}) \\ { N}_{\tau }^{\rm T} (\gamma_{\rm s} -\bar{{\gamma }}_{\rm s} ) \\ \end{array}} \right]{\rm d}\xi_{1} } \end{eqnarray}$ 剪切部分的刚度矩阵可以表示为
(33) $\begin{eqnarray} \label{eq33} { K}_{\rm s} =\left[ {{\begin{array}{c@{\ \ }c@{\ \ }c} {{ K}_{\rm s}^{\rm geo} }& {\bf 0}& {{ A}^{\rm T}} \\ {\bf 0}& {{ W}}& {-{ C}^{\rm T}}\\ {{ A}}& {-{ C}}& {\bf 0}\\ \end{array} }} \right] \end{eqnarray}$ 式中,
(34) $\left. \begin{array} \\ { K}_{\rm s}^{\rm geo} =\int_L {N_{I}^{\rm T} \varUpsilon ({ S}_{\rm s} )N_{I} {\rm d}\xi_{1} } , \ \ { A}=\int_L {{ N}_{\tau }^{\rm T}{\rm B}^{23}{\rm d}\xi_{1} } \\ { C}=\int_L {{ N}_{\tau }^{\rm T} { N}_{\gamma } {\rm d}\xi_{1}}, \ \ { W}=\int_L {{ N}_{\gamma }^{\rm T} { D}_{\rm s}{ N}_{\gamma } {\rm d}\xi_{1} } \\ \end{array} \right\}$ 离散方程组可以表示为
(35) $\begin{eqnarray} \label{eq35} \left[ {{\begin{array}{ccc} {{ K}_{\rm m} +{ K}_{\rm b} +{ K}_{\rm s}^{\rm geo} }& {\bf 0}& {{ A}^{\rm T}}\\ {\bf 0}& {{ W}}& {-{ C}^{\rm T}}\\ {{ A}}& {-{ C}}& {\bf 0}\\ \end{array} }} \right]\left[ {\begin{array}{l} \Delta { h}^{N} \\ \Delta \alpha_{e} \\ \Delta \beta_{e} \\ \end{array}} \right]=-{ R e s} \end{eqnarray}$ 式中, $ K_{\rm m}$和$ K_{\rm b}$分别为薄膜和弯曲部分的刚度矩阵, $\Delta h^{N}$为节点运动矢量的增量, $ R e s$为残差. 将上式静力缩聚化为单场形式
(36) $\begin{eqnarray} \label{eq36} \Big({ K}_{\rm m} +{ K}_{\rm b} +{ K}_{\rm s}^{\rm geo} +{ A}^{\rm T}{ C}^{\rm -T}{ W C}^{-1}{ A}\Big)\Delta { h}^{N}=-{ R e s} \end{eqnarray}$ 一阶插值的几何精确梁会产生严重的剪切闭锁, 二阶及二阶以上插值的几何精确梁剪切闭锁较轻. 因此, 本文关于几何精确梁剪切闭锁处理方法着重于一阶梁单元. Hu-Washizu三场变分原理中, 位移场、应变场和应力场均独立插值. 式(28)假设横向剪切应变为常数, 可避免产生过大的剪切应变能, 有效缓解了剪切闭锁.
此外, 对于梁单元, 缩减积分技术也是缓解剪切闭锁的有效手段. 数值算例展示了采用Hu-Washizu三场变分原理与缩减积分技术缓解梁单元剪切闭锁的效果. 相比较而言, 缩减积分技术由于其简便性在梁单元中应用较为广泛, 但用在板壳单元中会产生沙漏变形模式. 进一步, 上述Hu-Washizu三场变分原理的研究将为低阶几何精确板壳单元的闭锁缓解起到一定的借鉴作用.
3.2 基于应变分解法的泊松闭锁缓解方法 ANCF全参数梁单元的位置场和变形梯度可以写为
(37) $\begin{eqnarray} \label{eq37} \left. {\begin{array}{l} { r}={ r}^{\rm c}+Y{ r}_{Y} +Z{ r}_{Z} \\ { r}_{X} ={ r}_{X}^{\rm c} +Y{ r}_{YX} +Z{ r}_{ZX} \\ { r}_{Y} ={ r}_{Y} , \ \ { r}_{Z} ={ r}_{Z} \\ \end{array}} \right\} \end{eqnarray}$ 式中, 上标c代表中线. 于是, 位置矢量梯度矩阵可以写为
(38) $\begin{eqnarray} \label{eq38} { J}=[{ r}_{X}^{\rm c} +Y{ r}_{YX} +Z{ r}_{ZX}\ \ { r}_{Y}\ \ { r}_{Z} ]={ J}^{\rm c}(X)+{ J}^{k}(Y,Z) \end{eqnarray}$ 式中, ${ J}^{\rm c}=[{ r}_{X}^{\rm c}\ \ { r}_{Y}\ \ { r}_{Z} ]$仅与中线参数$X$相关, ${ J}^{k}=[Y{ r}_{YX}+Z{ r}_{ZX}\ \ {\bf 0}_{3\times 1}\ \ {\bf 0}_{3\times 1} ]$与单元的剪切、扭转和弯曲变形相关. 因此, Green-Lagrange应变张量可以写为
(39) $\begin{eqnarray} \label{eq39} &&{ \varepsilon }=\frac{1}{2}\Big({ J}^{\rm T}{ J}-{ I}\Big)=\frac{1}{2}\Big({ J}^{\rm cT}{ J}^{\rm c}+{ J}^{\rm cT}{ J}^{k}+{ J}^{k\rm T}{ J}^{\rm c}+\\&&\qquad{ J}^{k\rm T}{ J}^{k}-{ I}\Big) \end{eqnarray}$ 将$ J^{\rm c}$和$ J^{k}$代入上式得Green-Lagrange应变张量的分量为
(40) $\left. \begin{array} \\ { \varepsilon }=\dfrac{1}{2}\left[ {{\begin{array}{c@{\ \ \ }c@{\ \ \ }c} {\varepsilon_{11} }&& {\rm sym}\\ {{ r}_{X}^{\rm cT} { r}_{Y} +Y{ r}_{YX}^{\rm T} { r}_{Y} +Z{ r}_{ZX}^{\rm T} { r}_{Y} }& {{ r}_{Y}^{\rm T} { r}_{Y}-1}&\\ {{ r}_{X}^{\rm cT} { r}_{Z} +Y{ r}_{YX}^{\rm T} { r}_{Z} +Z{ r}_{ZX}^{\rm T} { r}_{Z} }& {{ r}_{Y}^{\rm T} { r}_{Z} }& {{ r}_{Z}^{\rm T} { r}_{Z} -1}\\ \end{array} }} \right] \\ \varepsilon_{11} ={ r}_{X}^{\rm cT} { r}_{X}^{\rm c}+2Y{ r}_{X}^{\rm cT} { r}_{YX} +2Z{ r}_{X}^{\rm cT} { r}_{ZX} +2YZ{ r}_{YX}^{\rm T} { r}_{ZX}+\\\qquad Y^{2}{ r}_{YX}^{\rm T} { r}_{YX}+Z^{2}{ r}_{ZX}^{\rm T} { r}_{ZX} -1 \\ \end{array} \right\}$ Green-Lagrange应变张量可以分解为两部分$ \varepsilon = \varepsilon^{\rm c} + \varepsilon^{k}$
[30 ] , 式中,
(41) $\left. \begin{array} \\ { \varepsilon}^{\rm c}\!=\!\dfrac{1}{2}\Big({ J}^{\rm cT}{ J}^{\rm c}-{ I}\Big)\!=\!\dfrac{1}{2}\left[ {{\begin{array}{ccc} {{ r}_{X}^{\rm cT} { r}_{X}^{\rm c} -1}& {{ r}_{X}^{\rm cT} { r}_{Y} }& {{ r}_{X}^{\rm cT} { r}_{Z} }\\ & {{ r}_{Y}^{\rm T} { r}_{Y} -1}& {{ r}_{Y}^{\rm T} { r}_{Z} }\\ {\rm sym}&& {{ r}_{Z}^{\rm T} { r}_{Z} -1}\\ \end{array} }} \right] \\ { \varepsilon}^{k}=\dfrac{1}{2}\Big({ J}^{\rm cT}{ J}^{k}+{ J}^{k\rm T}{ J}^{\rm c}+{ J}^{k\rm T}{ J}^{k}\Big) =\\\qquad\dfrac{1}{2}\left[ {{\begin{array}{c@{\ \ }c@{\ \ }c} {\varepsilon_{11}^{k}}& {Y{ r}_{YX}^{\rm T} { r}_{Y} +Z{ r}_{ZX}^{\rm T} { r}_{Y} }& {Y{ r}_{YX}^{\rm T} { r}_{Z}+Z{ r}_{ZX}^{\rm T} { r}_{Z} }\\ & 0& 0\\ {\rm sym}&& 0\\ \end{array} }} \right] \\ \varepsilon_{11}^{k} =2Y{ r}_{X}^{\rm cT} { r}_{YX} +2Z{ r}_{X}^{\rm cT} { r}_{ZX} +2YZ{ r}_{YX}^{\rm T} { r}_{ZX} +\\\qquad Y^{2}{ r}_{YX}^{\rm T} { r}_{YX} +Z^{2}{ r}_{ZX}^{\rm T} { r}_{ZX} \\ \end{array} \right\}$ $ \varepsilon^{\rm c}$中只有低阶应变, 且$ \varepsilon^{\rm c}$仅与梁中线参数$X$相关; $ \varepsilon^{k}$中包含高阶应变, 这些高阶应变会导致截面变形和弯曲变形. 为了缓解闭锁, 泊松耦合只用在正应变的低阶项中, 即仅在$ \varepsilon^{\rm c}$中考虑泊松效应, 忽略正应变的高阶项的泊松效应. 因此, 梁截面不会变形为曲面, 缓解了泊松闭锁. 第二类Piola-Kirchhoff应力的Voigt形式可以写为
(42) $\begin{eqnarray} \label{eq42} { \sigma }_{v} ={ D}^{\rm c}{ \varepsilon }_{v}^{\rm c} +{ D}^{k}{ \varepsilon }_{v}^{k} \end{eqnarray}$ 式中, 应变也是Voigt形式. 基于平面应变假设, 弹性系数矩阵可以被定义为
(43) $\left. \begin{array} \\ { D}^{\rm c}=\left[ {{\begin{array}{*{20}c} {\lambda +2\mu }& \lambda& \lambda& 0& 0 & 0\\ & {\lambda +2\mu }& \lambda& 0& 0& 0 \\ && {\lambda +2\mu }& 0& 0& 0\\ &&& \mu& 0& 0\\ & {\rm sym}&&& {\mu k_{\rm s2} }& 0\\ &&&&& {\mu k_{\rm s3} }\\ \end{array} }} \right] \\ { D}^{k}=\mbox{diag}(E,E,E,\mu ,\mu k_{\rm s2} ,\mu k_{\rm s3} ) \\ \end{array} \right\}$ 式中, $\lambda=E\nu /[(1+\nu )(1-2\nu )]$和$\mu =E/[2(1+\nu )]$为拉梅常数, $E$和$\nu $分别为杨氏模量和泊松比. 应变能可以表示为
(44) $\begin{eqnarray} \label{eq44} U=\frac{1}{2}\int_V {\left( {{ \varepsilon }_{v}^{\rm cT} { D}^{\rm c}{ \varepsilon }_{v}^{\rm c} +{ \varepsilon }_{v}^{k\rm T} D^{k}{ \varepsilon }_{v}^{k} } \right)} {\rm d}{V} \end{eqnarray}$ 弹性力可以表示为
(45) $\begin{eqnarray} \label{eq45} F^{\rm int }=\int_{V}\left[\left(\frac{\partial \varepsilon _{v}^{\rm c}}{\partial e }\right)^{\rm T} D^{\rm c} \varepsilon _{v}^{\rm c}+\left(\frac{\partial \varepsilon _{v}^{k}}{\partial e }\right)^{\rm T} D^{k} \varepsilon _{v}^{k}\right]{\rm d} V \end{eqnarray}$ 材料和几何刚度可以表示为
(46) $\left. \begin{array} \\ K _{\rm mat}=\int_{V}\left[\left(\dfrac{\partial \varepsilon _{v}^{\rm c}}{\partial e }\right)^{\rm T} D ^{\rm c} \dfrac{\partial \varepsilon _{v}^{\rm c}}{\partial e }+\left(\dfrac{\partial \varepsilon _{v}^{k}}{\partial e }\right)^{\rm T } D ^{k} \dfrac{\partial \varepsilon _{v}^{k}}{\partial e }\right] {\rm d} V \\ K _{\rm geo}=\int_{V}\left( \sigma _{v}^{\rm cT} \dfrac{\partial^{2} \varepsilon _{v}^{\rm c}}{\partial e ^{2}}+ \sigma _{v}^{k\rm T } \dfrac{\partial^{2} \varepsilon _{v}^{k}}{\partial e ^{2}}\right) {\rm d} V \end{array} \right\}$ 4 三类常用的有限元插值方法 本节简要介绍三类有限元插值方法, 包括Lagrange插值、Hermite插值以及非均匀有理B-样条(NURBS). 其中, Lagrange插值在有限元中应用最为广泛; 绝对节点坐标方法则采用Hermite插值, 以保证位移场的C$^{1}$连续; NURBS的基函数可以构造任意阶连续的近似函数, 建立了计算机辅助设计(CAD)与计算机辅助工程(CAE)间的桥梁, 该方法被称为等几何分析.
为对比分析不同插值方法的优劣势, 本文涉及一阶与三阶Lagrange插值、三阶Hermite插值以及三阶NURBS插值.
4.1 Lagrange插值 给定函数$f(x)$在$n+1$个互不相同的点$x_{i}$ ($i=0$, 1, 2, $\cdots$, $n)$上的函数值$f(x_{i}) = y_{i}$, $n$次Lagrange插值多项式可以表示为
(47) $\begin{eqnarray} \label{eq47} P_{n} (x)=\sum\limits_{i=0}^n {\left( {\prod\limits_{j=1,j\ne i}^n {\frac{x-x_{j} }{x_{i} -x_{j} }} } \right)} y_{i} \end{eqnarray}$ 4.2 Hermite插值 给定函数$f(x)$在$n+1$个互不相同的点$x_{i}$ ($i = 0$, 1, 2, $\cdots$, $n)$上的函数值$y_{i}$及一阶导数值${y}'_{i} $, Hermite插值多项式可以表示为
(48) $\begin{eqnarray} \label{eq48} H_{2n+1} (x)=\sum\limits_{i=0}^n {\Big[y_{i} A_{i} (x)+{y}'_{i} B_{i} (x)\Big]} \end{eqnarray}$ Hermite插值基函数$A_{i}(x)$和$B_{i}(x)$为
(49) $\left. \begin{array} \\ A_{i} (x)=\Big[1-2(x-x_{i} ){l}'_{i} (x_{i} )\Big]l_{i}^{2} (x) \\ B_{i} (x)=(x-x_{i} )l_{i}^{2} (x) \\ \end{array} \right\}$ 式中, $l_{i} (x)=\prod\limits_{j=0,j\ne i}^n {(x-x_{j} )/(x_{i} -x_{j} )} $.
4.3 非均匀有理B-样条 在区间[$a$, $b$]内, 设$U =\{u_{0},u_{1}, u_{2}, \cdots, u_{m}\}$ 是一不减的实数序列. 以$U$作为样条结点, $p$次的第$i$个B-样条基函数$N_{i,p}(u)$定义为
(50) $\left. \begin{array} \\ N_{i,0} (u)=\left\{\begin{array}{l@{\quad }l} 1, & u_{i} \leqslant u<u_{i+1} \\ 0, & {\rm else}\\ \end{array} \right. \\ N_{i,p} (u)=\dfrac{u-u_{i} }{u_{i+p} -u_{i} }N_{i,p-1} (u)+\\ \qquad \dfrac{u_{i+p+1} -u}{u_{i+p+1} -u_{i+1} }N_{i+1,p-1} (u) \\ \end{array} \right\}$ $p$次NURBS曲线定义为
(51) $\begin{eqnarray} \label{eq51} { C}(u)=\frac{\sum\limits_{i=0}^n {N_{i,p} (u)\omega_{i} { P}_{i} } }{\sum\limits_{i=0}^n {N_{i,p} (u)\omega_{i} } }, \ \ a\leqslant u\leqslant b \end{eqnarray}$ 式中, $ P_{i}$是控制点, $\omega_{i}$是权系数.
5 直/曲梁结构大变形静力学分析 本节考察$XOY$平面内的直梁与曲梁模型, 验证上述闭锁处理方法的有效性. 其中, 悬臂直梁长$L = 1$ m; 悬臂曲梁半径$R = 1$ m, 初始构型为1/4圆. 梁的截面均为矩形, 宽和高为0.01 m; 材料的杨氏模量为$E = 2.1$ $\times$ 10$^{11}$ Pa, 泊松比为$\nu = 0.3$. 首先, 在直梁和曲梁末端分别施加集中力(30, 20, 400) N (直梁)、(30, 20, 100) N (曲梁), 用ABAQUS计算得到梁端点沿$z$方向的位移分别为0.5143 m (直梁)、0.6138 m (曲梁). 计算的过程中, 上述两类载荷均分为10个加载步均匀加载. 此外, 对直梁与曲梁进行纯弯曲测试. 在直梁和曲梁末端均施加$z$方向的力矩$M_{z}=-\pi EI/2L$, 由材料力学解答知, 梁端沿$z$方向转角位移的解析解为$\theta_{z}=-\pi /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$方向的位移$u_{z}$, 计算相对误差, 将结果列于
表1 和
表2 中. 同时, 采用ABAQUS软件中B31梁单元(该单元为三维二节点一阶梁单元, 每个节点6个广义坐标, 分别是3个平移和3个转动参数)对本算例进行仿真计算.
图2 新窗口打开 |
下载原图ZIP |
生成PPT 图2(a)直梁和(b)曲梁变形前后示意图 Fig.2Initial and deformed configurations of (a) the straight beam and (b) the curved beam Table 1 表1 表1 悬臂直梁末端$z$方向位移和相对误差
Table 1
End displacements of $z$ direction and relative error of the straight cantilever beam 新窗口打开 |
下载CSV Table 2 表2 表2 悬臂曲梁末端$z$方向位移和相对误差
Table 2
End displacements of $z$ direction and relative error of the curved cantilever beam 新窗口打开 |
下载CSV 表1 和
表2 的数值结果表明: 采用缩减积分技术的SE(3)几何精确梁与ABAQUS中B31单元收敛结果基本一致. 采用Hu-Washizu三场变分原理缓解闭锁的SE(3)几何精确梁与其他两类方法的收敛结果不同, 其中, 直梁模型收敛值偏大, 曲梁模型收敛值偏小. 在三维复杂应力状态下, 一阶直梁和曲梁均有剪切闭锁问题. 缓解闭锁后, 直梁的计算精度提高了98%以上, 曲梁的计算精度提高了52%以上, 效果十分显著. 对于一阶直梁, 以上两种缓解闭锁的方法计算精度无明显差别. 对于一阶曲梁, 缩减积分技术的计算精度略高于Hu-Washizu三场变分原理.
5.2 悬臂直梁和曲梁纯弯曲测试 本节考查SE(3)几何精确一阶直梁与曲梁单元在纯弯曲应力状态下, 缩减积分技术和Hu-Washizu三场变分原理缓解剪切闭锁的能力, 同时评价单元描述纯弯曲的能力.
在悬臂直梁和曲梁末端施加集中弯矩, 此时, 悬臂梁发生纯弯曲, 处于简单应力状态. 单元位移场采用一阶Lagrange插值. 对比分析完全积分、缩减积分和Hu-Washizu三场变分原理三种情况下, 悬臂梁末端点$z$方向的转角位移$\theta_{z}$, 计算相对误差, 将结果列于
表3 和
表4 中. 同时, 采用ABAQUS软件中B31梁单元对本算例进行仿真计算.
Table 3 表3 表3 悬臂直梁末端$z$方向角位移和相对误差
Table 3
Angular end displacements of $z$ direction and relative error of the straight cantilever beam 新窗口打开 |
下载CSV Table 4 表4 表4 悬臂曲梁末端$z$方向角位移和相对误差
Table 4
Angular displacements of $z$ direction and relative error of the curved cantilever beam 新窗口打开 |
下载CSV 表3 和
表4 的数值结果表明: 在纯弯曲应力状态下, 一阶直梁和曲梁均有剪切闭锁问题. 缓解闭锁后, 直梁和曲梁模型的计算精度均提高了98%以上, 效果十分显著. 而且, 以上两种缓解闭锁的方法和ABAQUS中B31单元的计算精度完全一致.
表3 中, 一阶直梁在任意单元数目下, 转角与解析解的相对误差均为零. 这是由于SE(3)几何精确直梁模型能够精确地构造一个常弯曲单元, 因此仅用1个单元就能够精确模拟纯弯曲状态下的转动场.
表4 中, 一阶曲梁模型的转角渐近收敛, 这是由于一阶插值函数无法精确地描述一段圆弧.
同时, 采用三阶Lagrange插值构造了高阶几何精确直梁与曲梁单元. 通过上述悬臂直梁与曲梁算例, 测试了单元收敛精度, 其中同样采用Hu-Washizu三场变分原理处理剪切闭锁. 数值结果表明: 处理闭锁不会提升单元收敛性能, 进一步说明了高阶直梁与曲梁单元不存在剪切闭锁. 在5.3节中, 将进一步对比处理闭锁后的一阶单元与高阶单元的收敛性.
5.3 几类局部标架梁单元收敛性对比分析 本节对比一阶和三阶SE(3)几何精确梁单元、一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁、ANCF缩减(slope deficiency)梁、ANCF全参数梁与ABAQUS软件中B31梁单元的收敛性. 其中, 一阶几何精确梁与ANCF全参数梁单元均已处理闭锁.
在悬臂直梁和曲梁末端施加集中力. 计算以上几类建模方法在不同自由度下悬臂梁末端点$z$方向的位移$u_{z}$和相对误差, 将结果列于
表5 和
表6 中. ANCF缩减曲梁单元无算进行本算例的计算。
Table 5 表5 表5 悬臂直梁末端$z$方向位移和相对误差
Table 5
End displacements of $z$ direction and relative error of straight cantilever beams 新窗口打开 |
下载CSV Table 6 表6 表6 悬臂曲梁末端$z$方向位移和相对误差
Table 6
End displacements of $z$ direction and relative error of curved cantilever beams 新窗口打开 |
下载CSV 表5 和
表6 的数值结果表明: 直梁模型的五种建模方法收敛结果基本一致. 对于曲梁模型, SE(3)几何精确梁、$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁和ABAQUS中B31单元收敛结果基本一致, ANCF全参数梁单元收敛结果略大于其他三类方法的收敛值. 需要指出的是,
表5 中三阶NURBS插值的SE(3)几何精确梁单元在24个自由度时, 虽然计算结果的相对误差为6.6 $\times$ 10$^{-4}$, 但此时单元还没有收敛, 60个自由度时, 单元才收敛. 对于高阶单元, 三阶Lagrange比三阶NURBS插值的SE(3)几何精确梁的计算精度高, 几类几何精确梁均比ANCF梁单元的计算精度高. ANCF缩减梁比ANCF全参数梁单元的计算精度高.
对于一阶直梁单元, 两类几何精确梁具有完全相同的计算精度. 对于一阶曲梁单元, $\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁的计算精度略高于SE(3)几何精确梁. 对于一阶直梁和曲梁单元, 两类几何精确方法均比ABAQUS中B31单元的计算精度高. 显然高阶单元比低阶单元精度高, 但是高阶单元计算效率低, 即自由度相同时, 高阶单元刚度矩阵内非零元素的个数多.
6 高转速/大变形柔性多体系统动力学验证算例 本节通过高转速曲柄滑块机构以及空间双摆的动力学仿真, 分析局部标架下的梁单元在动力学仿真中的计算精度差异, 以及消除刚体运动带来的几何非线性后梁单元的数值特性. 由于商业软件无法直接进行大变形柔性多体系统动力学仿真, 因此本节不再采用商业软件进行仿真对比. 为了设计特定的小变形或大变形的柔性多体系统动力学模型, 本节算例中部分材料参数未与真实材料对应.
本节数值算例共涉及如下3类建模方法: (1)一阶和三阶Lagrange插值的SE(3)几何精确梁单元; (2)一阶Lagrange插值的$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁单元; (3) ANCF全参数梁单元. 需要指出的是, 除了三阶SE(3)几何精确梁单元, 其他单元均已处理闭锁.
6.1 发生小变形和中等变形的高速旋转的平面曲柄滑块机构 本节对一个典型的多体系统————曲柄滑块机构进行动力学仿真. 如
图3 所示, 平面曲柄滑块机构由两根杆与一个滑块组成, 杆I一端与地面采用球铰相连接, 另一端与杆II一端也通过球铰连接, 杆II末端与一个仅能在$X$轴上运动的刚体滑块铰接.
图3 新窗口打开 |
下载原图ZIP |
生成PPT 图3曲柄滑块机构本 Fig.3A crank-slider mechanism 算例中, 通过杆II中点$C$到其首尾连线的距离$d$度量杆II的弹性变形, 变形率为$d$与杆原长之间的比值. 同时, 以1944个自由度的$\mathbb{R}^{3}$空间ANCF全参数梁单元所得到的数值结果作为参考值, 建模方法数值结果收敛标准设为变形场最大值与参考解的相对误差在1%以内.
6.1.1 模型I
杆I与杆II的长度分别为0.3 m和0.5 m, 两杆的截面均为矩形, 宽与高为0.05 m, 材料参数设为: $\rho = 2000$ kg/m$^{3}$, $E = 8.2$ $\times$ 10$^{12}$ Pa, $\nu = 0.3$. 杆II末端滑块的质量为0.1 kg. 系统初始静止放置于$XOY$水平面内, 杆I与$X$轴夹角为0$^\circ$. 在曲柄上施加$Z$方向力矩900 N$\cdot$m, 持续时间为0.2 s, 杆I的最大转速约为1033 rad/s. 时间积分算法采用广义$\alpha $方法, 仿真时间为0.2 s, 时间步长设为1.0 $\times$ 10$^{-4}$ s, 算法参数谱半径设为0.8.
通过试算, 三阶SE(3)几何精确梁84个自由度收敛, 采用缩减积分技术和Hu-Washizu三场变分原理的一阶SE(3)几何精确梁与一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁均在156个自由度收敛, ANCF全参数梁312个自由度收敛.
图4 给出了上述三类建模方法收敛时杆II变形率随时间的变化曲线. 随着转速增大, 杆II的变形逐渐增大, 最大变形能够达到杆长度的0.04%.
图4 新窗口打开 |
下载原图ZIP |
生成PPT 图4杆II变形率随时间变化曲线 Fig.4Deformation rate of the rod II 6.1.2 模型II
杆I与杆II的长度分别为0.3 m和0.5 m, 两杆的截面均为矩形, 宽与高为0.03 m, 材料参数设为: $\rho = 5600$ kg/m$^{3}$, $E = 4.1$ $\times$ 10$^{11}$ Pa, $\nu = 0.3$. 杆II末端滑块质量为0.1 kg. 系统初始静止放置于$XOY$水平面内, 杆I与$X$轴夹角为0$^\circ$. 在曲柄上施加$Z$方向力矩500 N$\cdot$m, 持续时间为0.2 s, 杆I的最大转速约为518 rad/s. 时间积分算法采用广义$\alpha$方法, 仿真时间为0.2 s, 时间步长设为3.0 $\times$ 10$^{-5}$ s, 算法参数谱半径设为0.8.
通过试算, 三阶SE(3)几何精确梁84个自由度收敛, 采用缩减积分技术和Hu-Washizu三场变分原理的一阶SE(3)几何精确梁与一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁均在300个自由度收敛, ANCF全参数梁504个自由度收敛.
图5 给出了上述三类建模方法收敛时杆II变形率随时间的变化曲线. 随着转速增大, 杆II的变形逐渐增大, 最大变形能够达到杆长度的4%.
图5 新窗口打开 |
下载原图ZIP |
生成PPT 图5杆II变形率随时间变化曲线 Fig.5Deformation rate of the rod II 由
图4 和
图5 可知, SE(3)几何精确梁、$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁和ANCF全参数梁均能达到收敛的数值结果, 验证了这三类建模方法在描述高转速、小变形以及中等变形多体系统动力学问题时的正确性. 采用缩减积分技术与Hu-Washizu三场变分原理可以缓解动力学问题中的剪切闭锁, 两种方法的收敛速度无明显差别. 一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁与一阶SE(3)几何精确梁收敛速度无明显差别.
在一个牛顿迭代步中, 最耗时的两个部分包括计算系统刚度矩阵与迭代矩阵的LU分解/回代. 由于本文中的梁单元均采用局部标架建模方法, 刚度矩阵更新次数大幅降低. 因此, 可通过计算迭代矩阵中非零元素的个数来评价一个牛顿迭代步的计算速度. 对于模型I, 几类一阶几何精确梁迭代矩阵非零元素的个数均为2808, 三阶SE(3)几何精确梁和ANCF全参数梁迭代矩阵非零元素的个数分别为3528和11 232. 对于模型II, 几类一阶几何精确梁迭代矩阵非零元素的个数均为5400, 三阶SE(3)几何精确梁和ANCF全参数梁迭代矩阵非零元素的个数分别为3528和18 144. 由以上数据可知, 在达到相同的计算精度时, 几类一阶几何精确梁计算效率相同, SE(3)和$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁的计算效率明显高于ANCF全参数梁.
其次, 通过分析系统质量矩阵和刚度矩阵更新情况, 验证基于局部标架的SE(3)几何精确方法能够有效地减轻刚体运动带来的非线性问题. 若忽略外力对应的Jacobian矩阵, 系统迭代矩阵的表达式为
(52) $\begin{eqnarray} \label{eq52} { J}_{\rm iter} =\left[\begin{array}{c@{\quad }c} h^{2}\beta ({ M}+{ K}) & \varPhi_{q}^{\rm T}\\ \varPhi_{q} & {\bf 0} \\ \end{array} \right] \end{eqnarray}$ 式中, $h$为积分步长, $\beta $为广义$\alpha $方法算法参数, $ M$为广义质量矩阵, $ K$为切线刚度矩阵, $ \varPhi_{q}$为约束方程的Jacobian矩阵. 本算例中, 牛顿迭代收敛标准设为广义坐标修正量$||\varDelta||^{2}< 1.0$ $\times$ 10$^{-6}$.
表7 给出了一阶SE(3)几何精确梁迭代矩阵的主要部分$ J_{q} = h^{2}\beta ( M+ K)$的更新次数(number of update, NOU)以及仿真过程中牛顿迭代总次数(number of iterations, NOI).
表7 中"NNI $\geqslant i$"表示一个时间步内, 牛顿迭代次数大于等于$i$步时, $ J_{q}$强制更新.
Table 7 表7 表7 不同模型$ J_{q}$的更新次数及仿真过程牛顿迭代总次数
Table 7
Number of updating $ J_{q}$ and number of Newton iterations about different models 新窗口打开 |
下载CSV 模型I的数值结果表明: 对于基于局部标架的SE(3)几何精确梁模型, 若每个牛顿迭代步中$ J_{q}$均更新, 则整个仿真过程共需要进行4713次牛顿迭代步, 同时$ J_{q}$也需要更新4712次. 当NNI $\geqslant 3$时(即每个时间步中牛顿迭代达到三步时更新一次$ J_{q})$, $ J_{q}$需要更新357次, 总迭代次数增加了30%左右. 当NNI $\geqslant 4$时(即每个时间步中牛顿迭代达到四步时更新一次$ J_{q})$, $ J_{q}$需要更新10次, 总迭代次数增加了37%左右. 模型II与模型I情况类似, 这里不再赘述. 由此可得, 对于此类高速转动的多体动力学问题, 基于局部标架的SE(3)几何精确梁能够极大地降低刚体运动带来的几何非线性.
6.2 发生中等变形和超大变形的空间双摆 本节对空间双摆模型进行动力学仿真. 如
图6 所示, 空间双摆模型初始放置于$XOY$水平面上, 并在重力作用下开始运动. 与算例6.1一致, 本算例通过杆II中点$C$到其首尾连线的距离$d$度量杆II的弹性变形, 变形率为$d$与杆原长之间的比值. 同时, 以1944个自由度的$\mathbb{R}^{3}$空间ANCF全参数梁单元所得到的数值结果作为参考值, 建模方法数值结果收敛标准设为变形场最大值与参考解的相对误差在1%以内.
图6 新窗口打开 |
下载原图ZIP |
生成PPT 图6空间双摆 Fig.6A spatial double pendulum 6.2.1 模型III
双摆两个摆臂的长度分别为0.3 m和0.5 m, 两个摆臂的截面均为矩形, 宽与高为0.5 mm, 材料参数设为: $\rho = 2700$ kg/m$^{3}$, $E = 2.1$ $\times$ 10$^{10}$ Pa, $\nu = 0.3$, 重力加速度为$g = 9.81$ m/s$^{2}$. 时间积分算法采用广义$\alpha $方法, 仿真时间1 s, 时间步长设为5.0 $\times$ 10$^{-4}$ s, 算法参数谱半径设为0.8.
通过试算, 三阶SE(3)几何精确梁84个自由度收敛, 采用缩减积分技术和Hu-Washizu三场变分原理的一阶SE(3)几何精确梁均在300个自由度收敛, ANCF全参数梁1464个自由度收敛. 300个自由度的一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁在$t = 0.457 5$ s发散.
图7 给出了上述三类建模方法收敛时杆II变形率随时间的变化曲线, 最大变形能够达到杆长度的5.5%.
图7 新窗口打开 |
下载原图ZIP |
生成PPT 图7杆II变形率随时间变化曲线 Fig.7Deformation rate of the rod II 6.2.2 模型IV
双摆两个摆臂的长度分别为0.3 m和0.5 m, 两个摆臂的截面均为矩形, 宽与高为0.3 mm, 材料参数设为: $\rho = 2700$ kg/m$^{3}$, $E = 1.5$ $\times$ 10$^{10}$ Pa, $\nu = 0.3$, 重力加速度为$g = 9.81$ m/s$^{2}$. 时间积分算法采用广义$\alpha $方法, 仿真时间0.5 s, 时间步长设为5.0 $\times$ 10$^{-4}$s, 算法参数谱半径设为0.8.
通过试算, 三阶SE(3)几何精确梁84个自由度收敛, 采用缩减积分技术和Hu-Washizu三场变分原理的一阶SE(3)几何精确梁156个自由度收敛, ANCF全参数梁1752个自由度收敛. 156个自由度的一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁在$t = 0.444 5$ s发散.
图8 给出了上述三类建模方法收敛时杆II变形率随时间的变化曲线, 最大变形能够达到杆长度的15%.
图8 新窗口打开 |
下载原图ZIP |
生成PPT 图8杆II变形率随时间变化曲线 Fig.8Deformation rate of the rod II 由
图7 和
图8 可知, SE(3)几何精确梁、$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁和ANCF全参数梁均能达到收敛的数值结果, 验证了这三类建模方法在描述中等变形和超大变形多体系统动力学问题时的正确性. 采用缩减积分技术与Hu-Washizu三场变分原理均可有效缓解动力学中的剪切闭锁问题, 两种方法的收敛速度无明显差别.
下面通过计算迭代矩阵中非零元素的个数来评价一个牛顿迭代步的计算速度. 对于模型III, 几类一阶几何精确梁迭代矩阵非零元素的个数均为5400, 三阶SE(3)几何精确梁和ANCF全参数梁迭代矩阵非零元素的个数分别为3528和52 704. 对于模型IV, 几类一阶几何精确梁迭代矩阵非零元素的个数均为2808, 三阶SE(3)几何精确梁和ANCF全参数梁迭代矩阵非零元素的个数分别为3528和63 072. 由以上数据可知, 在达到相同的计算精度时, 几类一阶几何精确梁计算效率相同, SE(3)几何精确梁的计算效率明显高于ANCF全参数梁.
通过与6.1节中等变形的平面多体系统对比, 可以发现, 随着柔性体变形的增大, SE(3)几何精确梁比ANCF全参数梁单元优势更明显. 一阶$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁在仿真计算的过程中发散, 是由于本文仅采用Rescaling技术
[32 ] 处理转动参数, 没有进一步采用其他方法避免转动参数的奇异性. 相比较而言, 本文提出的局部标架建模方法则可自然地避免转动参数的奇异性问题.
与6.1节统计方法类似, 通过分析系统质量矩阵和刚度矩阵更新情况, 以说明基于局部标架的SE(3)几何精确建模方法描述大转动、大变形非线性运动的能力. 本算例中, 牛顿迭代收敛标准设为广义坐标修正量$||\varDelta||^{2}<1.0$ $\times$ 10$^{-6}$.
表8 中给出了一阶SE(3)几何精确梁$ J_{q}$的更新次数以及仿真过程中的牛顿迭代总次数.
Table 8 表8 表8 不同模型$ J_{q}$的更新次数及仿真过程牛顿迭代总次数
Table 8
Number of updating $ J_{q}$ and number of Newton iterations about different models 新窗口打开 |
下载CSV 模型IV的数值结果表明: 对于基于局部标架的SE(3)几何精确梁模型, 若每个牛顿迭代步中$ J_{q}$均更新, 则整个仿真过程共需要进行2479次牛顿迭代步, 同时$ J_{q}$也需要更新2478次. 当NNI $\geqslant 3$时(即每个时间步中牛顿迭代达到三步时更新一次$ J_{q})$, $ J_{q}$需更新440次, 总迭代次数增加了40%左右. 当NNI $\geqslant4$时(即每个时间步中牛顿迭代达到四步时更新一次$ J_{q})$, $ J_{q}$需要更新283次, 总迭代次数增加了54%左右. 模型III与模型IV情况类似, 这里不再赘述. 由此可得, 对于此类大变形多体动力学问题, 基于局部标架的SE(3)几何精确梁仍然能够明显地降低刚体运动带来的几何非线性.
7 结论 本文基于李群SE(3)局部标架, 研究几类梁单元的闭锁缓解方法. 采用Hu-Washizu三场变分原理缓解了SE(3)几何精确梁单元中的剪切闭锁, 采用应变分解法缓解了ANCF全参数梁单元中的泊松闭锁, 并对SE(3)几何精确梁、$\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁和ANCF全参数梁的单元性能进行了算例对比. 静力学与动力学算例结果均表明, 缓解闭锁后的SE(3)几何精确梁计算精度高. 在进行大转动、大变形动力学计算时, $\mathbb{R}^{3}$ $\times$ $\mathbb{R}^{3}$几何精确梁单元会出现转动参数奇异问题, 若处理不当则会不收敛, 而SE(3)几何精确梁单元可避免转动参数奇异性问题. SE(3)几何精确梁的广义质量矩阵主项为常数, 与转动参数相关的切线刚度矩阵满足刚体转动的不变性, 在计算中无须更新, 可极大地减少迭代矩阵的更新次数, 有效地提高计算效率. 基于SE(3)群局部标架的几何精确建模方法更适合用于描述梁结构的大转动、大变形运动.
[1] Meier C Popp A Wall WA Archives of Computational Methods in Engineering 2019 ,26(1 ):163 -243 [本文引用: 1] [2] 孙加亮 , 田强 , 胡海岩 . 多柔体系统动力学建模与优化研究进展力学学报 , 2019 ,51(6 ):1565 -1586 [本文引用: 1] ( Sun Jialiang Tian Qiang Hu Haiyan Chinese Journal of Theoretical and Applied Mechanics 2019 ,51(6 ):1565 -1586 (in Chinese)) [本文引用: 1] [3] Schulz M B?l M Computational Mechanics 2019 : 1 -21 [本文引用: 1] [4] Choi MJ Cho S Computer Methods in Applied Mechanics and Engineering 2019 ,351:153 -183 [本文引用: 1] [5] Duan LP Zhao JC Computers and Structures 2019 ,218:32 -59 [本文引用: 1] [6] Lestringant C Audoly B Kochmann DM Computer Methods in Applied Mechanics and Engineering 2020 ,361:112741 [本文引用: 1] [7] Shabana AA Chicago: University of Illinois at Chicago , 1996 [本文引用: 1] [8] 范纪华 , 章定国 , 谌宏 . 基于绝对节点坐标法的弹性线方法研究力学学报 , 2019 ,51(5 ):1455 -1465 [本文引用: 1] ( Fan Jihua Zhang Dingguo Shen Hong Chinese Journal of Theoretical and Applied Mechanics 2019 ,51(5 ):1455 -1465 (in Chinese)) [本文引用: 1] [9] 吴吉 , 章定国 , 黎亮 等 . 带集中质量的旋转柔性曲梁动力学特性分析力学学报 , 2019 ,51(4 ):1134 -1147 [本文引用: 1] ( Wu Ji Zhang Dingguo Li Liang et al . Dynamic characteristics analysis of a rotating flexible curved beam with a concentrated massChinese Journal of Theoretical and Applied Mechanics 2019 ,51(4 ):1134 -1147 (in Chinese)) [本文引用: 1] [10] 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 ANCFComputer Methods in Applied Mechanics and Engineering 2013 ,258:81 -95 [本文引用: 1] [11] Rong JL Wu ZP Liu C et al . Geometrically exact thin-walled beam including warping formulated on the special Euclidean group SE(3)Computer Methods in Applied Mechanics and Engineering 2020 ,369:113062 [本文引用: 3] [12] Zienkiewicz OC Taylor R Too JM International Journal for Numerical Methods in Engineering 1971 ,3(2 ):275 -290 [本文引用: 1] [13] Zienkiewicz OC Owen DRJ Lee KN International Journal for Numerical Methods in Engineering 1974 ,8(2 ):341 -358 [本文引用: 1] [14] Hughes TJR Taylor RL Kanoknukulchai W International Journal for Numerical Methods in Engineering 1977 ,11(10 ):1529 -1543 [本文引用: 1] [15] Simo JC Vu-Quoc L Computer Methods in Applied Mechanics and Engineering 1986 ,58(1 ):79 -116 [本文引用: 1] [16] Malkus DS Hughes TJR Computer Methods in Applied Mechanics and Engineering 1978 ,15(1 ):63 -81 [本文引用: 1] [17] Noor AK Peters JM International Journal for Numerical Methods in Engineering 1981 ,17(4 ):615 -631 [本文引用: 1] [18] Pian THH Finite Elements in Analysis and Design 1985 ,1(2 ):135 -140 [本文引用: 1] [19] Liu WK Belytschko T Chen J Computer Methods in Applied Mechanics and Engineering 1988 ,71(3 ):241 -258 [本文引用: 1] [20] Dorfi HR Busby HR Computers and Structures 1994 ,53(1 ):43 -52 [本文引用: 1] [21] Simo JC Rifai MS International Journal for Numerical Methods in Engineering 1990 ,29(8 ):1595 -1638 [本文引用: 1] [22] Simo JC Armero F International Journal for Numerical Methods in Engineering 1992 ,33(7 ):1413 -1449 [本文引用: 1] [23] Andelfinger U Ramm E International Journal for Numerical Methods in Engineering 1993 ,36(8 ):1311 -1337 [本文引用: 1] [24] Pian THH Sumihara K International Journal for Numerical Methods in Engineering 1984 ,20(9 ):1685 -1695 [本文引用: 1] [25] Gerstmayr J Shabana AA Nonlinear Dynamics 2006 ,45(1-2 ):109 -130 [本文引用: 1] [26] Gerstmayr J Matikainen MK Mechanics Based Design of Structures and Machines 2006 ,34(4 ):409 -430 [本文引用: 1] [27] Kerkkaenen KS Sopanen JT Mikkola AM Journal of Mechanical Design 2005 ,127(4 ):621 -630 [本文引用: 1] [28] Gerstmayr J Matikainen MK Mikkola A Multibody System Dynamics 2008 ,20(4 ):359 -384 [本文引用: 1] [29] Matikainen MK Dmitrochenko O Mikkola A Greece , 2010 : 1266 -1270 [本文引用: 1] [30] Patel M Shabana AA Acta Mechanica 2018 ,229(7 ):2923 -2946 [本文引用: 2] [31] 刘铖 , 胡海岩 . 基于李群局部标架的多柔体系统动力学建模与计算力学学报 , 2021 ,53(1 ), doi: 10.6052/0459-1879-20-292 [本文引用: 1] ( Liu Cheng Hu Haiyan Chinese Journal of Theoretical and Applied Mechanics 2021 ,53(1 ), doi: 10.6052/0459-1879-20-292 (in Chinese)) [本文引用: 1] [32] Olivier AB Springer Netherlands , 2011 : 534 -537 [本文引用: 1] 
,*,2), 刘绍奎?
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
新窗口打开|下载原图ZIP|生成PPT
