范纪华^*,?, 章定国^,**,2), 谌宏^?* 江苏科技大学机电与动力工程学院, 江苏张家港 215600
? 江苏科技大学苏州理工学院, 江苏张家港 215600
** 南京理工大学理学院, 南京 210094

RESEARCH ON ELASTIC LINE METHOD BASED ON ABSOLUTE NODAL COORDINATE METHOD¹⁾

Fan Jihua^*,?, Zhang Dingguo^,**,2), Shen Hong^?* School of Mechatronics and Power Engineering, Jiangsu University of Science and Technology, Zhangjiagang 215600, Jiangsu China
? Suzhou Institute of Technology, Jiangsu University of Science and Technology, Zhangjiagang 215600, Jiangsu China
** School of Sciences, Nanjing University of Science and Technology, Nanjing 210094, China

通讯作者: 2) 章定国, 教授, 主要研究方向：多体系统动力学与控制.E-mail:zhangdg419@mail.njust.edu.cn

收稿日期:2019-03-28网络出版日期:2019-09-18

基金资助:

1) 国家自然科学基金项目.11502098 国家自然科学基金项目.11772158 江苏省高校自然科学研究面上项目.15KJB130003 江苏科技大学博士科研启动基金项目资助.120140003

Received:2019-03-28Online:2019-09-18
作者简介 About authors


摘要
相比于浮动坐标系法, 绝对节点坐标法(absolute nodal coordinateformulation, ANCF)在处理柔性体非线性大变形问题上具有显著优势,ANCF将单元节点坐标定义在全局坐标系下,采用斜率矢量代替节点转角坐标, 具有常数质量阵,不存在科氏离心力等优点, 然而弹性力阵为非线性项,其求解将比较耗时且占用资源. 据此, 在弹性力求解方法中,引入弹性线方法(elastic line method, ELM),该方法将格林--拉格朗日应变张量定义在中心线上,采用曲率公式来定义弯曲应变, 转角公式来定义扭转应变.同时采用有限元法对三维柔性梁位移场进行离散,求解梁单元常数质量阵、广义刚度阵、广义力阵,进而得到单元的动力学方程, 通过转换矩阵得到三维梁的动力学方程.接着从理论上指出连续介质力学方法(continuum mechanics method,CMM)和弹性线方法在求解弹性力上的不同点, 并编制动力学仿真软件.最后分别采用连续介质力学方法和弹性线方法对柔性单摆以及履带式车辆的动力学问题进行仿真分析,结果表明：弹性线方法能在保证精度的前提下有效提高计算效率.
关键词： 绝对节点坐标法;连续介质力学方法;弹性线方法;动力学仿真;履带式车辆

Abstract
Compared with the floating frame of reference formulation, the absolute nodal coordinate formulation (ANCF) has significant advantages in dealing with the nonlinear large deformation problem of flexible bodies. ANCF defines the nodal co-ordinates in a global co-ordinate system, uses the global slopes instead of angles to define the orientation of the elements, has a constant mass matrix, and does not have the Coriolis centrifugal force. However,the elastic force matrix is a nonlinear term, and the solution will be time consuming and occupied resources. According to this,the elastic line method is introduced for solving the elastic force, the Green-Lagrangian strain tensor is defined on the centerline, the curvature formula is used to define the bending strain, and the torsion angle formula is used to define the torsional strain. At the same time, the finite element method is used to describe the displacement field of the three-dimensional flexible beam, and the constant mass matrix, the generalized stiffness matrix and the generalized force matrix of the beam element are solved, and then the dynamic equations of the element are obtained. The dynamic equations of the three-dimensional beam are obtained by the transformation matrix. Then the differences between the continuum mechanics method and the elastic line method are pointed out theoretically, and the dynamic simulation software is compiled. Finally, the dynamics behaviors of a flexible pendulum and a tracked vehicle are numerical investigated with the continuum mechanics method and the elastic line method. The results show that the elastic line method can effectively improve the calculation efficiency under the premise of ensuring accuracy.
Keywords：absolute nodal coordinate formulation;continuum mechanics method;elastic line method;dynamic simulation;tracked vehicle


PDF (1864KB)元数据多维度评价相关文章导出EndNote|Ris|Bibtex收藏本文
本文引用格式
范纪华, 章定国, 谌宏. 基于绝对节点坐标法的弹性线方法研究¹⁾. 力学学报[J], 2019, 51(5): 1455-1465 DOI:10.6052/0459-1879-19-076
Fan Jihua, Zhang Dingguo, Shen Hong. RESEARCH ON ELASTIC LINE METHOD BASED ON ABSOLUTE NODAL COORDINATE METHOD¹⁾. Chinese Journal of Theoretical and Applied Mechanics[J], 2019, 51(5): 1455-1465 DOI:10.6052/0459-1879-19-076

引言

随着航天技术的飞速发展,以大跨度、低刚度、轻质、高速的柔性机械臂、卫星天线、太阳能帆板等为代表的柔性构件越来越多地出现在空间飞行器中.这类柔性构件在工作过程中往往会由构件本身的大位移运动引起自身的弹性变形,大位移运动和自身的弹性变形存在着非线性强耦合效应,而传统的多刚体系统模型已无法正确描述该类系统的动力学性态^[1-2]. 在此背景下各国****从不同的参考坐标系出发,建立了目前运用较广泛的有浮动坐标系 法^[3-15]和绝对节点坐标法^[16-38].浮动坐标系法起源于多刚体系统动力学,该方法将柔性体的运动分解为随浮动坐标系的大位移运动和相对于浮动坐标系的弹性变形运动,该建模方法有利于变形描述的线性化,通常对于柔性体的小变形问题具有较高的计算精度和效率,但其在动力学方程建立中必须经历复杂的坐标转换,同时在处理柔性体大变形或复杂变形问题时, 会显得力不从心;绝对节点坐标法基于连续介质力学理论和非线性有限元法^[16-17],其将柔性体的大位移及其弹性变形都用相对于绝对坐标系的单元节点坐标表示,直接建立柔性体的应变--位移关系^[18],该方法具有能有效避免坐标转换的繁琐过程,并具有程式化好、编程方便等优点, 不过其也具有计算效率较低的缺点,因此合理优化刚度矩阵是绝对节点坐标法在多体系统动力学研究中的重点以及难点之一.

