引言

绝对节点坐标法(absolute nodal coordinate formulation, ANCF)于1996年由Shabana提出^[1]. 由于采用了全局坐标系下定义的位置和梯度向量作为单元节点坐标, 从而使得梁^[2]、 板、壳等在传统有限元描述下的非等参单元都具有了等参单元的性质, 可以用更少数目的单元刻画建模对象的复杂几何形状. 同时, 导出的动力学方程具有常值质量阵、方程中不含有科氏加速度和离心加速度项等优点. 非线性的格林‒拉格朗日应变张量的采用使得绝对节点坐标法可以描述柔性体的大范围变形, 并且能保证在任意的刚体转动下不产生应变, 因此绝对节点坐标法适用于大转动、大变形的非线性动力学问题^[3]. 目前, 绝对节点坐标法已经有了较为完备的单元库, 在解决工程实际问题中得到了日益广泛的应用^[4-7].

现有绝对节点坐标法单元分为全参数单元与梯度不完备单元, 其中全参数单元节点处均采用位置和沿长、宽、高3个方向的梯度向量作为节点坐标. 因此, 单元内部任意点处均可求取完备的梯度向量, 从而利用一般连续介质力学方法计算单元弹性力^[8-9], 同时能够系统地消除初始弯曲构型中的初应变^[10]. 然而, 由于在各个方向上采用了不同阶数的插值多项式, 全参数梁单元^[11]与板单元^[12]存在的截面闭锁问题会降低动力学仿真计算效率. 而梯度不完备单元如薄板单元^[13]和缆索单元等, 在计算某些特定结构时具有较高的仿真效率. 可同时, 由于在单元内部点处不具备完备的梯度向量, 导致梯度向量矩阵不是方阵, 无法采用求逆的方式消除初始弯曲结构中的初始应变. 为此, 本文从薄板单元基本假定出发, 改进其弹性力模型推导完备化的薄板单元, 将其与一般连续介质力学描述方法整合, 从而消除初始应变.

钢板弹簧是一种广泛应用于多种客、货运车辆悬架系统中的弹性支撑部件, 用于连接汽车底盘与车轴, 在车辆运行时为车架提供稳定的支撑^[12]. 它的动力学响应特性决定了车辆的操控性、舒适性、振动特性以及稳定性等关键性能. 然而, 由于钢板弹簧较复杂的结构、簧片本身的材料非线性^[14]、复杂的板间接触和摩擦等问题的存在, 钢板弹簧是一个高度非线性的系统^[15]. 对典型的串联悬挂系统的实验结果分析也表明, 钢板弹簧除了具有传统弹簧有效刚度的变化特性, 还具有簧片间库伦摩擦力导致的能量衰减和滞回现象. 因此, 钢板弹簧的动力学仿真是典型的大变形几何非线性动力学问题. 常用的分析模型有简化物理模型和传统浮动坐标法(floating frame reference, FFR)有限元模型. 简化模型包括美国汽车工程学会提出的三连杆模型^[16]、 变截面梁单元模型^[17]以及离散梁模型^[18]. 简化物理模型一般通过系数调整以期达到较好的近似, 实际上, 无论是简化物理模型还是传统的FFR有限元模型, 都未能很好地刻画钢板弹簧复杂的动力学特性, 并且不具有较好的概括性和适应性^[19]. Omar等^[11]分别采用浮动坐标法梁单元与绝对节点坐标全参数梁单元建立了车辆钢板弹簧模型. 通过对比得出了浮动坐标法容易实现自由度缩减而绝对节点坐标法更适合刻画簧片的大变形结论. 同时采用了给定节点的接触模型. Yu等^[12]则建立了基于全参数板单元的钢板弹簧模型, 并通过引入参考节点的概念, 用一个网格建立起车架刚柔耦合模型. 而全参数单元在使用中由于存在截面闭锁问题, 会导致求解效率降低. 本文拟采用绝对节点坐标薄板单元建立钢板弹簧模型, 以提高仿真分析效率.

钢板弹簧系统中存在的接触与摩擦作为多柔体问题动力学建模中难度较大的一环, 吸引了很多相关研究. Wriggers和Zavarise^[20]首先提出了圆截面梁之间的逐点无摩擦接触模型之后, 又进一步给出了基于库伦摩擦模型计算摩擦力的回归映射算法^[21]. 基于Wriggers和Zavarise提出的模型, Litewka和Wriggers ^[22]进一步给出了方截面接触模型. Litewka提出了一种基于艾尔米特多项式的曲线光滑方法^[23]并用在了Wriggers和Zavarise提出的梁接触模型上, 紧接着又提出了类比塑性的摩擦力计算模型^[24]. Konyukhov和Schweizerhof ^[25]提出了曲线与曲线接触计算的协变形式, 在曲线坐标系中完成投影距离最近点的计算. 这些研究都是建立在假设曲线间的接触是逐点发生的基础之上的. Litewka ^[26]在投影最近点两侧增加两对计算点建立了多点接触计算方法后, Durville 用中间几何体的方法建立了圆截面空间梁的有摩擦接触模型^[27]. 在中间几何体概念的基础上, Wang^[28-29] 等先后建立了基于绝对节点坐标缆索单元的单接触域和多接触域的接触与摩擦施加方法. 针对接触这一钢板弹簧建模中最为复杂的问题^[30], 本文将给出新的钢板弹簧板间接触模型, 通过建立簧片局部坐标系给出基于主从方法的接触副跨单元搜索策略, 并利用惩罚函数法与平滑化的库伦摩擦模型推导较为精确的接触力与摩擦力计算格式. 还将通过动力学仿真算例给出的结果, 分析过程中的动力学现象.

1 单元运动学描述

一般连续介质力学方法采用非线性的格林拉格朗日应变张量描述柔性体的变形, 定义为: $\textbf{ε}=({J}^{\rm T}{J}-{I})/2$. 其中梯度向量矩阵${J}=\partial {r}/\partial {X}$ 为当前构型下位置坐标对参考构型下点坐标取导数所得的雅可比矩阵, 如图1所示. 由于绝对节点坐标方法中所有的直接求导计算只能针对母构型进行^[31], 因此在参考构型带有初始弯曲时, 梯度向量的计算必须通过母构型利用链式求导法则进行

\begin{equation*}{J}=\dfrac{\partial {r}}{\partial X}=\dfrac{\partial {r}}{\partial {x}}\cdot\dfrac{\partial{x}}{\partial X}={J}_e\cdot{J}_0^{-1}\tag*{(1)}\end{equation*}

参考构型带有初始弯曲时, 格林拉格朗日应变张量的定义为

