2015, Vol. 47
2. 上海交通大学工程力学系, 上海200240
柔性体的变形场离散问题是柔性多体系统动力学中的研究热点和难点之一. 合理的离散方法不仅能正确描述柔性体的变形场，还能有效节 省时间和资源. 目前运用较广的主要有假设模态法(assumed mode method,AMM)和有限元法(finite element method,FEM)[1, 2]. 假设模态法源于结构力学中梁的固有振型，其突出优点是能用较少的模态获得较好的逼近结果，计算效 率高[3]. 但当柔性体为不规则形状时，其振动振型如何表示本身就是难题. 有限元法是通用性很强的一种离散方法，并在工程分析 中得到广泛的应用. 然而由于单元或网格的存在，使其存在许多固有缺陷[4]. 有限元法存在的问题主要有：前处理繁琐，需要花费 大量时间划分网格、不易构造出高阶连续的场函数，通常只是位移连续而应力应变不连续；另外，复杂三维结构的网格生成和重分也是相 当困难和费时的[5]. 国内外已有学者开始寻找新的变形场离散方法[6, 7]. 最近，范纪华等[8, 9]将B样条插值和贝齐儿插 值方法应用于柔性体变形场的离散，并将该方法用于旋转柔性梁的动力学研究中.
近年来，无网格法(meshfree method)作为一种较新的离散方法得到了迅速发展，成为了科学和计算方法研究的热点[10]. 无网格 法在建立整个问题域的系统代数方程时，不需利用预定义的网格信息进行离散. 该方法只需节点信息，无需划分网格，克服了有限元法 前处理复杂的缺点. 在构造形函数的过程中，采用更多的节点插值，通常具有高阶连续性，从而提高了计算精度. 现有的无网格法有多种，如无网格点插值法(point interpolation method,PIM)[11, 12]、无网格径向基点插值法(radial point interpolation method,RPIM)[13]、无单元迦辽金法(element-free Galerkin method,EFG)[14, 15]、 无网格局部Petrov-Galerkin法 (meshless local Petrov-Galerkin method,MLPG)[16]、再生核粒子法(reproducing kernel particle method,RKPM) [17]等. 径向基点插值法是在点插值法的基础上，为避免采用多项式基的点插值法在计算插值函数时矩阵易于奇异的问题所提出 的[18]，并成功应用于一维和二维弹性力学的静力学及频率分析中，获得了与解析解吻合很好的结果. 目前，将无网格法应用于柔性多体系统动力学的研究鲜有报道.
本文采用径向基点插值法描述柔性梁变形，并在此基础上对柔性旋转悬臂梁进行动力学分析. 考虑柔性梁纵向拉伸变形和横 向弯曲变形，并计入由横向弯曲变形引起的纵向缩短，即非线性耦合项. 采用浮动坐标系描述系统运动，运用第二类拉格朗日方程建立系统的动力学方程，并编制了基于径向基点插值法的旋转柔性梁动 力学仿真软件，将仿真结果与假设模态法、有限元法[19]等传统离散法所得结果进行对比，表明无网格法应用于该领域的正确性.
1 旋转柔性梁动力学模型 1.1 系统的动能与势能图1为水平面内运动的中心刚体-柔性梁系统，中心刚体在平面内绕固定转轴旋转，其上以悬臂方式固结柔性梁，材料均匀且各向同性.
|
图 1 旋转柔性梁变形示意图 Fig.1 The deformation of rotating flexible beam |
以刚体转动中心$O$为原点建立惯性坐标系$Oij$,在柔性梁上建立浮动坐标系$O'i'j'$. 中心刚体的转动惯量和半径分别 为$J_{\rm oh} $和$a$；柔性梁的弹性模量、长度、密度、横截面积、截面惯性矩分别为$E$,$L$,$\rho$,$S$\linebreak 和$I$.
如图1所示，柔性梁上任意点$P$变形后的矢径在惯性坐标系$Oij$中可表示为
| $r = \Theta ({r_A} + {\rho _0} + u)$ | (1) |
式中
| ${ r}_{ A} = (a,0)^{\rm T}$ | (2) |
| ${\rho}_0 = (x,0)^{\rm T}$ | (3) |
| ${ u} = (u_x ,u_y )^{\rm T}$ | (4) |
| $ \Theta = \left( \begin{array}{l} \cos \theta {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - \sin \theta \\ \sin \theta {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cos \theta \end{array} \right) $ | (5) |
其中，$\Theta$为浮动坐标系相对于惯性坐标系的方向余弦矩阵. 变形矢量${ u}$在浮动基下的坐标为
| $ u = \left[\begin{array}{l} {u_x}\\ {u_y} \end{array} \right] = \left[\begin{array}{l} {w_1} + {w_{\rm{c}}}\\ {w_2} \end{array} \right]$ | (6) |
式中，$w_1 $为柔性梁轴向拉伸量，$w_2 $为柔性梁横向弯曲挠度，$w_{\rm c}$为柔性梁横向弯曲引起的纵向缩短量，即变形 位移的二次耦合项[20]，表达式为
| ${w_{\rm{c}}} = - \frac{1}{2}\int {_0^x} {\left( {\frac{{\partial {w_2}}}{{\partial \zeta }}} \right)^2}{\rm{ d}}\zeta {\rm{ }}$ | (7) |
传统的零次混合坐标建模方法中未考虑该项，当大范围运动为高速时，会出现动力柔化效应，梁末端横向变形不收敛，与实际不符.
对式(1)求时间的一阶导数，得到柔性梁上任意点$P$在惯性坐标系下的速度
| $\dot r = \dot \Theta ({r_A} + {\rho _0} + u) + \Theta \dot u$ | (8) |
系统动能由中心刚体动能和柔性梁的动能两部分组成，表示为
| $T = \frac{1}{2}{J_{{\rm{oh}}}}{\dot \theta ^2} + \frac{1}{2}\int {_{_V}} \rho {\dot r^{\rm{T}}}{\dot r^{\rm{T}}}\dot rdV$ | (9) |
系统变形势能表示为
| $U = \frac{1}{2}\int {_0^L} ES{\left( {\frac{{\partial {w_1}}}{{\partial x}}} \right)^2}{\rm{ d}}x + \frac{1}{2}\int {_0^L} EI{\left( {\frac{{{\partial ^2}{w_2}}}{{\partial {x^2}}}} \right)^2}{\rm{ d}}x$ | (10) |
在径向基点插值法中，计算域是用一系列点来离散的，某点处的位移是通过对该点支持域内其他节点位移进行插值而得到的. 设 定义在问题域$\Omega $中的一个连续函数$u({ x})$可用径向基和多项式基的线性组合表示. 在计算点${ x}$处$u({ x})$可近似表示为
| $ \begin{array}{l} u(x) = \sum\limits_{i = 1}^n {{R_i}(x){a_i}} + \sum\limits_{i = 1}^m {{p_i}(x){b_i}} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R^{\rm{T}}}(x)a + {p^{\rm{T}}}(x)b \end{array}$ | (11) |
式中，$R_i ({ x})$为径向基函数，$n$为计算点${ x}$处局部支持域中的节点数，$p_i ({ x})$为基函数在空间坐标${ x}^{\rm T} = [x,y]$上的给定单项式，$m$为多项式基函数的个数. 如$m = 0$，则为单纯的径向基函数，否则为添加了$m$个多项式基的径向基函数. $a_i $和$b_i $为待定系数.
在径向基函数$R_i ({ x})$中，仅有表示计算点${ x}$与节点${ x}_i $之间距离的变量，对于一维和二维空间，其变量分别表示如下
| $r = \sqrt {(x - x_i )^2} $ | (12) |
| $r = {\rm{ }}\sqrt {{{(x - {x_i})}^2} + {{(y - {y_i})}^2}} $ | (13) |
本文径向基采用复合2次函数[5]，表示为
| $R_i ({ x}) = [r_i^2 + (\alpha _c d_c )^2]^q$ | (14) |
式中，$r_i $为变量，$\alpha _c $和$q$为两个无量纲形状参数，具体数值由分析者确定，$d_c $为位于计算点附近的节点间距. 如节点是均匀分布的，$d_c $即可简单视为相邻两节点的距离；若节点分布不均匀时，可视为计算点支持域内的平均节点间距.
对于一维和二维空间，式(11)中的线性基函数${ p}^{\rm T}$分别为
| ${ p}^{\rm T}({ x}) = [1 \ \ x],\qquad m = 2$ | (15) |
| ${ p}^{\rm T}({ x}) = [1 \ \ x \ \ y],\qquad m = 3$ | (16) |
其二次基函数分别为
| ${ p}^{\rm T}({ x}) = [ 1 \ \ x \ \ {x^2}] ,\qquad m = 3 $ | (17) |
| ${ p}^{\rm T}({ x}) = [ 1 \ \ x \ \ y \ \ {x^2} \ \ {xy} \ \ {y^2}],\qquad m = 6$ | (18) |
完备的$p$阶多项式可由以下一般形式表示为
| ${ p}^{\rm T}({ x}) = [1 \ \ x \ \ {x^2} \ \ \cdots \ \ {x^{p - 1}} \ \ {x^p}]$ | (19) |
| $ { p}^{\rm T}({ x}) = [1 \ \ x \ \ y \ \ {x^2} \ \ {xy} \ \ {y^2} \ \ \cdots \ \ {x^p} \ \ {y^p}]$ | (20) |
为求解得到待定系数矩阵${ a}$和${ b}$，须形成计算点${ x}$处的支持域，如图2所示，其中包含$n$个场节点. 对问题 域中的任一计算点，其支持域尺寸$d_s $由下式确定
| $ d_s = \alpha _s d_c$ | (21) |
|
图 2 无网格法中计算点的支持域 Fig.2 The support domain of a point of interest |
式中，$\alpha _s $为该支持域的无量纲尺寸.
式(11)中待定系数矩阵${ a}$和${ b}$可通过令$u({ x})$等于该$n$个节点上的值来确定，即
| $ \left. {\begin{array}{*{20}{l}} \begin{array}{l} {u_1} = \sum\limits_{i = 1}^n {{R_i}({x_1}){a_i}} + \sum\limits_{i = 1}^m {{p_i}({x_1}){b_i}} \\ {u_2} = \sum\limits_{i = 1}^n {{R_i}({x_1}){a_i}} + \sum\limits_{i = 1}^m {{p_i}({x_1}){b_i}} \\ \vdots \\ {u_n} = \sum\limits_{i = 1}^n {{R_i}({x_n}){a_i}} + \sum\limits_{i = 1}^m {{p_i}({x_1}){b_i}} \end{array} \end{array}} \right\}$ | (22) |
上式可用如下矩阵表示为
| ${{\rm{U}}_{\rm{s}}} = {R_0}a + {P_m}b$ | (23) |
式中
| ${u_1 } \ \ {u_2 } \ \ {u_3 } \ \ \cdots \ \ {u_n } \}^{\rm T}$ | (24) |
是节点变形向量.
| $ {R_0}{\rm{ = }}{\left[\begin{array}{l} {R_1}({r_1}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_2}({r_1}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_n}({r_1})\\ {R_1}({r_2})\;{R_2}({r_2})\;\; \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_n}({r_2})\\ \vdots \;\;\;\;\;\;\;\;\;\;\;\; \vdots \;\;\;\;\;\;\;\; \ddots \;\;\;\;\; \vdots \;\\ {R_1}({r_n})\;{R_2}({r_n})\;\; \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {R_n}({r_n}) \end{array} \right]_{n \times n}}$ | (25) |
是径向基函数力矩矩阵.
| $ { a} = \left\{ {a_1 } \ \ {a_2 } \ \ {a_3 } \ \ \cdots \ \ {a_n } \right\}^{\rm T}$ | (26) |
是径向基未知系数向量.
| $ {p_m}{\rm{ = }}{\left[\begin{array}{l} {p_1}({x_1}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_2}({x_1}){\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_m}({x_1})\\ {p_1}({x_2})\;{p_2}({x_2})\;\; \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_m}({x_2})\\ \vdots \;\;\;\;\;\;\;\;\;\;\;\; \vdots \;\;\;\;\;\;\;\; \ddots \;\;\;\;\; \vdots \;\\ {p_1}({x_n})\;{p_2}({x_n})\;\; \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {p_m}({x_n}) \end{array} \right]_{n \times m}}$ | (27) |
是多项式力矩矩阵.
| $ { b} = \left\{ {b_1 } \ \ {b_2 } \ \ {b_3 } \ \ \cdots \ \ {b_m } \right\}^{\rm T}$ | (28) |
是多项式未知系数向量.
式(25)中，$R_i (r_k )$的$r_k $在一维与二维中的定义分别为
| $ r_k = \sqrt {(x_k - x_i )^2}$ | (29) |
| ${r_k} = {\rm{ }}\sqrt {{{({x_k} - {x_i})}^2} + {{({y_k} - {y_i})}^2}} $ | (30) |
由于式(23)中有$n + m$个变量，可使用下面$m$个约束条件添加$m$个方程
| $\sum\limits_{i = 1}^m {{p_j}({x_i}){a_i}} = P_m^{\rm{T}}a = 0,\;\;\;j = 1,2,\cdots ,m$ | (31) |
联立式(23)与式(31)可得如下矩阵方程
| $ {\tilde U_s} = \left[{\begin{array}{*{20}{c}} \begin{array}{l} {U_s}\\ {\bf{0}} \end{array} \end{array}} \right] = \left[\begin{array}{l} {R_{0}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {P_m}\\ P_m^T{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 \end{array} \right]\left[{\begin{array}{*{20}{c}} \begin{array}{l} a\\ b \end{array} \end{array}} \right] = G{a_0}$ | (32) |
求解式(32)可得未知系数矩阵
| ${a_0} = {G^{ - 1}}{\tilde U_s}$ | (33) |
将式(33)回代到式(11)中得
| $ u(x) = \left\{ {{R^{\rm{T}}}(x)\;\;{p^{\rm{T}}}(x)} \right\}{G^{ - 1}}{\tilde U_s} = {\rm{ }}{\tilde \Phi ^{\rm{T}}}(x){\tilde U_s}$ | (34) |
式中，$\tilde \Phi (x)$为径向基点插值法形函数向量，定义为
| ${\tilde \Phi ^{\rm{T}}}(x) = \left\{ {{\phi _1}(x)\;\;{\phi _2}(x)\;\; \cdots \;{\phi _{n + m}}(x)} \right\}{\rm{ }}$ | (35) |
最终对应于节点位移向量的径向基点插值法形函数可表示为
| ${\Phi ^{\rm{T}}}(x) = \left\{ {{\phi _1}(x)\;\;{\phi _2}(x)\;\; \cdots \;{\phi _n}(x)} \right\}$ | (36) |
以上是径向基点插值法形函数的形成过程. 该形函数具有克罗内克函数性质，可表示为
| $ {\phi _i}(x = {x_j}) = \left\{ \begin{array}{l} 1,i = j,i,j = 1,2 \cdots ,n\\ 0,i \ne j,i,j = 1,2 \cdots ,n \end{array} \right.$ | (37) |
因此可以用同有限元法一样的方式处理边界条件. 本文中柔性悬臂梁对应的边界条件为梁固定端纵向变形、横向变形及转角为0.
需要注意的是，径向基点插值法在积分过程中，不同积分点对应不同的支持域，即不同的积分点其形成的形函数矩阵可能不同，这与有限元法是不同的. 在有限元法中，一个单元内所有的积分点均采用相同的节点进行插值.
如图3所示，在梁轴上沿$x$方向将问题域用$N$个场节点表示，即
| $0 = x_1 ＜ x_2 ＜ \cdots ＜ x_N = L$ | (38) |
|
图 3 无网格法中柔性梁的离散 Fig.3 The discrete of the flexible beam |
取梁轴向和横向变形位移函数
| ${w_1} = \sum\limits_{i = 1}^n {{\Phi _{xi}}} (x){d_i}(t) = {\Phi _x}d\\ {w_2} = \sum\limits_{i = 1}^{2n} {{\Phi _{yi}}(x){d_i}(t) = {\Phi _y}d} $ | (39) |
式中，${\Phi}_{xi} (x)$和${\Phi}_{yi}(x)$分别为梁轴向和横向变形的形函数矩阵；${ d}_i(t)$为节点$i$轴向变形、横向变形和转角随时间变化量的列矢量. 分别为
| ${\Phi _{xi}}(x) = \left( {{\phi _{xi}}\;\;0\;\;0} \right){\rm{ }}\\ {\Phi _{yi}}(x) = \left( {0\;\;{\phi _{yi}}\;\;{\phi _{\theta i}}} \right){\rm{ }}\\ {d_i} = {\left( {{u_{xi}}\;\;{u_{yi}}} \right)^{\rm{T}}}$ | (40) |
式中，$u_{xi} $为第$i$个节点的轴向变形，${ u}_{yi} $为第$i$个节点的横向变形和转角组成的行阵.
将梁的变形位移函数代入式(6)中，得到梁上任意点的变形位移为
| $u = \left[ {\begin{array}{*{20}{l}} \begin{array}{l} {u_x}\\ {u_y} \end{array} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} \begin{array}{l} {\Phi _x}d - \frac{1}{2}{d^{\rm{T}}}Hd\\ {\Phi _y}d \end{array} \end{array}} \right]$ | (41) |
其中，${ H}(x)$为耦合形函数，表达式如下
| $ H(x) = {\rm{ }}\int {_0^x} {\Phi ^\prime }_y^{\rm{T}}(\zeta ){\Phi ^\prime }_y(\zeta )\Phi \zeta {\rm{ }}$ | (42) |
其中，${\Phi}'_y (\zeta )$表示${\Phi}_y^ (\zeta )$对$\zeta $求一阶导数.
对式(31)求时间$t$的一阶导数，可得到梁上任意点的变形速度为
| $\dot { u} = \left[\!\!\begin{array}{l} \dot {u}_x \\ \dot {u}_y \\ \end{array} \!\!\right] = \left[\!\!\begin{array}{l} {\Phi}_x \dot{ d} - { d}^{\rm T}{ H}\dot{ d} {\Phi}_y \dot{ d} \end{array} \!\!\right]$ | (43) |
将柔性梁轴向和横向变形位移函数代入系统的动能和变形势能表达式中，取广义坐标${ q}= (\theta ,{ d}^{\rm T})^{\rm T}$,运用第二类拉格朗日方程
| $\frac{{\rm{d}}}{{dt}}\left( {\frac{{\partial T}}{{\partial \dot q}}} \right) - \frac{{\partial T}}{{\partial q}} = - \frac{{\partial U}}{{\partial q}} + {F_q}$ | (44) |
式中，${ F}_q = (\tau ,{\bf 0})^{\rm T}$为广义力列阵，$\tau $为刚体上合外力关于转动中心$O$的主矩. 由式(7)可以看出横向弯曲引起的纵向缩短量$w_c $是横向变形$w_2 $的二阶小量，作适当的简化，舍去与$w_c $相关的一些高阶小量，如$w_c ^2$，$w_1 w_c $，$\dot {w}_2 w_c $等，可得系统的动力学方程为
| $\left( \begin{array}{l} {M_{11}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {M_{12}}{\kern 1pt} \\ {M_{21}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {M_{22}} \end{array} \right)\left( {\begin{array}{*{20}{c}} \begin{array}{l} {\ddot \theta }\\ {\ddot d} \end{array} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} \begin{array}{l} {Q_\theta }\\ {Q_d} \end{array} \end{array}} \right)$ | (45) |
式中
| $\begin{array}{l} {M_{11}} = {J_{{\rm{oh}}}} + {J_{{\rm{ob}}}} + 2{S_x}d + {d^{\rm{T}}}({M_1} + {M_2})d - \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \underline {{d^{\rm{T}}}(aC + D){d^{\rm{T}}}} \end{array}$ | (46) |
| ${ M}_{21} = { M}_{12} ^{\rm T} = ({ M}_3 ^{\rm T} - { M}_3 ){ d}-{ S}_y ^{\rm T}$ | (47) |
| $\begin{array}{l} {M_{22}} = {M_1} + {M_2} = \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \int {_0^L} \rho S\Phi _x^{\rm{T}}{\Phi _x}dx + \int {_0^L} \rho S\Phi _y^{\rm{T}}{\Phi _y}dx \end{array}$ | (48) |
| $ \begin{array}{l} {Q_\theta } = \tau - 2\dot \theta [{S_x}\dot d + {d^{\rm{T}}}({M_1} + {M_2})\dot d] + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \underline {2\dot \theta {d^{\rm{T}}}(aC + D)\dot d}] \end{array}$ | (49) |
| $ \begin{array}{l} {Q_d} = {{\dot \theta }^2}S_x^{\rm{T}} + 2\dot \theta ({M_3} - M_3^{\rm{T}})\dot d + \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} [{{\dot \theta }^2}({M_{22}}\underline { - aC - D} ) - K]d \end{array}$ | (50) |
以上各式中 ${ S}_x $,${ S}_y $,${ M}_1 $,${ M}_2 $,${ M}_3 $,${ C}$,${ D}$和${ K}$为常系数矩阵，带下划线的项是附加耦合项，由于在变形位移中考虑了横向变形引起的纵向缩短，这些附加的耦合项出现在广义质量阵和广义力阵中. 以上各式为定义在全局问题域上的，为计算上式中的积分需将整个问题域离散成一组不相互重叠的积分背景网格. 取柔性梁相邻两节点的间距为积分域，总体积分可表示为
| $\int {_\Omega } G{\rm{d}}\Omega = {\rm{ }}\sum\limits_k^{{n_c}} {\int {_{\Omega k}} } GD\Omega $ | (51) |
式中，$n_c $为积分背景网格个数，如图2所示，${ G}$表示被积函数，$\Omega _k $为第$k$个积分背景网格的域. 采用高斯积分法求解数值积分，在每个积分背景网格中使用$n_g $个高斯点，式(51)变为
| $\int {_\Omega } G{\rm{d}}\Omega = {\rm{ }}\sum\limits_k^{{n_c}} {\int {_{\Omega k}} } GD\Omega = \sum\limits_k^{{n_c}} {\sum\limits_{i = 1}^{{n_g}} {{{\mathop w\limits^\frown }_i}G({x_{Qi}})\left| {{J_{ik}}} \right|{\rm{ }}} } $ | (52) |
式中，${\mathop{w}\limits^{\frown}}_i $为第$i$个高斯积分点${ x}_{Qi} $的高斯加权因子，$J_{ik} $为对积分背景网格$k$在积分点${ x}_{Qi} $处进行积分的雅克比矩阵.
为得到各常数阵的数值解，将整个问题域中的场节点从1到$N$顺序编号，以方阵${ K}$为例，可得
| $\begin{array}{l} {K_{IJ}} = \sum\limits_k^{{n_c}} {\sum\limits_{i = 1}^{{n_g}} {\underbrace {ES{{\mathop w\limits^\frown }_i}\phi _{xI}^{{\rm{,T}}}({x_{Qi}})\phi _{xJ}^,({x_{Qi}})|{J_{ik}}|}_{K_{IJ}^{ik}} + } } \qquad \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_k^{{n_c}} {\sum\limits_{i = 1}^{{n_g}} {\underbrace {EI{{\mathop w\limits^\frown }_i}\phi _{yI}^{{\rm{,,T}}}({x_{Qi}})\phi _{yJ}^{,,}({x_{Qi}})|{J_{il}}|}_{K_{IJ}^{il}}} } \end{array}$ | (53) |
式(53)表示节点矩阵${ K}_{IJ} $通过局部支持域中包含节点$I$和节点$J$的所用积分点的贡献之和得到的数值结果. 如果节点$I$和节点$J$不同时位于积分点${ x}_{Qi} $的局部支持域中，则${ K}_{IJ}^{ik} $和${ K}_{IJ}^{il} $为0. ${ K}$的形式为
| $ K = \left[\begin{array}{l} {K_{11}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{12}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{1N}}\\ {K_{21}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{22}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{2N}}\\ {K_{N1}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{N2}}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \cdots {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {K_{NN}} \end{array} \right]$ | (54) |
类似的，以列阵${ S}_y $为例，可得
| $\begin{array}{l} {S_y} = \sum\limits_k^{{n_c}} {\sum\limits_{i = 1}^{{n_g}} {\underbrace {\rho S{{\mathop w\limits^\frown }_i}\phi _{xI}^{'{\rm{T}}}(a + {x_{Qi}})\phi _{yI}^{\rm{T}}({x_{Qi}})|{J_{ik}}|}_{s_I^{ik}}} } = \qquad \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \sum\limits_k^{{n_c}} {\sum\limits_{i = 1}^{{n_g}} {S_l^{ik}} } ,\;\;\;I = 1,2, \cdots ,N \end{array}$ | (55) |
按照以上方阵与列阵的求法，可得到式(45)中各矩阵的形式.
2 大范围运动已知的动力学仿真大范围运动已知时，式(45)转化为
| ${M_{22}}\ddot d = {Q_d} - {M_{21}}\ddot \theta $ | (56) |
为了验证无网格径向基点插值法的正确性及分析大范围运动变形位移，令中心刚体半径$a = 0$，梁的参数取与文献[21]中相同，分别为：长度$L = 10 $m，横截面积$S = 4\times 10^{ - 4}$m$^2$，截面 惯性矩$I = 2\times 10^{ - 7}$m$^4$，体积密度$\rho = 3000$kg/m$^3$，弹性模量$E = 70$GPa.
假定悬臂梁由静止开始作大范围旋转运动，大范围运动展开角速度规律为
| $ \dot \theta = \left\{ \begin{array}{l} \frac{{{\Omega _0}}}{T}t - \frac{{{\Omega _0}}}{{2\pi }}\sin \left( {\frac{{2\pi }}{T}t} \right),0 \le t \le T\\ {\Omega _0}{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} t{\rm{ > }}T \end{array} \right.$ | (57) |
式中，$T = 15$s，15s后转速到达$\Omega _0 $，然后以该速度匀速旋转. $\Omega _0 $分别取为$\Omega _0 = 4$rad/s和$\Omega _0 = 20$rad/s.
图4表示$\Omega _0 = 4$rad/s时柔性梁末端纵向变形，即$u_x $. 如图4所示，径向基点插值法、假设模态法和有限元法基本重合. $u_x $由两部分组成：柔性梁轴向拉伸量$w_1 $和柔性梁横向弯曲引起的纵向缩短量$w_c $. 末端纵向变形为负值，说明$w_1 $相对于$w_c $为小量，即纵向变形主要由横向弯曲导致的纵向缩短引起. 图5表示$\Omega _0 = 4$rad/s时柔性梁末端横向变形，即$u_y $. 如图5所示，径向基点插值法、假设模态法和有限元法的振幅及频率基本一致. 图6表示$\Omega_0 = 4$rad/s时柔 性梁末端横向变形速度，如图6所示，3种离散方法的仿真结果基本一致. 图7表示$\Omega_0 = 20$rad/s时柔性梁末端 纵向变形，为了体现在变形较大的情况下假设模态法的不适用性，将弹性模量减小到$E = 30$GPa. 从图中可明显看出径向基 点插值法和有限元法仿真结果基本一致，而假设模态法则出现了较大的误差. 图8、图9分别表示$\Omega _0 = 20$rad/s，$E = 30$GPa时柔性梁末端横向弯曲变形和变形速度，从图中可更明显的看出，径向基点插值法和有限元法的仿真结果基本一致，但假设模态法的仿真结果已经与另两种方法出现较大的差别. 15s后的局部放大图说明假设模态法的振动相位与振幅误差很大，进一步说明在变形较大的情况下，假设模态法的精度会降低. 图10表示$\Omega _0 = 4$rad/s，$E = 5$GPa时柔性梁末端横向弯曲变 形. 如图10所示，梁的最大变形超过了5m，属于大变形，径向基点插值法与有限元法结果仍然基本一致，而假设模态法已经发散，说明源 于结构力学中固有振型的假设模态法，适用范围局限于小变形情况，不能处理大变形问题，而本文方法和有限元法适用于大变形问题.
|
图 4 柔性梁末端纵向变形($\Omega _0 = 4$rad/s) Fig.4 Tip deformation of the beam along $x$ direction ($\Omega _0 = 4$rad/s) |
|
图 5 柔性梁末端横向变形($\Omega _0 = 4$rad/s) Fig.5 Tip deformation of the beam along $y$ direction ($\Omega _0 = 4$rad/s) |
|
图 6 柔性梁末端横向变形速度($\Omega _0 = 4$rad/s) Fig.6 Tip deformation rate of the beam along $y$ direction |
|
图 7 柔性梁末端纵向变形($\Omega _0 = 20$rad/s) Fig.7 Tip deformation of the beam along $x$ direction |
|
图 8 柔性梁末端横向变形($\Omega _0 = 20$rad/s) Fig.8 Tip deformation of the beam along $y$ direction |
|
图 9 柔性梁末端横向变形速度($\Omega _0 = 20$rad/s) Fig.9 Tip deformation rate of the beam along $y$ direction |
|
图 10 柔性梁末端横向变形($\Omega _0 = 4$rad/s, $E=5$GPa) Fig.10 Tip deformation of the beam along $y$ direction |
表1和表2分别表示$\Omega _0 = 4$rad/s和$\Omega _0 = 20$rad/s时3种离散方法的计算相对时间、相对误差及大范围旋转角 速度恒定时的响应振幅. 其中，假设模态法纵、横模态截断数分别各取3$\sim $7阶. 计算相对时间以3阶模态截断数为标准，相对误差以有限元法为标准. 有限元法将梁离散为10个单元，径向基点插值法将梁离散为11个节点. 从表中可以看出，模态截断数为3阶的假设模态法计算效率最高，但随着转动速度的增加，大范围旋转角速度恒定时的振幅与有限元法 对比误差越来越大，当$\Omega _0 = 20$rad/s时，误差达到48.42{\%}；同一转动角速度下，假设模态法的模态截断数越多，计算效率越低，但精度并没有明显提高； 径向基点插值法比有限元法计算效率稍高，无论低速转动还是高速转动，大范围恒定时的振幅两者基本一致.
|
表 1 3种离散方法的计算相对时间、相对误差、旋转角速度恒定时响应振幅 ($\Omega _0 = 4$rad/s) Table 1 The relative time, the relative error and the amplitude of these three discrete methods ($\Omega _0 = 4$rad/s) |
|
|
表 2 3种离散方法的计算相对时间、相对误差、旋转角速度恒定时响应振幅 ($\Omega _0 = 20$rad/s) Table 2 The relative time, the relative error and the amplitude of these three discrete methods ($\Omega _0 = 20$rad/s) |
本文采用常用的3种时间积分方法对动力学方程(56)进行求解，分别为4阶龙格库塔法、亚当姆斯预报校正法和纽马克法. 其中前两 种方法为显式方法，后一种为隐式方法，具体的迭代格式文献[22]中有详细阐述. 表3表示$\Omega_0 = 4$rad/s时，径向基点插值法和有限元法在上述3种时间积分方法下的计算时间. 从表中可以看出，3种方法的计算时间不在一个数量级上，纽马克法计算效率最高，最快达到0.47s. 4阶龙格库塔法次之，而亚当姆斯预报校正法最低. 这是由于时间步长引起的，隐式方法如纽马克法可以将时间步长取的较大，而另外两种显式方法的时间步长相比而言就要取的 小很多结果才能收敛. 对于显式方法，径向基点插值法比有限元计算效率略高，而隐式方法则比有限元低. 这是因为径向基点插值法计算形函数时更加复杂，纽马克法中时间主要消耗在形函数的计算中.
|
|
表 3 3种时间积分方法的计算时间 ($\Omega _0 = 4$rad/s) Table 3 The computation time of three time integration methods ($\Omega _0 = 4$rad/s) |
在求解径向基点插值法形函数时用到了两个无量纲形状参数$\alpha _c $和$q$，通常没有最佳的理论值来确定这两个参数. 文献[5]中给出结论：$q$为0.98和1.03时能取得较为精确的结果. 根据上述结 论，本文研究另一个形状参数$\alpha_c $对仿真结果的影响. 如图11，横坐标为形状参数$\alpha _c $的取值，纵坐标为旋转角速度恒定时响应振幅与有限元法的相对误差. 从图中可明显看出，$\alpha _c =1$时，误差最大，甚至超过了50{\%}，当$\alpha _c $取$2 \sim 7$时，误差就很小了，几乎与有限元法结果一致. 当$\alpha _c $取值再大时，将会导致径向基点插值法力矩矩阵奇异而计算失败. 因此$\alpha _c $取$2 \sim 7$可得到令人满意的结果.
|
图 11 径向基点插值法形状参数对仿真结果的影响 Fig.11 The influence of shape parameter of RPIM on results |
(1)模态截断数取3阶的假设模态法计算效率最高，模态截断数越多，计算效率越低；径向基点插值法比有限元法计算效率稍高. 随着 变形的增大，假设模态法精度降低，即使增加模态阶数，精度也不会明显提高，但计算效率会降低明显，而径向基点插值法与有限元法 仿真结果基本一致.
(2)在处理大变形问题的研究中，假设模态法结果发散，不能用于求解大变形问题. 径向基点插值法与有限元法结果收敛且基本一致，说 明这两种方法适用于求解大变形问题.
(3)求解动力学方程时，隐式方法纽马克方法效率最高，4阶龙格库塔法次之，而亚当姆斯预报校正法最低. 因此，对于柔性悬臂梁动力 学方程的求解优先选用隐式方法.
(4)径向基点插值法形状参数$\alpha _c =1$时，误差很大，当$\alpha _c $取值超过7时，将会导致径向基点插值法力矩矩阵奇异而计算失败. $\alpha _c $取$2 \sim 7$可得到令人满意的结果.
| [1] | 洪嘉振,尤超蓝. 刚柔耦合系统动力学研究进展. 动力学与控制学报,2004,2(2):1-6 (Hong Jiazhen, You Chaolan. Advances in dynamics of rigid-flexible coupling system. Journal of Dynamics and Control, 2004, 2(2):1-6 (in Chinese)) |
| [2] | Dwivedy SK, Eberhard P. Dynamic analysis of flexible manipulators, a literature review. Mechanism and Machine Theory , 2006, 41(7): 749-777 |
| [3] | 吴胜宝,章定国,康新. 刚体-微梁系统的动力学特性. 机械工程学报, 2010,46(3): 76-82 (Wu Shengbao, Zhang Dingguo, Kang Xin. Dynamic properties of hub-microbeam system. Journal of Mechanical Engineering , 2010, 46(3):76-82 (in Chinese)) |
| [4] | 张雄,刘岩,马上. 无网格法的理论及应用. 力学进展,2009, 39(1): 1-36 (Zhang Xiong, Liu Yan, Ma Shang. Meshfree methods and their applications. Advanes in Mechanics, 2009: 39(1): 1-36 (in Chinese)) |
| [5] | 刘桂荣,顾元通.无网格法理论及程序设计. 王建明,周学军 译. 济南:山东大学出版社,2007:45-133 (Liu Guirong, Gu Yuantong. An Introduction to Meshfree Methods and Their Programming. Shandong: Shandong University Press, 2007: 45-133 (in Chinese)) |
| [6] | Sanborn GG, Shanana AA. On the integration of computer aided design and analysis using the finite element absolute nodal coordinate formulation. Multibody Syst Dyn , 2009, 22: 181-197 |
| [7] | Sanborn GG, Shanana AA. A rational finite element method based on the absolute nodal coordinate formulation. Nonlinear Dyn , 2009, 58: 565-572 |
| [8] | 范纪华,章定国. 旋转悬臂梁动力学的B样条差值方法. 机械工程学报,2012, 48(23): 59-64 (Fan Jihua, Zhang Dingguo. B-spline interpolation method for the dynamics of rotating cantilever beam. Journal of Mechanical Engineering , 2012, 48(23): 59-64 (in Chinese)) |
| [9] | 范纪华,章定国. 旋转柔性悬臂梁动力学的 Bezier 插值离散方法研究. 物理学报,2014, 63(15): 154501 (Fan Jihua, Zhang Dingguo. Bezier interpolation method for the dynamics of rotating flexible cantilever beam. Acta Phys Sin , 2014, 63(15): 154501 (in Chinese)) |
| [10] | Belytschko T, Krongauz Y, Organ D, et al. Meshless methods: an overview and recent developments. Comput Methods Appl Mech Eng, 1996, 139: 3-47 |
| [11] | Liu GR, Gu Y, A point interpolation method for two-dimensional solids. Int J Muner Meth Engng, 2001, 50: 937-951 |
| [12] | Liu GR, Dai KY, Lim KM, et al. A point interpolation mesh free method for static and frequency analysis of two-dimensional piezoelectric structures. Computational Mechanics, 2002, 29(6): 510-519 |
| [13] | Wang JG, Liu GR. A point interpolation meshless method based on radial basis functions. Int J Numer Meth Eng, 2002, 54(11): 1623-1648 |
| [14] | Belytschko T, Lu YY, Gu L. Element free Galerkin methods. International Journal for Numerical Methods in Engineering , 1994, 37(2): 229-256 |
| [15] | Cheng RJ, Cheng YM, Ge HX. Element-free Galerkin (EFG) method for a kind of two-dimensional linear hyperbolic equation. Chin Phys B, 2009, 18(10): 4059-4064 |
| [16] | Atluri SN, Zhu T. A new meshless local Petrov-!-Galerkin (MLPG) approach in computational mechanics. Computational Mehcanics, 1998, 22(2): 117-127 |
| [17] | Liu WK, Jun S, Zhang YF. Reproducing kernel particle methods. Int J Numer Meth Fluids , 1995, 20: 1081-1106 |
| [18] | 夏茂辉,贾延,刘广君. 弹性力学问题的径向点插值法. 青岛科技大学学报, 2006, 27(2): 155-158 (Xia Maohui, Jia Yan, Liu Guangjun. A radial point interpolation method for the elasticity problem. Journal of Qing dao University of Science and Technology, 2006, 27(2): 155-158 (in Chinese)) |
| [19] | Chung J, Yoo HH. Dynamic analysis of a rotating cantilever beam by using the finite element method. Journal of Sound and Vibration, 2002, 249(1): 147-164 |
| [20] | 陈思佳,章定国,洪嘉振. 大变形旋转柔性梁的一种高次刚柔耦合动力学模型. 力学学报,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(2): 251-256 (in Chinese)) |
| [21] | 章定国,余纪邦. 作大范围运动的柔性梁的动力学分析. 振动工程学报, 2006,19(4): 475-479 (Zhang Dingguo, Yu Jibang. Dynamical analysis of a flexible cantilever beam with large overall motions. Journal of Vibration Engineering, 2006,19(4): 475-479 (in Chinese)) |
| [22] | 马振华. 现代应用数学手册. 北京:清华大学出版社,2005. 567-611 (Ma Zhenhua. The Application of Modern Mathematics Handbook. Beijing: Tsinghua University Press, 2005. 567-611 (in Chinese)) |
2. Dept. of Engineering Mechanics, Shanghai Jiaotong University, Shanghai 200240, China