Berzeri等^[19]采用绝对节点坐标法对二维欧拉伯努利梁进行研究,并建立了简单精确的四种弹性力求解模型,其模型考虑了弯曲变形和轴向变形,以及弯曲变形和轴向变形的相互耦合作用. Omar等^[20]采用有限元法和连续介质力学理论对二维铁木辛柯梁进行研究,并通过非线性应变位移关系来求解弹性力,该文献中对于形函数的定义和刚度矩阵的求解存在一些错误.Dufva等^[21]采用绝对节点坐标法对二维铁木辛柯梁进行研究,并将仿真结果与Ansys软件仿真结果进行对比.随后Shabana和 Yakoub^[22-23]对三维梁模型的建模理论、仿真应用进行了研究,通过非线性的格林--拉格朗日应变张量和连续介质力学理论对三维梁弹性力进行求解.田强等^[24-25]利用绝对节点坐标法建立了含间隙铰关节、曲梁和圆柱壳单元的动力学模型,对绝对节点坐标法的单元模型进行了补充和拓展.随后文献^[26,27]将绝对节点坐标法应用到含芯拧绞绳非线性弯曲动力学特性分析和车辆钢板弹簧模型动力学建模中.郑彤等^[28]对三维大变形柔性梁系统的动力学建模与仿真进行研究,将仿真结果和ADAMS软件计算的结果进行对比,指出ADAMS在计算大变形物体动力学时具有局限性.赵春璋等^[29]对复杂变截面柔性梁结构进行研究,提出在不影响梁刚度的情况可通过改变截面参数有效降低梁的质量.文献^[30]将传热与连续介质力学相结合,提出了对位移和温度场的统一描述,构建了一个热集成的绝对节点坐标公式全参数化梁单元,同时航天器中刚柔热耦合效应被重点关注^[31-32]. 为了优化刚度矩阵, 提升绝对节点坐标法计算效率,Garcia-Vallejo等^[33]对三维梁模型弹性力的求解进行简化,提出不变矩阵法.章孝顺等^[34]运用拉格朗日恒等式给出了柔性梁横向弯曲变形能新的表达式,通过新的变形能表达式得到新的弹性力模型,进而提高计算效率.陈渊钊等^[35]为了消除或减弱传统ANCF中缩减单元的"失真现象",提出了柔性梁基于无网格径向基点插值的ANCF法^[36], 其将有效提高计算精度和计算效率.文献^[37]在ANCF低阶单元上开发了一种数学分析方法,进而提高计算效率.文献^[38]提出了一种基于正交分解和Galerkin投影的模型阶数约简的系统方法来提高计算效率,并通过刚性和柔性多体系统仿真实例验证方法的有效性.由上述研究可知,计算效率仍然是绝对节点坐标法目前急需解决的难点之一.

本文采用有限元法对三维梁单元位移场进行离散,通过求解梁单元常数质量阵、广义刚度阵、广义力阵,进而建立单元的动力学方程, 通过转换矩阵得到三维梁的动力学方程.接着重点给出求解弹性力的两种方法：连续介质力学方法与弹性线方法,并从理论上对这两种方法进行分析比较. 最后,基于上述理论对绝对节点坐标法下柔性单摆的动力学问题进行研究,并将连续介质力学方法和弹性线方法用于履带式车辆动力学模型中,进而通过仿真实例说明弹性线方法在计算效率上的突出提升.

1 三维梁单元位移场离散

如图1所示三维梁单元, 其上任意一点$P$的全局位移矢量$r$可表示为

(1)$$r = {\left[ {\begin{array}{l} r_1 \\ r_2 \\ r_3 \\ \end{array}} \right] = \left[ {\begin{array}{l} X \\ Y \\ Z \\ \end{array}} \right] }= \\ \left[\!\!\! {\begin{array}{l} a_0 \!+\! a_1 x \!+\! a_2 y \!+\! a_3 z \!+\! a_4 xy \!+\! a_5 xz \!+\! a_6 x^2 \!+\! a_7 x^3 \\ b_0 \!+\! b_1 x \!+\! b_2 y \!+\! b_3 z \!+\! b_4 xy \!+\! b_5 xz \!+\! b_6 x^2 \!+\! b_7 x^3 \\ c_0 \!+\! c_1 x \!+\! c_2 y \!+\! c_3 z \!+\! c_4 xy \!+\! c_5 xz \!+\! c_6 x^2 \!+\! c_7 x^3 \\ \end{array}} \!\!\!\right]\!= \! \\ {Se} $$

图1

新窗口打开|下载原图ZIP|生成PPT
图1三维梁单元模型

Fig. 1The model of 3D beam element



其中, $S$是单元形函数, $x$为沿着轴线的单元局部坐标,取值范围为$\left[ {0, l} \right]$. $l$为单元的初始长度, $e$为单元节点坐标的广义坐标矢量

(2)\begin{equation}\label{eq2} e = [ e^{1\rm T} \quad e^{2\rm T}]^{\rm T}\end{equation}

式中, 上标$k = 1,2$代表单元的节点编号,单元节点广义坐标矢量在节点$k$包括节点的三个全局位移以及九个全局斜率,具体形式如下

(3)\begin{equation} \label{eq3} e^k = \left[ {r_1^k r_2^k {\kern 1pt} r_3^k {\kern 1pt} \frac{\partial r_1^k }{\partial x}{\kern 1pt} {\kern 1pt} \frac{\partial r_2^k }{\partial x}{\kern 1pt} {\kern 1pt} \frac{\partial r_3^k }{\partial x}{\kern 1pt} {\kern 1pt} \frac{\partial r_1^k }{\partial y}{\kern 1pt} {\kern 1pt} \frac{\partial r_2^k }{\partial y}{\kern 1pt} {\kern 1pt} \frac{\partial r_3^k }{\partial y}{\kern 1pt} {\kern 1pt} \frac{\partial r_1^k }{\partial z}{\kern 1pt} {\kern 1pt} \frac{\partial r_2^k }{\partial z}{\kern 1pt} {\kern 1pt} \frac{\partial r_3^k }{\partial z}} \right]^{\rm T} \end{equation}

通过单元函数表达式和单元始端、末端的函数值进行推导,可以得到单元形函数$S$的表达式为

(4)$$S = \left[ {s_1 I_3~~s_2 I_3 ~~ s_3 I_3 ~~ s_4 I_3 ~~s_5 I_3 ~~s_6 I_3 ~~ s_7 I_3 ~~s_8 I_3 }\right] $$

其中, $I_3 $为3$\times $3单位矩阵, 形函数$s_i = s_i (\xi ,\eta,\zeta ),~ i = 1,2, \cdots, 8$形式为

(5)\begin{equation} \label{eq4} \left. \begin{array}{l} s_1 = 1-3\xi ^2 + 2\xi ^3, \quad s_2 = l(\xi-2\xi ^2 + \xi ^3) \\ s_3 = l(\eta-\xi \eta ), \quad s_4 = l(\zeta-\zeta \xi ) \\ s_5 = 3\xi ^2-2\xi ^3, \quad s_6 = l(-2\xi ^2 + \xi ^3) \\ s_7 = l\xi \eta , \quad s_8 = l\zeta \xi \\ \end{array} \right\} \end{equation}

其中, $\xi = x / l$, $\eta = y / l$, $\zeta = z / l$,$x$的取值范围为$\left[ {0,l} \right]$, $y$的取值范围为$\left[ {-b / 2, b / 2} \right]$, $z$的取值范围为$\left[ {-h / 2, h / 2}\right]$, $l$代表单元的长度, $b$代表单元的宽度,$h$代表单元的厚度.

2 绝对节点坐标下的动力学方程