\begin{equation*}\textbf{ε}=\dfrac{1}{2}({J}^{\rm T}{J}-{I})=\dfrac{1}{2} ({J}_0^{-\rm T}{J}_e^{\rm T}\cdot{J}_e{J}_0^{-1}-{I})\tag*{(2)}\end{equation*}

其中${J}_e=\partial({S}\cdot{e})/\partial{x}$, 而${J}_0=\partial({S}\cdot{e}_0)/\partial{x}$只由初始时刻的节点坐标决定. 可以看出, 在仿真开始的时刻${e}={e}_0$, 则将${J}_e={J}_0$代入将得到全零的应变张量, 消除了参考构型带有初始弯曲时原本会引入的预应力. ${J}_0$必须是方阵才能求逆, 所以这就要求单元必须是全参数的, 节点具有完备的3个方向上的梯度向量, 薄板等非全参数单元则无法直接使用该方法.

另一方面, 绝对节点坐标全参数板单元在长度和宽度方向采用三次多项式插值, 而厚度方向则采用了一次多项式. 由此造成单元发生弯曲变形时, 在厚度方向产生不符合实际的附加变形. 而各向同性材料本构关系的应用会在厚度方向产生高频震荡, 造成单元闭锁, 从而迫使系统动力学方程采用较小的积分时间步长, 降低了求解效率.

薄板单元假定单元截面厚度保持不变, 并且时刻与中面保持垂直, 单元包含中面上的4个节点, 每个节点选择位置向量和长、宽两个方向的梯度向量作为节点坐标向量^[13]

\begin{equation*}\left.\begin{split}&{r}={S}\cdot{e}\\&{e}=\left\lbrack {r}_1^{\rm T}, {r}_{x1}^{\rm T}, {r}_{y1}^{\rm T}, {r}_2^{\rm T}, {r}_{x2}^{\rm T}, {r}_{y2}^{\rm T}, {r}_3^{\rm T}, {r}_{x3}^{\rm T}, {r}_{y3}^{\rm T}, {r}_4^{\rm T}, {r}_{x4}^{\rm T}, {r}_{y4}^{\rm T}\right\rbrack^{\rm T}\end{split}\right\}\tag*{(3)}\end{equation*}

薄板单元消除了全参数板单元在厚度方向上产生的闭锁现象, 在描述薄壳状柔性体变形时相比全参数板单元具有更高的计算效率^[13]. 本文采用薄板单元建立的钢板弹簧模型如图2所示.

现有薄板单元应变能由单元中面的拉压、剪切和中面的弯曲两部分构成

\begin{equation*}U=\dfrac{1}{2}\int_V \textbf{ε}^{\rm T} {E} \textbf{ε} {\rm d}V+\dfrac{1}{2} \int_V {\kappa}^{\rm T}{E}\kappa{\rm d}V\tag*{(4)}\end{equation*}

其中$ E$为单元刚度系数阵, 格林-拉格朗日应变张量表征单元中面的拉压和剪切, 向量形式$\textbf{ε}$记做

\begin{equation*}\textbf{ε}=\left\lbrack \varepsilon_{xx}\varepsilon_{yy}\varepsilon_{xy}\right\rbrack^{\rm T}\tag*{(5)}\end{equation*}

${\kappa}$表示的是用单元中面曲率变化所表征的弯曲变形

\begin{equation*}\left.\begin{split}&{\kappa}=z\left\lbrack \kappa_{xx}, \kappa_{yy}, \kappa_{xy}\right\rbrack^{\rm T}\\&\kappa_{xx}=\dfrac{{r}_{xx}^{\rm T}\cdot{n}}{\big|\big|{n}\big|\big|^3}, \kappa_{yy}=\dfrac{{r}_{yy}^{\rm T}\cdot{n}}{\big|\big|{n}\big|\big|^3}, \kappa_{xy} \dfrac{{r}_{xy}^{\rm T}\cdot{n}}{\big|\big|{n}\big|\big|^3}\end{split}\right\}\tag*{(6)}\end{equation*}

其中, 中面单位法向量$ n$的定义式为

\begin{equation*}{n}=\dfrac{{r}_x\times{r}_y}{\big|\big|{r}_x\times{r}_y\big|\big|}\tag*{(7)}\end{equation*}

通过应变能对于节点坐标求导, 得到单元的弹性力${Q}_{\rm e}$, 由代表单元中面的拉压和剪切的${Q}_{\rm s}$ 和表示单元中面弯曲的${Q}_\text{κ}$两项组成

\begin{equation*}\begin{split}{Q}_{\rm e}=&\dfrac{\partial U}{\partial{e}}={Q}_{\rm s}+{Q}_\text{κ}=\\&\int_V\dfrac{\partial\textbf{ε}^{\rm T}}{\partial{e}}{E}\textbf{ε}{\rm d}V+\int_V\dfrac{\partial{\kappa}^{\rm T}}{\partial{e}}{E\kappa}{\rm d}V\end{split}\tag*{(8)}\end{equation*}

可以看出, 现有薄板单元缺少沿厚度方向的梯度向量$\partial{r}/\partial z$, 导致梯度向量矩阵不是方阵, 无法对带有初始弯曲的参考构型的梯度向量矩阵求逆, 弹性力模型无法将参考构型的初始弯曲引入的初应变消除掉. 若将其直接应用于钢板弹簧的动力学仿真, 会出现簧片自动回弹至平直状态的现象. 这就需要提出系统性解决这一问题的方法.

2 完备化的薄板单元弹性力模型

由前文的分析可知, 由于薄板单元直接建立初始弯曲构型时带有初始应变. 为此, 将单元中面单位法向量$ n$ 作为厚度方向的梯度向量写入梯度向量矩阵

\begin{equation*}{J}=\left\lbrack\dfrac{\partial{r}}{\partial x}\dfrac{\partial{r}}{\partial y}{n}\right\rbrack\tag*{(9)}\end{equation*}

可知其始终保持与单元中面所在点的其他两个梯度向量垂直, 并且模长始终为1. 这样就可以采用处理全参数单元的方式对薄板单元进行去除初始应变的操作, 得到的格林-拉格朗日应变张量记为

\begin{equation*}\begin{split}\textbf{ε}_v=&\left\lbrack \varepsilon_{XX}\varepsilon_{YY}\varepsilon_{ZZ}\varepsilon_{XY}\varepsilon_{XZ}\varepsilon_{YZ}\right\rbrack^{\rm T}=\\&\left\lbrack \dfrac{1}{2}({r}_X^{\rm T}{r}_X-1)\dfrac{1}{2}({r}_Y^{\rm T}{r}_Y-1)\dfrac{1}{2}({n}^{\rm T}{n}-1)\right.\\&\left.\dfrac{1}{2}({r}_X^{\rm T}{r}_Y)\dfrac{1}{2}({r}_X^{\rm T}{n})\dfrac{1}{2}({r}_Y^{\rm T}{n})\right\rbrack^{\rm T}\end{split}\tag*{(10)}\end{equation*}

由于单位法向量$ n$的定义方式, 应变张量中与厚度相关的分量$\varepsilon_{ZZ}$, $\varepsilon_{XZ}$ 和$\varepsilon_{YZ}$始终为零.

完备化的薄板单元弹性力公式为

\begin{equation*}\begin{split}{Q}_{\rm e}=&{Q}_{\rm s}+{Q}_\text{κ}=\\&\int_V \dfrac{\partial \textbf{ε}_v^{\rm T}}{\partial{e}}{E}_v \textbf{ε}_v{\rm d}V+\int_V\dfrac{\partial \kappa_v^{\rm T}}{\partial {e}}{E} \kappa_v{\rm d}V\end{split}\tag*{(11)}\end{equation*}

式中${E}_v$为对应于格林-拉格朗日应变张量$\textbf{ε}_v$的单元刚度系数阵. 格林-拉格朗日应变张量中厚度方向有关项对于节点坐标的导数为

\begin{equation*}\left.\begin{split}&\dfrac{\partial \textbf{ε}_{XZ}}{\partial e_i}=\dfrac{1}{2}\dfrac{\partial ({r}_X^{\rm T}\cdot{n})}{\partial e_i}=\dfrac{1}{2}\left\lgroup {S}_{X,i}^{\rm T}\cdot {n}+{r}_X^{\rm T}\cdot \dfrac{\partial {n}}{\partial e_i}\right\rgroup\\&\dfrac{\partial \textbf{ε}_{YZ}}{\partial e_i}=\dfrac{1}{2}\dfrac{\partial ({r}_Y^{\rm T}\cdot {n})}{\partial e_i}=\dfrac{1}{2} \left\lgroup {S}_{Y,i}^{\rm T}\cdot{n}+{r}_Y^{\rm T}\cdot\dfrac{\partial {n}}{\partial e_i}\right\rgroup\\&\dfrac{\partial \textbf{ε}_{ZZ}}{\partial e_i}=\dfrac{1}{2}\dfrac{\partial ({n}^{\rm T}{n}-1)}{\partial e_i}=\dfrac{\partial{n}^{\rm T}}{\partial e_i}\cdot{n}\end{split}\right\}\tag*{(12)}\end{equation*}

这里令${g}={r}_x\times{r}_y$, 则单位法向量$ n$对于节点坐标的导数项表示为

\begin{equation*}\begin{split}\dfrac{\partial {n}}{\partial e_i}=&\partial\left\lbrack {g}\cdot\left({g}^{\rm T}\cdot{g}\right)^{-\tfrac{1}{2}}\right\rbrack\bigg/\partial e_i=\\&\dfrac{\partial {g}}{\partial e_i}\cdot\left({g}^{\rm T}\cdot {g}\right)^{-\tfrac{1}{2}}-{g}\cdot\left({g}^{\rm T}\cdot {g}\right)^{-\tfrac{3}{2}}\cdot\dfrac{\partial {g}^{\rm T}}{\partial e_i}\cdot{g}\quad\end{split}\tag*{(13)}\end{equation*}

式(12)中, ${S}_{X,i}$代表薄板单元形函数$ S$对物质坐标$X$的导数${S}_X$的第$i$列. 值得注意的是, 虽然通过添加单位法向量的方法将薄板单元纳入到一般连续介质力学的表达范式中, 但薄板单元不存在厚度方向的物质点坐标, 求取单元弹性力时仍然只做中面上的积分, 这样就无法引入抗弯刚度. 因此, 现有薄板单元弹性力中的弯曲项仍然需要保留, 并用当前构型下的弯曲向量减去参考构型下弯曲向量以消除参考构型中的初始弯曲应变

\begin{equation*}{\kappa}_v=z\left\lbrack \kappa_{xx}-\kappa_{xx0}, \kappa_{yy}-\kappa_{yy0}, \kappa_{xy}-\kappa_{xy0}\right\rbrack^{\rm T}\tag*{(14)}\end{equation*}

为采用隐式积分器求解系统动力学方程, 需要用到单元弹性力对节点坐标求导的雅可比矩阵$\textbf{JQ}$, 故而进一步给出其表达式

\begin{equation*}\begin{split}{JQ}_{i,j}=&\dfrac{\partial{Q}_{\rm e}}{\partial e_j}=\dfrac{\partial({Q}_{\rm s}+{Q}_\text{κ})} {\partial e_j}=\\&\int_V\left\lgroup\dfrac{\partial^2\textbf{ε}_v^{\rm T}}{\partial e_i\partial e_j}{E}_v\textbf{ε}_v+\dfrac{\partial\textbf{ε}_v^{\rm T}}{\partial e_i} {E}_v\dfrac{\partial\textbf{ε}_v}{\partial e_j}\right\rgroup{\rm d}V+\\&\int_V\left\lgroup\dfrac{\partial^2\kappa_v^{\rm T}}{\partial e_i\partial e_j}{E}\kappa_v+\dfrac{\partial \kappa_v^{\rm T}}{\partial e_i}{E}\dfrac{\partial \kappa_v}{\partial e_j}\right\rgroup{\rm d}V\end{split}\tag*{(15)}\end{equation*}

关于格林-拉格朗日应变张量$\textbf{ε}_v$与弯曲应变${\kappa}_v$对于节点坐标的二阶导数项推导可参考文献^[32], 本文给出梯度向量矩阵增加的单位法向量$ n$对于节点坐标的二阶导数表达式

\begin{equation*}\begin{split}\dfrac{\partial^2{n}}{\partial e_i\partial e_j}=&\dfrac{\partial}{\partial e_j}\left\lbrack \dfrac{\partial{g}}{\partial e_i}\cdot\left({g}^{\rm T}\cdot{g}\right)^{-\tfrac{1}{2}} -{g}\cdot\left({g}^{\rm T}\cdot{g}\right)^{-\tfrac{3}{2}}\cdot\right.\\&\left.\dfrac{\partial{g}^{\rm T}}{\partial e_i}\cdot{g}\right\rbrack=\dfrac{\dfrac{\partial^2{g}}{\partial e_i\partial e_j}}{\left({g}^{\rm T}\cdot{g}\right)^{\tfrac{1}{2}}}+\dfrac{3{g}\dfrac{\partial{g}^{\rm T}}{\partial e_j}{g}\cdot\dfrac{\partial{g}^{\rm T}}{\partial e_i}\cdot{g}}{\left({g}^{\rm T}\cdot{g}\right)^{\tfrac{5}{2}}}-\\&\dfrac{\dfrac{\partial {g}}{\partial e_i}\cdot\dfrac{\partial {g}^{\rm T}}{\partial e_j}\cdot{g}}{\left({g}^{\rm T}\cdot{g }\right)^{\tfrac{3}{2}}}-\dfrac{\dfrac{\partial{g}}{\partial e_j}\cdot\dfrac{\partial{g}^{\rm T}}{\partial e_i}\cdot{g}}{\left({g}^{\rm T}\cdot{g }\right)^{\tfrac{3}{2}}}-\\&\dfrac{{g}\cdot\dfrac{\partial^2{g}^{\rm T}}{\partial e_i\partial e_j}\cdot{g}+{g}\cdot\dfrac{\partial{g}^{\rm T}}{\partial e_i}\cdot\dfrac{\partial {g}}{\partial e_j}}{\left({g}^{\rm T}\cdot{g }\right)^{\tfrac{3}{2}}}\end{split}\tag*{(16)}\end{equation*}