采用有限元法对柔性梁位移场进行描述,柔性梁单元上任意点的位移可由广义坐标矢量$e$与单元形函数$S$表示,这样就使连续无限自由度的动力学模型转化为有限自由度的动力学模型,进而模型的动力学方程就可数值求解.绝对节点坐标下三维梁单元的动力学方程形式为

(6)\begin{equation}\label{eq5} M\ddot { e} + K( e) e = { Q}_e\end{equation}

式中, $M$为单元常数质量阵; 弹性力${ Q}_k = K( e) e$, 其中$K( e)$为广义刚度阵, ${ Q}_e $为广义力阵.

2.1 常数质量阵

梁单元上任意点的全局位移矢量已经通过方程(1)得到,将方程对时间求导数, 可得到全局位移速度矢量

(7)\begin{equation}\label{eq6} \dot { r} = S\dot { e} \end{equation}

通过单元的全局节点速度矢量, 可得到单元的动能为

(8)$$T = \frac{1}{2}\int_V {\rho \dot { r}^{\rm T} \dot{ r}\mbox{d}V} = \\ \frac{1}{2}\dot { e}^{\rm T} \left(\int_V {\rho S^{\rm T} S{\rm d}V} \right)\dot { e} = \frac{1}{2}\dot { e}^{\rm T} M\dot { e} $$

其中, $V$为梁单元的体积, $\rho $为梁的材料密度, 则$M$即为单元的广义质量阵, 上述方程的广义质量阵为常对称矩阵,通过方程(4)给出的形函数, 可以得到广义质量阵.

2.2 广义刚度阵

由非线性Green-Lagrange应力--位移关系计算任意质点的变形梯度为

(9)$$J = \frac{\partial r}{\partial X} = \frac{\partial r}{\partial \overline { X}} \left[\frac{\partial X}{\partial \overline { X}}\right]^{ -1} = \frac{\partial r}{\partial \overline { X}} J_0^{-1} $$

(10)$$J_0 = \frac{\partial X}{\partial \overline { X}} $$

其中, $\overline { X}$为单元坐标系, $\overline { X} =[x,y,z]$定义单元坐标系下的某个点的位置矢量.当梁的初始位置是任意的$J_0 $必须考虑, 而当梁初始位置是水平时,$J_0 $为单位矩阵.

通过变形梯度方程可得到应变张量为

(11)$$\varepsilon = \frac{1}{2}( J^{\rm T} J- I) = \left[ {\begin{array}{ccc} \varepsilon _{11} \ \varepsilon _{12} \ \varepsilon _{13} \\ \varepsilon _{12} \ \varepsilon _{22} \ \varepsilon _{23} \\ \varepsilon _{13} \ \varepsilon _{32} \ \varepsilon _{33} \\ \end{array}} \right] = \\ \left[ {\begin{array}{ccc} \varepsilon _{11} \ \dfrac{1}{2}\gamma _{12} \ \dfrac{1}{2}\gamma _{13} \\[3mm] \dfrac{1}{2}\gamma _{12} \ \varepsilon _{22}\ \dfrac{1}{2}\gamma _{23} \\[3mm] \dfrac{1}{2}\gamma _{13} \ \dfrac{1}{2}\gamma _{23} \ \varepsilon _{33} \\ \end{array}} \right] = \\ \dfrac{1}{2}\left[ {\begin{array}{ccc} r_x^{\rm T} r_x-1 \ r_x^{\rm T} r_y \ r_x^{\rm T} r_z \\ r_x^{\rm T} r_y \ r_y^{\rm T} r_y- 1 \ r_y^{\rm T} r_z \\ r_x^{\rm T} r_z \ r_y^{\rm T} r_z \ r_z^{\rm T} r_z-1 \\ \end{array}} \right] $$

其中, $r_\alpha = \dfrac{\partial r}{\partial \alpha },~\alpha = x,y,z$. 结合位移公式$r = S(x,y,z) e$, 则

(12)\begin{equation} \label{eq11} r_x = S_x e = S_1 e,\quad r_y = S_y e = S_2 e,\quad r_z = S_z e = S_3 e \end{equation}

则应变张量方程(11)的各个分量表达式

(13)\begin{equation} \label{eq12} \left. \begin{array}{l} \varepsilon _{11} = \dfrac{1}{2}( r_x^{\rm T} r_x-1) = \dfrac{1}{2}( e^{\rm T} S_x^{\rm T} S_x e -\mbox{1}) \\[3mm] \varepsilon _{22} = \dfrac{1}{2}( r_y^{\rm T} r_y-1) = \dfrac{1}{2}( e^{\rm T} S_y^{\rm T} S_y e -1) \\[3mm] \varepsilon _{33} = \dfrac{1}{2}( r_z^{\rm T} r_z-1) = \dfrac{1}{2}( e^{\rm T} S_z^{\rm T} S_z e-1) \\[3mm] \varepsilon _{12} = \dfrac{1}{2} r_x^{\rm T} r_y = \dfrac{1}{2} e^{\rm T} S_x^{\rm T} S_y e \\[3mm] \varepsilon _{13} = \dfrac{1}{2} r_x^{\rm T} r_z = \dfrac{1}{2} e^{\rm T} S_x^{\rm T} S_z e \\[3mm] \varepsilon _{23} = \dfrac{1}{2} r_y^{\rm T} r_z = \dfrac{1}{2} e^{\rm T} S_y^{\rm T} S_z e \\ \end{array} \right\}\end{equation}

通过应变表示式以及应力应变关系式可得到单元应变能为

(14)$$U_e = \frac{1}{2}\int_V { \sigma \varepsilon } \mbox{d}V = \frac{1}{2}\int_V (\sigma _{11} \varepsilon _{11} + \sigma _{12} \varepsilon _{12} +\\ \sigma _{13} \varepsilon _{13} + \sigma _{22} \varepsilon _{22} + \sigma _{23} \varepsilon _{23} + \sigma _{33} \varepsilon _{33} ) \mbox{d}V =\\ \int_V \bigg[\frac{\lambda + 2G}{2} (\varepsilon _{11}^2 + \varepsilon _{22}^2 + \varepsilon _{33}^2 ) + \lambda (\varepsilon _{11} \varepsilon _{22} + \\ \varepsilon _{11} \varepsilon _{33} + \varepsilon _{22} \varepsilon _{33} ) + G(\varepsilon _{12}^2 + \varepsilon _{13}^2 + \varepsilon _{23}^2 )\bigg]\mbox{d}V $$

其中, $V$为单元的体积, $\lambda $和$G$为材料的模量常数,则应变能产生的弹性广义力矢量为