由此, 本文给出了一种既能利用薄板单元不存在截面闭锁的优势, 又能将钢板弹簧这一参考构型带有初始弯曲的变形体中的初应变消除的弹性力描述方法. 由于该方法通过引入第3个梯度向量, 将原本梯度不完备的薄板单元纳入了只有全参数单元才可以应用的一般连续介质力学方法, 故称其为完备化的薄板单元弹性力模型. 该方法不仅可以将簧片中的初始应变消除, 还可以通过建立更大曲率的簧片未变形构型, 在钢板弹簧中引入可以人为控制的预应力.

3 钢板弹簧建模

3.1 板间接触搜索

钢板弹簧是多个簧片由中心螺栓紧固而成的, 如图3所示. 簧片装配前都有各自的曲率, 越短的簧片弯曲程度越大. 装配后由于中心螺栓的作用, 所有簧片在初始状态都具有相同的曲率.

由于实际装配前各簧片自由状态的曲率与簧片长度负相关, 因此簧片之间作用力最大的区域位于相邻两簧片较短簧片端头棱线位置, 这部分区域的正压力最大, 也是摩擦作用最为显著的部位. 为了将板间接触问题简化, 本文假设接触发生在下片的两端上表面棱线和上片的下表面之间, 如图4所示. 为此, 在下片棱线上布置接触点, 并对上片的下表面进行接触搜索. 设$Pc$点为下片中面上的棱线点, 上片中面的距离最近投影点为$Qc$, 二者需满足垂直条件

\begin{equation*}\left.\begin{split}&\left({r}_{Qc}-{r}_{Pc}\right)^{\rm T}\cdot{r}_{xQc}=0\\&\left({r}_{Qc}-{r}_{Pc}\right)^{\rm T}\cdot{r}_{yQc}=0\end{split}\right\}\tag*{(17)}\end{equation*}

其中${r}_{xQc}$与${r}_{yQc}$分别代表上片下表面的实际接触点$Q$点在单元中心面映射点$Qc$处$x$和$y$两个梯度向量.

文章在采用Newton-Raphson法对式(17)求解的过程中, 需要不断更新$Qc$点位置, 需要首先确定$Qc$所处的单元. 而薄板单元内部物质点坐标取值范围均为[0,1], 处理这一类跨单元搜索问题较为困难. 为此, 在上片板簧中面建立了簧片局部坐标系, 并定义物质点坐标$\xi$为

\begin{equation*}\xi=\bar{\xi}+\hat{\xi},~\xi\in[0,n_l]\tag*{(18)}\end{equation*}

其中, $n_l$为簧片长度方向单元数. $\xi$由整数部分$\bar{\xi}$与小数部分$\hat{\xi}$组成, $\bar{\xi}$ 表示当前$Qc$点所在的单元编号, $\hat{\xi}$表示当前$Qc$点所在该单元内部的物质坐标, 如图4所示.

3.2 接触力和摩擦力计算与施加

在使用主从方法完成接触点物质坐标查找后, 是否施加接触力以及摩擦力需要进行接触判定. 接触点副搜索对于上下片单元中面进行, 接触判定需要由实际接触表面点的嵌入距离$T$计算判断^[20]

\begin{equation*}T=\left({r}_Q-{r}_P\right)^{\rm T}\cdot{N},~{N}=\dfrac{{r}_{Qc}-{r}_{Pc}}{\big| {r}_{Qc}-{r}_{Pc}\big|}\tag*{(19)}\end{equation*}

其中表面点$P$与$Q$的全局坐标为

\begin{equation*}{r}_P={r}_{Pc}+\dfrac{H}{2}\cdot{N}_{Pc},~{r}_Q={r}_{Qc}-\dfrac{H}{2}\cdot{N}_{Qc}\tag*{(20)}\end{equation*}

式中, $H$表示单元厚度, ${N}_{Pc}$与${N}_{Qc}$分别为单元中面$Pc$与$Qc$点处的单位法向量, 如图5所示.

如果$T<0$, 则认为上下片单元在该点处发生了接触, 该点嵌入深度$\delta$为

\begin{equation*}\delta=\big|T\big|\tag*{(21)}\end{equation*}

嵌入速度的计算不仅涉及到法向接触力的施加, 也进一步参与摩擦力计算. 首先需要计算接触点副映射的表面点$Q$相对于$P$的速度项, 不同于已有文献中采用针对计算步长的相对位置差分^[28-29], 文章选择实际接触面上的点在各迭代步得到的速度实际值计算相对速度$V_{\rm relative}$

\begin{equation*}\begin{split}V_{\rm relative}=&\dot{ r}_Q-\dot{ r}_P=\\&\left\lgroup\dot{ r}_{Qc}-\dfrac{H}{2}\cdot\dfrac{\partial{N}_{Qc}}{\partial t}\right\rgroup-\left\lgroup\dot{ r}_{Pc}+\dfrac{H}{2}\cdot\dfrac{\partial{N}_{Pc}}{\partial t}\right\rgroup\quad\end{split}\tag*{(22)}\end{equation*}

将$V_{\rm relative}$在接触点副$Pc$与$Qc$点连线${N}$上投影得到嵌入速度$D_\delta$

\begin{equation*}D_{\delta}=V_{\rm relative}\cdot N\tag*{(23)}\end{equation*}

板间法向接触力大小$F_{Cn}$由接触刚度$K$与接触阻尼系数$C$采用惩罚函数法计算得到^[33]

\begin{equation*}F_{Cn}=K\cdot\delta-C\cdot D_{\delta}\cdot\delta\tag*{(24)}\end{equation*}

摩擦力的计算需要表面点$Q$相对于$P$在垂直于$Pc$与$Qc$点连线的平面内的滑移速度${V}_{\rm tan}$

\begin{equation*}{V}_{\rm tan}={V}_{\rm relative}-D_\delta\cdot {N}\tag*{(25)}\end{equation*}

对切向摩擦力$F_{Ct}$, 由于接触时存在的脱离、再接触、振动和爬行, 会导致直接使用摩擦系数$\mu$计算得到方向和大小都很不稳定的动摩擦力, 对于求解步长影响很大, 本文摩擦力计算采用计及临界滑动速度$V_t$ 的平滑化库伦摩擦计算方法

从而得到下片点$P$对上片点$Q$的接触力和摩擦力的合力矢量${F}_C$

\begin{equation*}{F}_C=F_{Cn}\cdot{N}_{Qc}+F_{Ct}\cdot\dfrac{{V}_{\rm tan}}{\big|\big|{V}_{\rm tan}\big|\big|}\tag*{(27)}\end{equation*}

最后, 给出该接触点副所对应的广义力${Q}_c$

\begin{equation*}{Q}_c=\dfrac{\partial{e}}{\partial{e}^i}\cdot{S}^{i^{\rm T}}\cdot {F}_C+\dfrac{\partial{e}}{\partial{e}^j}\cdot{S}^{j^{\rm T}}\cdot\left(-{F}_C\right)\tag*{(28)}\end{equation*}

并将${Q}_c$作为系统外力施加于系统动力学方程中. 式(28)中${e}^i$和${e}^j$分别对应为点$Q$和点$P$所在单元的节点坐标.

3.3 车身刚柔耦合建模}

在对钢板弹簧减震系统进行动力学建模仿真时, 为了较真实地反映系统的动态响应, 需要将与钢板弹簧连接的车身、吊耳等构件及其机构运动关系在模型中加以描述. 文章借助绝对节点坐标法中参考节点的概念, 在模型中刻画了吊耳、车身等刚体构件, 实现了机构运动关系的描述, 尝试了参考节点与梯度不完备单元混合建模, 获得了可以实现与整车模型整合进行该类减震系统刚柔耦合动力学仿真的建模方法.

文章选择了最为常见的钢板弹簧减震系统构成形式进行模拟, 如图6所示. 钢板弹簧两端均有铰接卷耳, 一端卷耳直接与车架通过柱销铰接, 另外一端先与吊耳铰接, 再通过吊耳铰接在车架上. 钢板弹簧的簧片中部通过U型扣与车轴相连, 钢板弹簧所负载的车身重量通过卷耳与吊耳施加给钢板弹簧系统. 通过设置实际的刚体质量与各轴惯量特性, 能较真实地再现钢板弹簧系统的动力学行为.

参考节点(reference nodes for ANCF)是在绝对节点坐标方法中引入用于描述与多柔体系统有机构连接关系的刚体的一种方法. 在使用绝对节点坐标方法建立的单元节点网格中, 参考节点是单独存在的, 不属于任何一个单元. 参考节点具有完备的梯度向量, 能够通过线性约束将参考节点的某些自由度与单元网格关联起来, 并通过非线性约束满足自身的刚体条件. 同时, 参考节点也能够通过释放特定的约束实现对运动副的描述. 参考节点已经成功地应用于全参数单元网格的建立与求解^[12], 本文将实现参考节点与薄板单元这种梯度不完备单元的混合建模.

将参考节点代表的刚体与代表柔体系统的单元网格节点关联需要使用线性约束, 将普通的单元网格节点的自由度相对于参考节点完全约束的形式为^[34]

\begin{equation*}\left.\begin{split}&{r}^{ik}={r}^{ir}+{J}^{ir}{r}_0^{ik}\\&{J}^{ik}={J}^{ir}{J}_0^{ikr}\end{split}\right\}\tag*{(29)}\end{equation*}

其中, ${r}^{ik}$和${r}^{ir}$分别代表一个普通单元网格节点$k$与一个ANCF参考节点的位置向量. ${J}^{ik}$和${J}^{ir}$分别表示普通节点$k$与参考节点的梯度向量矩阵. ${r}_0^{ik}$表示在初始时刻普通节点$k$相对于参考节点的位置向量, ${J}_0^{ikr}$表示初始时刻普通节点$k$相对于参考节点的梯度向量矩阵. 参照刚体位移的概念, 这里的位置向量${r}^{ir}$表示参考节点所代表的刚体的平动, 而梯度向量矩阵${J}^{ir}$则描述了刚体的转动, 如图7 所示. 式(29) 包含1 个位置向量与3个梯度向量共4 个向量方程, 相当于12 个标量方程, 从而将普通节点$k$ 的自由度全部限制. 此外, 只需将式(29)中关于转动的约束释放, 即可实现对铰接关系的描述. 值得注意的是, 式(29)中约束方程对于节点坐标的雅可比矩阵为常值矩阵, 故其为线性约束, 可以在动力学方程的程序前处理部分一次性施加, 消去非独立坐标, 减少系统的求解自由度个数.

ANCF中引入参考节点概念的目的是为了使用一个节点描述一个刚体, 为此必须对其梯度向量矩阵施加约束, 使之始终为正交阵, 约束方程为

\begin{equation*}\left.\begin{split}&\big|{r}_x^{ir}\big|=1,~\big|{r}_y^{ir}\big|=1,~\big|{r}_z^{ir}\big|=1\\&{r}_x^{ir{\rm T}}\cdot{r}_y^{ir}=0,~{r}_x^{ir{\rm T}}\cdot{r}_z^{ir}=0,~ {r}_y^{ir{\rm T}} \cdot{r}_z^{ir}=0\end{split}\right\}\tag*{(30)}\end{equation*}

可以看出, 式(30)中的约束方程是系统节点坐标的非线性函数, 其对节点坐标的雅可比矩阵时变, 需要在系统动力学方程中引入拉格朗日乘子进行求解.

由于参考节点不属于任何一个单元, 而系统的质量阵是对单元的积分计算得到, 因此需要给出参考节点的质量阵计算方法. 由式(29)可以看出, 参考节点采用了多刚体动力学描述中的自然坐标法描述其运动, 即刚体上任意点的坐标${r}={r}^r+{J}^r\bar{{u}}$, 其中${r}^r$相当于是参考节点的全局坐标系下位置坐标, ${J}^r=[{r}_x^r{r}_y^r{r}_z^r]$是参考节点的梯度向量矩阵, 且在式(30) 的作用下保证为正交阵. $\bar{{u}}=[x,y,z]^{\rm T}$ 则表示了该任意点在参考节点局部坐标系下的相对位置. 由此, 该点的位置向量可以写做用参考节点坐标描述的形式

\begin{equation*}{r}=[{I}, x{I}, y{I}, z{I}]\cdot\left\lbrack{r}^{r\rm T}, {r}_x^{r\rm T}, {r}_y^{r\rm T}, {r}_z^{r\rm T}\right\rbrack^{\rm T}\tag*{(31)}\end{equation*}