(15)$$Q_k =-\frac{\partial U_e }{\partial e} =-\int_V \bigg[\frac{\lambda + 2G}{2}\bigg(2\varepsilon _{11} \frac{\partial \varepsilon _{11} }{\partial e} + \\ 2\varepsilon _{22} \frac{\partial \varepsilon _{22} }{\partial e} + 2\varepsilon _{33} \frac{\partial \varepsilon _{33} }{\partial e}\bigg) + \\ \lambda \bigg(\varepsilon _{11} \frac{\partial \varepsilon _{22} }{\partial e} + \frac{\partial \varepsilon _{11} }{\partial e}\varepsilon _{22} + \varepsilon _{11} \frac{\partial \varepsilon _{33} }{\partial e} + \frac{\partial \varepsilon _{11} }{\partial e}\varepsilon _{33} + \\ \varepsilon _{22} \frac{\partial \varepsilon _{33} }{\partial e} + \frac{\partial \varepsilon _{22} }{\partial e}\varepsilon _{33} \bigg) + \\ G \bigg(2\varepsilon _{12} \frac{\partial \varepsilon _{12} }{\partial e} + 2\varepsilon _{13} \frac{\partial \varepsilon _{13} }{\partial e} + 2\varepsilon _{23} \frac{\partial \varepsilon _{23} }{\partial e}\bigg)\bigg]\mbox{d}V $$

其中

(16)\begin{equation} \label{eq15} \frac{\partial \varepsilon _{ij} }{\partial e} = \frac{1}{2}( S_i^{\rm T} S_j + S_i^{\rm T} S_j ) e,~~i,j = 1,2,3,~j \ge i \end{equation}

则方程(15)转化为

(17)$$Q_k =-\frac{\partial U_e }{\partial e} =\\-\sum\limits_{\alpha = 1}^3 {\frac{\lambda + 2G}{2}\int_V ( e^{\rm T} S_\alpha ^{\rm T} S_\alpha } e e^{\rm T} S_\alpha ^{\rm T} S_\alpha- e^{\rm T} S_\alpha ^{\rm T} S_\alpha )\mbox{d}V-\\ \sum\limits_{\alpha = 1}^3 {\sum\limits_{ \beta = 1 \atop \beta \ne \alpha }^3 {\frac{\lambda }{2}} } \int_V {( e^{\rm T} S_\alpha ^{\rm T} S_\alpha } e e^{\rm T} S_\beta ^{\rm T} S_\beta - e^{\rm T} S_\alpha ^{\rm T} S_\alpha )\mbox{d}V-\\ \sum\limits_{\alpha = 1}^3 {\sum\limits_{ \beta = 1 \atop \beta \ne \alpha }^3 G } \int_V ( e^{\rm T} S_\alpha ^{\rm T} S_\beta e e^{\rm T} S_\alpha ^{\rm T} S_\beta )\mbox{d}V $$

对上述方程进一步简化得

(18)\begin{equation} \label{eq17} Q_k =-\frac{\partial U_e }{\partial e} =- K( e) e \end{equation}

其中

(19)$$K( e) = \sum\limits_{\alpha = 1}^3 {\frac{\lambda + 2G}{2}\int_V ( S_\alpha ^{\rm T} S_\alpha } e e^{\rm T} S_\alpha ^{\rm T} S_\alpha- S_\alpha ^{\rm T} S_\alpha )\mbox{d}V +\\ \sum\limits_{\alpha = 1}^3 {\sum\limits_{ \beta = 1 \atop \beta \ne \alpha }^3 {\frac{\lambda }{2}} } \int_V ( S_\alpha ^{\rm T} S_\alpha e e^{\rm T} S_\beta ^{\rm T} S_\beta- S_\alpha ^{\rm T} S_\alpha )\mbox{d}V +\\ \sum\limits_{\alpha = 1}^3 {\sum\limits_{ \beta = 1 \atop \beta \ne \alpha }^3 G } \int_V ( S_\alpha ^{\rm T} S_\beta e e^{\rm T} S_\alpha ^{\rm T} S_\beta )\mbox{d}V $$

文献^[33]提出不变矩阵法, 通过该方法将广义刚度矩阵$K( e)$分解为一常数刚度矩阵和一个非线性刚度矩阵,其将大大提高仿真的计算效率.

2.3 广义外力阵

采用虚功原理对广义外力矢量进行求解, 作用在单元上任意一点的外力$F$所做的虚功为

(20)\begin{equation} \label{eq18} \delta up W_e = F^{\rm T}\delta up r = F^{\rm T} S\delta up e = Q_e^{\rm T} \delta up e \end{equation}

其中, $r$为外力作用点的位置矢量, $Q_e^{\rm T} = F^{\rm T} S$代表等效到该单元节点的等效广义力.考虑重力作用下的广义外力, 假设单元的截面面积在变形过程中恒为常数,梁单元材料为均匀材质, 即密度为常量.单元质心的全局位移矢量可参考文献^[19]得到

(21)\begin{equation} \label{eq19} r_G = \left[ {\frac{1}{2} I~~ \frac{l}{12} I~~ 0 I~~ 0 I~~ \frac{1}{2} I~~ - \frac{l}{12} I~~0 I~~ 0 I } \right] e \end{equation}

其中, $I$为3$\times $3单位矩阵, $y$方向为重力加速度作用方向,则重力的虚功为

(22)$$\delta up W_g = \delta up (m[0-\!g~~ 0] r_G ) = m[0 -\!g~~ 0]\delta up r_G =\\ -mg \bigg[0~~\frac{1}{2}~~ 0~~ 0~~ \frac{l}{12}~~ 0~~ 0~~ 0~~ 0~~0~~ 0~~ 0~~ 0 \\ \frac{1}{2}~~ 0~~ 0 ~~ -\frac{l}{12}~~ 0 ~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0 \bigg]\delta up e $$

则可知单元上的重力为

(23)$$Q_g =-mg\bigg[0~~\frac{1}{2} ~~ 0~~ 0 ~~ \frac{l}{12}~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0 \\ 0 ~~ 0~~ \frac{1}{2}~~ 0~~0 ~~ - \frac{l}{12}~~ 0~~ 0~~ 0 ~~0~~ 0~~ 0~~ 0 \bigg]^{\rm T} $$

考虑外力矩作用在梁截面上, 该外力矩绕$Z$轴旋转(如图1),则该外力矩的虚功为

(24)\begin{equation} \label{eq22} \delta up W_{Mz} = M_z \delta up \gamma \end{equation}

其中, $\gamma $为单元梁截面绕$Z$轴的旋转角,则通过坐标旋转矩阵可得到

(25)$$\left[ {\begin{array}{cc} \cos \gamma & -\sin \gamma \\ \sin \gamma & \cos \gamma \\ \end{array}} \right] = \dfrac{1}{f}\left[ {\begin{array}{cc} \dfrac{\partial r_2 }{\partial y} \ \dfrac{\partial r_1 }{\partial y} \\ -\dfrac{\partial r_1 }{\partial y} \ \dfrac{\partial r_2 }{\partial y} \\ \end{array}} \right] \\ f = \sqrt {\left(\frac{\partial r_1 }{\partial y}\right)^2 + \left(\frac{\partial r_2 }{\partial y}\right)^2} $$

则可知$\sin \gamma , \cos \gamma $为

(26)\begin{equation} \label{eq24} \sin \gamma =-\frac{1}{f}\frac{\partial r_1 }{\partial y} , \quad\cos \gamma = \frac{1}{f}\frac{\partial r_2 }{\partial y} \end{equation}

通过方程(25)和(26), 可以得到外力矩产生的虚拟旋转角$\gamma $为

(27)\begin{equation} \label{eq25} \delta up \gamma = \frac{\left(\dfrac{\partial r_2 }{\partial y}\right)\delta up \left(\dfrac{\partial r_1 }{\partial y}\right)-\left(\dfrac{\partial r_1 }{\partial y}\right)\delta up \left(\dfrac{\partial r_2 }{\partial y}\right)}{f^2} \end{equation}

若外力矩作用在图1中单元的节点$A$处, 则外力矩$M_z $产生的广义力矢量为

(28)$$Q_{Mz} = [0~~ 0~~ 0~~ 0~~ 0~~0~~ \frac{M_z e_8 }{f_{Ay}^2 } -\frac{M_z e_7 }{f_{Ay}^2 } 0~~ 0 \\ 0~~ 0~~ 0~~ 0~~ 0~~0~~ 0~~ 0~~ 0~~ 0~~ 0~~0 ~~0~~0 ]^{\rm T} $$

同理, 可以得到梁截面受到绕$Y$轴旋转的外力矩$M_y $产生的作用在图1中单元节点$A$处的广义力 矢量

(29)$$Q_{My} = [0~~0~~0~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0~~ \frac{M_y e_{12} }{f_{Az}^2 } ~~ 0~~ -\frac{M_y e_{10} }{f_{Az}^2 } \\ 0~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0~~ 0~~0~~ 0~~ 0 ]^{\rm T} $$

三维单元模型上受到的集中力、均匀分布力以及外力矩的具体表达方程已经给出,则可根据动力学模型所受广义外力得出广义外力矢量形式.通过求解得到梁单元常数质量阵、广义刚度阵、广义力阵,进而根据梁单元节点坐标列阵对于总体节点坐标列阵的转换矩阵^[28,34],就能得到梁整体的动力学方程.

3 弹性力求解的连续介质力学方法与弹性线方法

基于绝对节点坐标法的质量矩阵是常数, 这也增加了求解弹性力的繁琐性.本节详细介绍求解弹性力的两种方法：一般连续介质力学方法和弹性线方法,并从理论上分析比较两种方法的优缺点.

3.1 连续介质力学方法

非线性格林--拉格朗日应变张量如方程(11)所示, 使用虚功的原理,广义的弹性力的矢量$Q_k $被定义为

(30)\begin{equation} \label{eq28} Q_k =-\int_v {\left(\frac{\partial \varepsilon }{\partial e}\right)} ^{\rm T} E \varepsilon \mbox{d}v \end{equation}

其中, $E$为弹性系数矩阵, $e$为单元节点坐标矢量, $Q_k $的具体求解方式详见2.2节.

3.2 弹性线方法

在弹性线方法中格林--拉格朗日应变张量被定义在中心线上,而采用曲率公式来定义弯曲应变, 转角公式来定义扭转应变.三维梁单元中心线上的位移梯度矢量表达式为

(31)\begin{equation} \label{eq29} r_{0x} = r_x (x,0,0),\quad r_y = r_y (x,0,0), \quad r_z = r_z (x,0,0) \end{equation}

其中, 矢量$r_{0x} $位于中心线的切线方向, 而$r_y $和$r_z $定义在中心线任意点$x$处的横截面上.由于位移场被定义为空间坐标$y$和$z$的线性关系, 因此$r_y $和$r_z $不需要下标"0".三维梁单元的中心线的应变张量可以通过拉格朗日应变张量和曲率表示

(32)$$\varepsilon _{011} = \frac{1}{2}( r_{0x}^{\rm T} r_{0x}-1),\quad \varepsilon _{22} = \frac{1}{2}( r_y^{\rm T} r_y-1)\\ \ \varepsilon _{33} = \frac{1}{2}( r_z^{\rm T} r_z-1) \\ \ \gamma _{012} = r_{0x}^{\rm T} r_y ,\quad \gamma _{013} = r_{0x}^{\rm T} r_z ,\quad \gamma _{23} = r_y^{\rm T} r_z \\ \ \theta = \frac{1}{2}({ r}'^{\rm T}_y { r}_z-{ r}_y^{\rm T} { r}'_z ) ,\quad \bar {\kappa }_{0y} = \frac{ r_z^{\rm T} r_{0xx} }{\left| { r_z } \right|\left| { r_{0xx} } \right|} \\ \ \bar {\kappa }_{0z} =-\frac{ r_y^{\rm T} r_{0xx} }{\left| { r_y } \right|\left| { r_{0xx} } \right|} $$