则在该点处的质量微元$\rho{\rm d}V$的惯性力虚功为: $\text{δ} W=(\rho {\rm d}V)\ddot{{r}}^{\rm T}\delta{r}$. 对该点的位置向量对时间求取二阶导数得到的加速度表示为$\ddot{{r}}=[{I}, x{I}, y{I},z{I}]^{\rm T}\cdot\left\lbrack \ddot{{r}}^{r\rm T}, \ddot{{r}}_x^{r\rm T}, \ddot{r}_y^{r\rm T}, \ddot{{r}}_z^{r\rm T}\right\rbrack^{\rm T}$,

进而可以得到参考节点所表示的刚体的质量阵为

\begin{equation*}{M}=\int_V\rho[{I}, x{I}, y{I}, z{I}]^{\rm T}\cdot[{I}, x{I}, y{I}, z{I}]{\rm d}V\tag*{(32)}\end{equation*}

由此可见, 参考节点的质量阵也是常值矩阵, 通过对绝对节点坐标单元网格质量阵进行增广即可实现系统组集.

4 系统动力学方程建立与求解

多体系统动力学方程的通用形式为

\begin{equation*}\left.\begin{split}&{M}\ddot{{q}}={Q}\\&{C}({q})={0}\end{split}\right\}\tag*{(33)}\end{equation*}

其中广义力项${Q}$包括系统弹性力${Q}_F$和广义重力${G}$以及接触力${Q}_c$等系统外力, ${M}$为系统质量阵, ${C}({q})$为系统的非线性约束方程.

前文介绍可知, 参考节点的引入同时为系统带来了线性与非线性约束方程. 针对线性约束部分, 将约束方程对于时间求导, 得到$C_q\cdot\dot{{q}}=0$, 其中${C}_q$表示约束方程对于节点坐标的导数项. 将节点坐标按照独立坐标${q}_i$与非独立坐标${q}_d$分离后, 代入得到约束方程的微分形式: ${C}_qd\cdot\text{δ}{q}_d+{C}_qi\cdot\text{δ}{q}_i=0$. 进一步得到非独立坐标对于时间的一阶导数项的表达形式: $\dot{{q}}_d=-{C}_q^{-1}d\cdot{C}_q{i}\cdot\dot{{q}}_i$, 从而表示出全部节点坐标对于时间的一阶导数

\begin{align*}\dot{{q}}=\begin{bmatrix}\dot{{q}}_d\\\dot{{q}}_i\end{bmatrix}=\begin{bmatrix}{C}_d{i}\cdot\dot{{q}}_i\\\dot{{q}}_i\end{bmatrix}=\begin{bmatrix}{C}_d{i}\\{I}\end{bmatrix}\dot{{q}}_i={B}_d{i}\cdot\dot{{q}}_i\tag*{(34)}\end{align*}

进一步对于动力学方程中的节点坐标的二阶导数进行自由度缩减

\begin{equation*}\ddot{{q}}={B}_d{i}\cdot\ddot{{q}}_i\tag*{(35)}\end{equation*}

代入系统动力学方程中进行自由度缩减得到

\begin{equation*}{M}_{ii}\ddot{{q}}_i={Q}_{ii}\tag*{(36)}\end{equation*}

其中${M}_{ii}$与${Q}_{ii}$分别为自由度缩减后的质量阵和广义力^[35]

\begin{equation*}\left.\begin{split}&{B}_d^{\rm T}{i}\cdot {M}\cdot {B}_d{i}={M}_{ii}\\&{B}_d^{\rm T}{i}\cdot {Q}={Q}_{ii}\end{split}\right\}\tag*{(37)}\end{equation*}

可以看出, 线性约束由于具有定常的雅可比矩阵, 可以对其进行分解, 区分出系统的独立与非独立坐标, 定义二者间的线性变换关系, 从而将系统动力学方程改写成只由独立坐标描述的形式, 实现了系统自由度的缩减. 而非线性约束则需要引入拉格朗日乘子, 建立起如式(38)所示的微分‒代数方程组进行求解.

\begin{equation*}\left.\begin{split}&{M}\ddot{{q}}-{Q}+{C}_q^{\rm T}\cdot{\lambda}_k={0}\\&{C}({q})={0}\end{split}\right\}\tag*{(38)}\end{equation*}

本文拟采用基于隐式积分策略的广义$\alpha$法对系统动力学方程进行求解. 广义$\alpha$法可以在耗散系统高频响应的同时较好地保持系统低频特性^[36], 相对较大的时间步长也使得其可以在计算效率与精度之间保持较好的均衡.

5 数值算例

5.1 完备化的薄板单元弹性力验证

本文针对描述带有初始弯曲参考构型的柔性体变形问题, 对原有薄板单元弹性力模型进行了改进, 提出了完备化的薄板单元弹性力模型. 而在描述初始直构型的柔性体变形时, 两种方法应当得到相同的结果. 文章选择现有绝对节点坐标薄板单元文献^[13,37]中常用的$n\times n$单元网格柔性板算例. 表1所示为柔性板水平下摆算例模型参数. 如图8 所示单角点铰接的正方形板水平位置静止释放后在重力作用下下落, 比较其铰接点对侧角点在仿真时间$0\sim0.3$ s 内竖直位移曲线.

由图9可以看出, 完备化的薄板单元在采用无数值阻尼求解得到的系统总能量维持不变. 同时, 由图10可以看出, 对比现有薄板单元与完备化的薄板单元在相同仿真参数下的铰接点对侧角点竖直位移曲线, 两者结果在描述初始直构型柔性体时完全一致, 符合预测.

5.2 含有双片钢板弹簧的车架刚柔耦合模型

使用参考节点描述车身与吊耳, 建立刚柔耦合的双片钢板弹簧模型, 如图11所示. 将上片与下片各自的一对中间位置节点的位置坐标${r}$与梯度向量${r}_X$固定, 模拟车轴在簧片中部对于车辆的支撑. 簧片前端棱线铰接在车身上, 后端棱线通过吊耳铰接在车身上. 其中参考节点1代表车身, 质量取\mbox{500 kg}; 参考节点2 与参考节点3 构成钢板弹簧与车身之间连接的吊耳. 簧片厚度均为0.014 m, 表2中板间接触参数参照文献^[38] 选取.

为了描述实际钢板弹簧系统在未受载时由于簧片装配已经引入的预应力, 文章针对弹性力计算定义的未变形构型区别于初始构型. 这样就在初始弯曲的参考构型中加入了预应力, 能够较大程度上增加钢板弹簧的整体刚度并提升阻尼.