其中, ${(}')$表示对空间坐标$x$求导, $\varepsilon _{011} ,\varepsilon _{22} ,\varepsilon _{33} $表示正应变, $\gamma _{012},\gamma _{013} ,\gamma _{23} $表示剪切应变, $\theta $表示扭曲转角,$\bar {\kappa }_{0y} ,\bar {\kappa }_{0z} $表示弯曲曲率.梁中心线的格林--拉格朗日应变张量可以代入方程(30),但是该方法求解的弹性力不包含梁单元的弯曲和扭转, 弯曲曲率$\bar{\kappa }_{0y} ,\bar {\kappa }_{0z} $用于求解弯曲应变能,而扭转转角$\theta $用于表达扭转应变能,弯曲应变能和扭转变形能的表达式为

(33)$$U_{\rm b} = \frac{1}{2}\int_0^l {[EI_y (\bar {\kappa }_{0y} )^2 + EI_z (\bar {\kappa }_{0z} )^2]} {\rm d}x \\ \ U_{\rm t} = \frac{1}{2}\int_0^l {GI_p \theta ^2} {\rm d}x $$

则广义的弹性力的矢量$Q_k $为

(34)\begin{equation} \label{eq32} Q_k =-\int_v {\left(\frac{\partial \varepsilon _{0} }{\partial e}\right)} ^{\rm T} E\varepsilon _{0} \mbox{d }v-\left(\frac{\partial U_{\rm b} }{\partial e} + \frac{\partial U_{\rm t} }{\partial e}\right) \end{equation}

其中, $\varepsilon _{0} $为弹性线方法中格林--拉格朗日应变张量,具体见式(32).

3.3 两种方法的比较

为了比较一般连续介质力学方法和弹性线的应变张量, 全局位移矢量$r$可写为

(35)\begin{equation}\label{eq33} r = r_{0} + y r_y + z r_z \end{equation}

其中, $r_{0} $为中心线上任意一点的位移矢量, $\Delta r = y r_y + z r_z $表示梁中心线横截面上任意一点的位移矢量.则全局位移矢量对于空间坐标$x$的导数$r_x $为

(36)\begin{equation}\label{eq34} r_x = r_{0x} + y r_{yx} + z r_{zx}\end{equation}

采用如前所示的位移场和格林--拉格朗日应变张量,基于连续介质力学方法得到梁的应变分量为

(37)$$\varepsilon _{11} = \frac{1}{2}( r_x^{\rm T} r_x-1) = \frac{1}{2}( r_{0x}^{\rm T} r_{0x}-1) + y r_{0x}^{\rm T} r_{yx} +\\ z r_{0x}^{\rm T} r_{zx} + \frac{1}{2}(y r_{yx} + z r_{zx} )^{\rm T}(y r_{yx} + z r_{zx} ) \\ \ \varepsilon _{22} = \frac{1}{2}( r_y^{\rm T} r_y-1) ,\quad \varepsilon _{33} = \frac{1}{2}( r_z^{\rm T} r_z-1) \\ \ \gamma _{12} = r_x^{\rm T} r_y = r_{0x}^{\rm T} r_y + y r_{yx}^{\rm T} r_y + z r_{zx}^{\rm T} r_y \\ \ \gamma _{13} = r_x^{\rm T} r_z = r_{0x}^{\rm T} r_z + y r_{yx}^{\rm T} r_z + z r_{zx}^{\rm T} r_z \\ \ \gamma _{23} = r_y^{\rm T} r_z $$

应变张量可转化为

(38)$$\varepsilon _{11} = \frac{1}{2}( r_{0x}^{\rm T} r_{0x}-1)-y\bar {\kappa }_3 + z\bar {\kappa }_2 +\\ \frac{1}{2}(y r_{yx} + z r_{zx} )^{\rm T}(y r_{yx} + z r_{zx} ) \cong \varepsilon _{011}-y\bar {\kappa }_3 + z\bar {\kappa }_2 \\ \ \gamma _{12} = r_{0x}^{\rm T} r_y-z\bar {\kappa }_1 \equiv \gamma _{012}-z\bar {\kappa }_1 \\ \ \gamma _{13} = r_{0x}^{\rm T} r_z + y\bar {\kappa }_1 \equiv \gamma _{013} + y\bar {\kappa }_1 \\ \ \varepsilon _{22} = \frac{1}{2}( r_y^{\rm T} r_y-1) ,\quad \gamma _{23} = r_y^{\rm T} r_z ,\quad \varepsilon _{33} = \frac{1}{2}( r_z^{\rm T} r_z-1) \\ \ \bar {\kappa }_1 = \frac{1}{2}({ r}'^{\rm T}_y r_z- r_y^{\rm T} { r}'_z ) \equiv \theta \\ \ \bar {\kappa }_2 = \frac{ r_y^{\rm T} r_{0xx} }{\left| { r_y } \right|\left| { r_{0xx} } \right|} =-\bar {\kappa }_{0z} \approx r_y^{\rm T} r_{0xx} \\ \ \bar {\kappa }_3 = \frac{ r_z^{\rm T} r_{0xx} }{\left| { r_z } \right|\left| { r_{0xx} } \right|} = \bar {\kappa }_{0y} \approx r_z^{\rm T} r_{0xx} $$

比较式(37)和式(38),可知采用一般连续介质力学方法和弹性线方法在描述变形上的差异有：

(1) 在一般连续介质力学方法中, 弯曲应变是通过$y r_{0x}^{\rm T} r_{yx} + z r_{0x}^{\rm T} r_{zx} $来定义的;而弹性线方法的弯曲应变是通过弯曲曲率$\bar {\kappa }_{0y} $和$\bar{\kappa }_{0z} $来定义的. 这两种定义方式是不一样的,并可能带来不同的刚度特性.

(2) 在弹性线理论方法中, 忽\vspace{1mm}略了对弯曲变形运动影响相对很小的高阶正应变分量$\dfrac{1}{2}(y r_{yx} + z r_{zx} )^{\rm T}\cdot$$(y r_{yx} + z r_{zx} )$.

(3) 在一般连续介质力学理论方法中, 扭转应变通过$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 $来定义表达;而弹性线方法通过转角表达式$\theta $来表达定义.这两种定义方式是不一样的, 并可能带来不同的刚度特性.

4 动力学仿真

4.1 柔性单摆动力学仿真

物理模型如图2所示,柔性单摆分别受到集中力、均匀分布力以及外力矩作用,自由落体柔性单摆(重力加速度为$g = 9.81$m/s$^{2}$)、受外力矩作用柔性单摆($M = 10$Nm)以及受外力作用柔性单摆($F = 10$N). 柔性单摆物理参数如表1所示.

图2

新窗口打开|下载原图ZIP|生成PPT
图2柔性单摆物理模型

Fig. 2Physical model of the flexible pendulum



Table 1
表1
表1柔性单摆物理参数
Table 1Physical parameter of the flexible pendulum

新窗口打开|下载CSV

图3和图4分别表示自由落体柔性单摆末端$B$点$X,Y$方向位移;图5和图6分别表示柔性单摆在外力矩作用下末端$B$点$X$方向、$Y$方向的位移;图7和图8分别表示柔性单摆在外力作用下末端$B$点$X$方向、$Y$方向的位移.从图3$\sim $\!图8可知两种方法的仿真结果基本一致,说明弹性线方法能在简化的同时满足计算需要.表2给出两种方法的计算时间, 从表中可知,弹性线方法能有效提高绝对节点坐标法的计算效率.图9表示三种受力情况下柔性单摆某些时刻的整体位形,从图中可知, 柔性单摆的变形较大.

Table 2
表2
表2两种方法的计算时间(min)
Table 2The calculating time of two methods

新窗口打开|下载CSV

图3

新窗口打开|下载原图ZIP|生成PPT
图3自由落体柔性单摆末端$B$点$X$方向的位移

Fig. 3The displacement of the $X$direction at the end of the free-fall flexible pendulum

图4

新窗口打开|下载原图ZIP|生成PPT
图4自由落体柔性单摆末端$B$点$Y$方向的位移

Fig. 4The displacement of the $Y$direction at the end of the free-fall flexible pendulum

图5

新窗口打开|下载原图ZIP|生成PPT
图5受外力矩作用柔性单摆末端$B$点$X$方向位移

Fig. 5The displacement of the $X$direction at the end of the flexible pendulum which subject to external torque

图6

新窗口打开|下载原图ZIP|生成PPT
图6受外力矩作用柔性单摆末端$B$点$Y$方向位移

Fig. 6The displacement of the $Y$direction at the end of the flexible pendulum which subject to external torque

图7

新窗口打开|下载原图ZIP|生成PPT
图7受外力作用柔性单摆末端$B$点$X$方向位移

Fig. 7The displacement of the $X$direction at the end of the flexible pendulum which subject to external force

图8

新窗口打开|下载原图ZIP|生成PPT
图8受外力作用柔性单摆末端$B$点$Y$方向位移

Fig. 8The displacement of the $X$direction at the end of the flexible pendulum which subject to external force

图9

新窗口打开|下载原图ZIP|生成PPT
图9柔性单摆整体位置

Fig. 9The configuration of the flexible pendulum

4.2 履带式车辆动力学仿真

考虑如图10所示的履带式车辆,本仿真采用三维梁单元对该模型两条链条进行建模, 而其他结构设为刚体,链条的材料为不锈钢S15600, 其材料参数为：密度$\rho =15.6~\mbox{Mg/m}^{3}$, 弹性模量$E =200~\mbox{GPa}$, 泊松比$\mu =0.3$, 重力加速度$g = 9.81~\mbox{m/s}^{2}$.每条链条被分成64个三维梁单元, 由于链条首尾相连,即第一个节点亦为最后一个节点, 则每条链条的单元节点为64个,每个单元的长度$l = 0.151~624~\mbox{m}$, 单元宽度$b =0.381~\mbox{m}$, 厚度$h = 0.02~\mbox{m}$.全局坐标系为：履带式车辆的前进方向为$Y$轴, 沿地面向上为$Z$轴,则$X$轴沿图10向内. 地面为刚体, 链条与地面的接触采用Hertz接触.采用连续介质力学方法和弹性线方法对履带式车辆进行仿真,定义链轮的旋转角速度(单位： rad/s)为如下函数形式

(39)\begin{equation} \label{eq37} \omega \!=\! \left\{\!\!\! {\begin{array}{ll} 0, \hspace{-3mm}& 0 < t \le 1~\mbox{s} \\[2mm] \dfrac{80}{3}(t-1)^2 , \hspace{-3mm}& 1~\mbox{s} < t \le 1.25~\mbox{s} \\[4mm] \dfrac{40}{3}(t \!-\! 1.25) \!+\! \dfrac{5}{3} , \hspace{-3mm}& 1.25~\mbox{s} \!<\! t \!\le\! 1.75~{\rm s} \\[4mm] \dfrac{40}{3}(t \!-\! 1.75) \!-\! \dfrac{80}{3}(t \!-\! 1.75)^2 \!+\! \dfrac{25}{3}, \hspace{-3mm}& 1.75~\mbox{s} \!<\! t \!\le\! 2~\mbox{s} \\[2mm] 10 , \hspace{-3mm}& 2~\mbox{s} < t \le 5~\mbox{s} \\ \end{array}} \right. \end{equation}

图10

新窗口打开|下载原图ZIP|生成PPT
图10履带式车辆模型

Fig. 10The model of a tracked vehicle



图11表示履带式车辆运动$X,Y$和$Z$方向的位移,Z方向有位移波动是因为履带式车辆在Z方向受到重力作用.图12表示履带式车辆前进方向的运动速度. 从图11和图12可知,弹性图11表示履带式车辆运动$X,Y$和$Z$方向的位移,Z方向有位移波动是因为履带式车辆在Z方向受到重力作用.图12表示履带式车辆前进方向的运动速度.从图11和图12可知, 弹性线方法的仿真结果与连续介质力学方法的基本一致,说明采用弹性线方法来求解弹性力的方法是可行且正确的.同时从计算时间上考虑,采用连续介质力学方法求解弹性力的履带式车辆的计算时间是采用弹性线方法的9倍.仿真结果表明：采用弹性线方法求解弹性力的履带式车辆仿真在保证计算精度的情况下,计算效率将远高于连续介质力学方法的.因此弹性线方法值得在绝对节点坐标法中重点研究,其将有效弥补绝对节点坐标法计算效率较低的缺点.

图11

新窗口打开|下载原图ZIP|生成PPT
图11履带式车辆运动位移

Fig. 11The displacement of the tracked vehicle

图12

新窗口打开|下载原图ZIP|生成PPT
图12履带式车辆前进速度

Fig. 12The forward speed of the tracked vehicle

5 结论

(1)弹性线方法将弯曲变形和扭转应变通过曲率和转角来表示,将弹性力求解从三重体积分简化到沿梁长方向的一重积分,其大大节省仿真计算时间.

(2)通过柔性单摆的动力学仿真表明：弹性线的仿真结果与一般连续介质力学方法的一致,且仿真计算效率要优于一般连续介质力学方法的.

(3)采用弹性线方法能有效处理实际工程履带式车辆的动力学问题,并在保证精度的要求下提高计算效率.

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

[1] Rafiee M, Nitzsche F, Labrosse M . Dynamics, vibration and control of rotating composite beams and blades: A critical review
Thin-Walled Structures, 2017,119:795-819
[本文引用: 1]

[2] 曹登庆, 白坤朝, 丁虎 等. 大型柔性航天器动力学与振动控制研究进展
力学学报, 2019,51(1):1-13
[本文引用: 1]

( Cao Dengqing, Bai Kunchao, Ding Hu , et al. Advances in dynamics and vibration control of large-scale flexible spacecraft
Chinese Journal of Theoretical and Applied Mechanics, 2019,51(1):1-13 (in Chinese))
[本文引用: 1]

[3] Liu JY, Hong JZ . Geometric stiffening effect on rigid-flexible coupling dynamics of an elastic beam
Journal of Sound and Vibration, 2004,278(4-5):1147-1162
[本文引用: 1]

[4] 高晨彤, 黎亮, 章定国 等. 考虑剪切效应的旋转 FGM 楔形梁刚柔耦合动力学建模与仿真
力学学报, 2018,50(3):654-666


( Gao Chentong, Li Liang, Zhang Dingguo , et al. Dynamic modeling and simulation of rotating FGM tapered beams with shear effect
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(3):654-666 (in Chinese))


[5] Cai GP, Hong JZ, Yang SX . Dynamic analysis of a flexible hub-beam system with tip mass
Mechanics Research Communications, 2005,32(2):173-190


[6] 和兴锁, 李雪华, 邓峰岩 . 平面柔性梁的刚--柔耦合动力学特性分析与仿真
物理学报, 2011,60(2):024502


( He Xingsuo, Li Xuehua, Deng Fengyan . Analysis and imitation of dynamic properties for rigid-flexible coupling systems of a planar flexible beam
Acta Physica Sinica, 2011,60(2):024502 (in Chinese))


[7] 陈思佳, 章定国, 洪嘉振 . 大变形旋转柔性梁的一种高次刚柔耦合动力学模型
力学学报, 2013,45(2):251-256


( Chen Sijia, Zhang Dingguo, Hong Jiazhen . A high-order rigid-flexible coupling model of a rotating flexible beam under large deformation
Chinese Journal of Theoretical and Applied Mechanics, 2013,45(3):251-256 (in Chinese))


[8] 范纪华, 章定国 . 旋转柔性悬臂梁动力学的Bezier 插值离散方法研究
物理学报, 2014,63(15):154501


( Fan Jihua, Zhang Dingguo . Bezier interpolation method for the dynamics of rotating flexible cantilever beam
Acta Physica Sinica, 2014,63(15):154501 (in Chinese))


[9] 杜超凡, 章定国, 洪嘉振 . 径向基点插值法在旋转柔性梁动力学中的应用
力学学报, 2015,47(2):279-288


( Du Chaofan, Zhang Dingguo, Hong Jiazhen . A meshfree method based on radial point interpolation method for the dynamic analysis of rotating flexible beams
Chinese Journal of Theoretical and Applied Mechanics, 2015,47(2):279-288 (in Chinese))


[10] 范纪华, 章定国 . 基于变形场不同离散方法的柔性机器人动力学建模与仿真
力学学报, 2016,48(4):843-856


( Fan Jihua, Zhang Dingguo . Dynamic modeling and simulation of flexible robots based on different discretization methods
Chinese Journal of Theoretical and Applied Mechanics, 2016,48(4):843-856 (in Chinese))


[11] 章孝顺, 章定国, 洪嘉振 . 考虑曲率纵向变形效应的大变形柔性梁刚柔耦合动力学建模与仿真
力学学报, 2016,48(3):692-701


( Zhang Xiaoshun, Zhang Dingguo, Hong Jiazhen . Rigid-flexible coupling dynamic modeling and simulation with the longitudinal deformation induced curvature effect for a rotating flexible beam under large deformation
Chinese Journal of Theoretical and Applied Mechanics, 2016,48(3):692-701 (in Chinese))


[12] Zhao G, Wu Z . Coupling vibration analysis of rotating three-dimensional cantilever beam
Computers & Structures, 2017,179:64-74


[13] Sujash B, Debabrata D . Free vibration analysis of bidirectional-functionally graded and doubletapered rotating micro-beam in thermal environment using modified couple stress theory
Composite Structures, 2019,215:471-492


[14] Tian JJ, Zhang ZG, Hua HX . Free vibration analysis of rotating functionally graded double-tapered beam including porosities
International Journal of Mechanical Sciences, 2019,150:526-538


[15] Fang JS, Wang HW, Zhang XP . On size-dependent dynamic behavior of rotating functionally graded Kirchhoff microplates
International Journal of Mechanical Sciences, 2019,152:34-50
[本文引用: 1]

[16] Shabana AA. Computational Continuum Mechanics, New York: Cambridge University Press, 2008
[本文引用: 2]

[17] Bonet J, Wood RD. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge: Cambridge University Press, 1997
[本文引用: 1]

[18] Simo JC, Vu-Quoc L . On the dynamics in space of rods undergoing large motions the planar case: Part I-II. ASME Transactions,
Journal of Applied Mechanics, 1986,53:849-863
[本文引用: 1]

[19] Berzeri M, Shabana AA . Development of simple models for the elastic forces in the absolute nodal coordinate formulation
Journal of Sound and Vibration, 2000,235(4):539-565
[本文引用: 2]

[20] Omar MA, Shabana AA . A two-dimensional shear deformable beam for large rotation and deformation problems
Journal of Sound and Vibration, 2001,243(3):565-576
[本文引用: 1]

[21] Dufva KE, Sopanen JT, Mikkola AM . A two-dimensional shear deformation beam element based on the absolute nodal coordinate formulation
Journal of Sound and Vibration, 2005,280:719-738
[本文引用: 1]

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

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

[24] Tian Q, Zhang Y, Chen L , et al. Simulation of planar flexible multibody systems with clearance and lubricated revolute joints
Nonlinear Dynamics, 2010,60(4):489-511
[本文引用: 1]

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

[26] 过佳雯, 魏承, 谭春林 等. 含芯拧绞绳非线性弯曲动力学特性分析与研究
力学学报, 2018,50(2):373-384
[本文引用: 1]

( Guo Jiawen, Wei Cheng, Tan Chunlin , et al. Analysis of the cored stranded wire rope on the nonlinear bending dynamic characteristics
Chinese Journal of Theoretical and Applied Mechanics, 2018,50(2):373-384 (in Chinese))
[本文引用: 1]

[27] 兰朋, 崔雅琦, 於祖庆 . 完备化ANCF 薄板单元及在钢板弹簧动力学建模中的应用
力学学报, 2018,50(5):1156-1167
[本文引用: 1]

( 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))
[本文引用: 1]

[28] 郑彤, 章定国, 洪嘉振 . 三维大变形梁系统的动力学建模与仿真
机械工程学报, 2016,52(19):81-87
[本文引用: 2]

( Zheng Tong, Zhang Dingguo, Hong Jiazhen . Dynamic modeling and simulation for three dimensional flexible beam systems with large deformations
Journal of Mechanical Engineering, 2016,52(19):81-87 (in Chinese))
[本文引用: 2]

[29] 赵春璋, 余海东, 王皓 等. 基于绝对节点坐标法的变截面梁动力学建模与运动变形分析
机械工程学报, 2014,50(17):38-45
[本文引用: 1]

( Zhao Chunzhang, Yu Haidong, Wang Hao , et al. Dynamic modeling and kinematic behavior of variable cross-section beam based on the absolute nodal coordinate formulation
Journal of Mechanical Engineering, 2014,50(17):38-45 (in Chinese))
[本文引用: 1]

[30] Cui YQ, Yu ZQ, Lan P . A novel method of thermo-mechanical coupled analysis based on the unified description
Mechanism and Machine Theory, 2019,134:376-392
[本文引用: 1]

[31] Liu J, Pan K . Rigid-flexible-thermal coupling dynamic formulation for satellite and plate multibody system
Aerospace Science and Technology, 2016,52:102-114
[本文引用: 1]

[32] Shen Z, Li H, Liu X , et al. Thermal shock induced dynamics of a spacecraft with a flexible deploying boom
Acta Astronautica, 2017,141:123-131
[本文引用: 1]

[33] García-Vallejo D, Mayo J, Escalona JL , et al. Efficient evaluation of the elastic forces and the Jacobian in the absolute nodal coordinate formulation
Nonlinear Dynamics, 2004,35(4):313-329
[本文引用: 2]

[34] 章孝顺, 章定国, 陈思佳 等. 基于绝对节点坐标法的大变形柔性梁几种动力学模型研究
物理学报, 2016,65(9):094501
[本文引用: 2]

( Zhang Xiaoshun, Zhang Dingguo, Chen Sijia , et al. Several dynamic models of a large deformation flexible beam based on the absolute nodal coordinate formulation
Acta Physica Sinica, 2016,65(9):094501 (in Chinese))
[本文引用: 2]

[35] 陈渊钊, 章定国, 黎亮 . 平面细长梁基于无网格径向基点插值的绝对节点坐标法
振动工程学报, 2018,31(2):245-254
[本文引用: 1]

( Chen Yuanzhao, Zhang Dingguo, Li Liang . An absolute nodal coordinate formation based on radial point interpolation method for planar slender beams
Journal of Vibration Engineering, 2018,31(3):245-254 (in Chinese))
[本文引用: 1]

[36] Chen Y, Zhang D, Li L . Dynamic analysis of rotating curved beams by using absolute nodal coordinate formulation based on radial point interpolation method
Journal of Sound and Vibration, 2019,441:63-83
[本文引用: 1]

[37] Zemljari$\check{\rm o}$B, A$\check{\rm z}$be V . Analytically derived matrix end-form elastic-forces equations for a low-order cable element using the absolute nodal coordinate formulation
Journal of Sound and Vibration, 2019,446:263-272
[本文引用: 1]

[38] Luo K, Hu H, Liu C , et al. Model order reduction for dynamic simulation of a flexible multibody system via absolute nodal coordinate formulation
Computer Methods in Applied Mechanics and Engineering, 2017,324:573-594
[本文引用: 2]