簧片在自重与车身参考节点1以及吊耳参考节点2, 3的重力作用下变形, 上片下压与下片发生接触, 模拟出带有吊耳的钢板弹簧减震系统在整车重力作用下的响应. 仿真运行时间为$0\sim5$ s. 图12中给出车身参考节点在仿真过程中的竖直方向位移.

由图12看出, 由于系统中摩擦阻尼的存在, 系统动力学响应体现出了预期的阻尼特性. 相比于无预应力的仿真结果, 车身对应的参考节点竖直位移的波动幅值更小, 预应力的引入提高了钢板弹簧的整体刚度. 系统初始构型及车身最大竖直位移时刻构型如图13所示, 由$\kappa_{xx}$对应的弯曲应变云图可以看出下片初始时刻存在的预应变以及参考节点大幅下压后簧片内弯曲应变的分布.

本文直接采用节点速度参与接触力与摩擦力计算, 较为细腻地刻画了板间接触摩擦过程. 在只对簧片使用端部外力平滑加载时, 如图14所示, 整个系统产生高频响应, 局部放大后可以看出簧片在基频振动上叠加的高频振动. 这导致簧片接触后嵌入深度$\delta$方向上存在高频激振, 加之高频振动带来的嵌入速度$D_\delta$的高频换向波动, 由式(24)得出, 法向接触力的大小$F_{Cn}$将产生与之对应的较大幅度高频波动. 簧片表面上的接触点副在中面上的投影点之间相对速度$\dot{{r}}_{Qc}-\dot{{r}}_{Pc}$的大小波动与高频换向, 叠加上下片各自中面法向量在空间指向上出现的异步扭振, 由式(22)可以看出, 将引起摩擦力计算依赖的表面点相对速度${V}_{\rm relative}$高频换向波动的幅度进一步放大. 由此引起的簧片间摩擦力高频跃变会使得系统动力学仿真效率降低. 使用平滑化的库伦摩擦模型则可有效改善求解状况. 由于每次求解迭代均重新计算当前接触力, 各时间步计算收敛过程中的数值波动以及结构自身受外载冲击的高频响应、接触和摩擦对于速度项也存在影响, 故本文方法可较多计入接触与摩擦的变化对系统响应的影响. 模型仿真过程中, 求解器对于临界滑移速度的很大敏感性以及对于接触力向量的实时输出观察证实了相关猜测, 引入预应力增大系统刚度后求解时收敛步数的进一步增大也符合预测.

6 结论

为避免绝对节点坐标全参数单元在钢板弹簧建模中存在的截面闭锁问题, 本文建立了基于薄板单元的车辆钢板弹簧动力学模型. 为描述板簧的初始弯曲以及预应力, 本文在现有绝对节点坐标薄板单元基本位移假定的基础, 提出了完备化的薄板单元弹性力模型. 用对比算例证明了所提方法的正确性. 文章建立了带有初始弯曲及预应力的双片钢板弹簧模型. 提出了建立簧片局部坐标系进行接触点跨单元搜索的方法, 根据假设提出了计及表面速度的线‒面接触模型, 并进一步使用惩罚函数法和平滑化的库伦摩擦模型给出了接触力、摩擦力计算和加载方法. 文章利用参考节点进行了整合车身、吊耳与钢板弹簧的系统刚柔耦合动力学仿真. 结果表明, 本文提出的绝对节点坐标薄板单元的完备化弹性力模型以及钢板弹簧板间接触模型较好地捕捉到了钢板弹簧复杂的非线性动力学行为. 数值仿真结果同时也证明了绝对节点坐标参考节点与梯度不完备单元进行混合建模的可行性. 后续的工作中, 将继续研究钢板弹簧模型中预应力对系统动力学特性的影响以及更加复杂精细的面接触和摩擦描述方法.

The authors have declared that no competing interests exist.

---

参考文献

文献选项

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

[1]	Shabana AA. An absolute nodal coordinates formulation for the large rotation and deformation analysis of flexible bodies . Technical Report. No. MBS96-1-UIC, University of Illinois at Chicago, 1996 [本文引用: 1]
[2]	章孝顺, 章定国, 陈思佳等. 基于绝对节点坐标法的大变形柔性梁几种动力学模型研究 . 物理学报, 2016, 65(9): 094501 [本文引用: 1]      (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)) [本文引用: 1]
[3]	Gerstmayr J, Sugiyama H, Mikkola A. Review on the absolute nodal coordinate formulation for large deformation analysis of multi-body systems . Journal of Computational and Nonlinear Dynamics, 2013, 8(3): 031016 [本文引用: 1]
[4]	田强, 张云清, 陈立平等. 柔性多体系统动力学绝对节点坐标方法研究进展 . 力学进展, 2010, 40(2): 189-202 [本文引用: 1]      (Tian Qiang, Zhang Yunqing, Chen Liping, et al. Review of the advances in absolute nodal coordinate formulation in flexible multi-body system dynamics . Advances in Mechanics, 2010, 40(2): 189-202 (in Chinese)) [本文引用: 1]
[5]	张越, 赵阳, 谭春林等. ANCF 索梁单元应变耦合问题与模型解耦 . 力学学报, 2016, 48(6): 1406-1415 (Zhang Yue, Zhao Yang, Tan Chunlin, et al. The strain coupling problem and model decoupling of ANCF cable/beam element . Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(6): 1406-1415 (in Chinese))
[6]	岳宝增, 于嘉瑞, 吴文军. 多储液腔航天器刚液耦合动力学与复合控制 . 力学学报, 2017, 49(2): 390-396 (Yue Baozeng, Yu Jiarui, Wu Wenjun. Rigid and liquid coupling dynamics and hybrid control of spacecraft with multiple propellant tanks . Chinese Journal of Theoretical and Applied Mechanics, 2017, 49(2): 390-396 (in Chinese))
[7]	胡景晨, 王天舒. 一种 $O(n)$ 算法复杂度的递推绝对节点坐标法研究 . 力学学报, 2016, 48(5): 1172-1183 [本文引用: 1]      (Hu Jingchen, Wang Tianshu. A recursive absolute nodal coordinate formulation with $O(n)$ algorithm complexity . Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5): 1172-1183 (in Chinese)) [本文引用: 1]
[8]	刘铖, 田强, 胡海岩. 基于绝对节点坐标的多柔体系统动力学高效计算方法 . 力学学报, 2010, 42(6): 1197-1205 [本文引用: 1]      (Liu Cheng, Tian Qiang, Hu Haiyan. Efficient computational method for dynamics of flexible multibody systems based on absolute nodal coordinate . Chinese Journal of Theoretical and Applied Mechanics, 2010, 42(6): 1197-1205 (in Chinese)) [本文引用: 1]
[9]	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 Dynamic, 2004, 35: 313-329 [本文引用: 1]
[10]	Shabana AA. ANCF tire assembly model for multibody system applications . Journal of Computational and Nonlinear Dynamics, 2015, 10(2): 024504 [本文引用: 1]
[11]	Omar MA, Shabana AA, Mikkola A, et al. Multibody system modeling of leaf springs . Journal of Vibration and Control, 2004, 10(11): 1601-1638 [本文引用: 2]
[12]	Yu ZQ, Liu YG, Tinsley B, et al. Integration of geometry and analysis for vehicle system applications: Continuum-based leaf spring and tire assembly . Journal of Computational and Nonlinear Dynamics, 2015 11(3): 031011 [本文引用: 4]
[13]	Dufva K, Shabana AA. Analysis of thin plate structures using the absolute nodal coordinate formulation . Proceedings of the Institution of Mechanical Engineers Part K: Journal of Multibody Dynamics, 2005, 219(4): 345-355 [本文引用: 4]
[14]	Wang JS, Li ZK, Jiang QB. The analysis of composite leaf spring by finite element method and experimental measurements //FISITA 2012 World Automotive Congress Proceedings-Volume 7: Vehicle Design and Testing, Beijing, China, 2012-11-27, 803 [本文引用: 1]
[15]	Zhou ZF, Guo WF, Shen TJ, et al.Research and application on dynamic stiffness of leaf spring. FISITA 2012 World Automotive Congress Proceedings-Volume 10: Chassis Systems and Integration Technology Beijing, China, 2012-11-27, 2012, 111 [本文引用: 1]
[16]	Nunney MJ. Light and Heavy Vehicle Technology . Oxford: Society for Neuroscience, 2006: 352 [本文引用: 1]
[17]	Dunaevskii B. Calculation of solid beams and flat springs (leaf springs) of varying cross-section . Soviet Engineering Research, 1981, 61(4): 36-38 [本文引用: 1]
[18]	周孔亢, 陆建辉, 侯永涛等. 基于RecurDyn的钢板弹簧动力学模型的建立与参数辨识 . 机械工程学报, 2014, 50(4): 128-134 [本文引用: 1]      (Zhou Kongkang, Lu Jianhui, Hou Yongtao, et al. Dynamics modeling and parameter identification of leaf spring based on RecurDyn . Journal of Mechanical Engineering, 2014, 50(4): 128-134 (in Chinese)) [本文引用: 1]
[19]	段亮, 杨树凯, 宋传学等. 平衡悬架钢板弹簧动态特性的研究 . 机械工程学报, 2016, 52(6): 153-158 [本文引用: 1]      (Duan Liang, Yang Shukai, Song Chuanxue, et al. Research on dynamic characteristics of the tandem suspension leaf spring . Journal of Mechanical Engineering, 2016, 52(6): 153-158 (in Chinese)) [本文引用: 1]
[20]	Wriggers P, Zavarise G. On contact between three-dimensional beams undergoing large deflections . Communications in Numerical Methods in Engineering, 1997, 13(6): 429-438 [本文引用: 2]
[21]	Zavarise G, Wriggers P. Contact with friction between beams in 3-D space . International Journal for Numerical Methods in Engineering, 2000, 49: 977-1006 [本文引用: 1]
[22]	Litewka P, Wriggers P. Contact between 3D beams with rectangular cross-sections . International Journal for Numerical Methods in Engineering, 2002, 53: 2019-2041 [本文引用: 1]
[23]	Litewka P. Hermite polynomial smoothing in beam-tobema frictional contact . Computational Mechanics, 2007, 40: 815-826 [本文引用: 1]
[24]	Litewka P. Smooth frictional contact between beams in 3D//Nack-enhorst U, Wriggers P. (eds.) IUTAM Symposium on Computational Contact Mechanics. IUTAM Bookseries . Dordrecht: Springer, 2007, 157-176 [本文引用: 1]
[25]	Konyukhov A, Schweizerhof K. Geometrically exact covariant approach for contact between curves . Computer Methods in Applied Mechanics and Engineering, 2010, 199: 2510-2531 [本文引用: 1]
[26]	Litewka P. Multiple-point beam-to-beam contact finite element //19th International Conference on Computer Methods in Mechanics, Warsaw, Poland, 9-12, May, 2011 [本文引用: 1]
[27]	Durville D. Contact-friction modeling within elastic beam assemblies: an application to knot tightening . Computational Mechanics, 2012, 49(6): 687-707 [本文引用: 1]
[28]	Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional contacts of thin beams during large overall motions via absolute nodal coordinate formulation . Nonlinear Dynamics, 2014, 77(4): 1411-1425 [本文引用: 2]
[29]	Wang QT, Tian Q, Hu HY. Dynamic simulation of frictional multi-zone contacts of thin beams . Nonlinear Dynamics, 2016, 83(4): 1919-1937 [本文引用: 2]
[30]	Zahavi E. Analysis of a contact problem in leaf springs . Mechanics Research Communications, 1992, 19(1): 21-27 [本文引用: 1]
[31]	Shabana AA. Computational Continuum Mechanics.London: Cambridge University Press. 2012: 46-93 [本文引用: 1]
[32]	Liu C, Tian Q, Hu HY. New spatial curved beam and cylindrical shell elements of gradient-deficient absolute nodal coordinate formulation . Nonlinear Dynamics, 2012, 70(3): 1903-1918 [本文引用: 1]
[33]	Gantoi FM, Brown MA, Shabana AA. Finite element modeling of the contact geometry and deformation in biomechanics . Journal of Computational and Nonlinear Dynamics, 2013, 8(4): 041013 [本文引用: 1]
[34]	Shabana AA. ANCF reference node for multibody system analysis . Proceedings of the Institution of Mechanical Engineers Part K: Journal of Multibody Dynamics, 2015, 229(1): 109-112 [本文引用: 1]
[35]	Shabana AA.Dynamics of Multibody Systems. Cambridge: Cambridge University Press, 2005. [本文引用: 1]
[36]	田强. 基于绝对节点坐标方法的柔性多体系统动力学研究与应用 . [博士论文]. 武汉: 华中科技大学, 2009 [本文引用: 1]      (Tian Qiang, Flexible multibody dynamics research and application based on the absolute nodal coordinate method . [PhD Thesis]. Wuhan: Huazhong University of Science and Technology, 2009 (in Chinese)) [本文引用: 1]
[37]	Yu ZQ, Shabana AA. Mixed-coordinate ancf rectangular finite element . Journal of Computational and Nonlinear Dynamics, 2015, 10(6): 061003 [本文引用: 1]
[38]	Shabana AA, Zaazaa KE, Sugiyama H.Railroad Vehicle Dynamics: A Computational Approach. Boca Raton: CRC Press. 2007: 145-147 [本文引用: 1]