0.   引言

柔性多体系统(flexible multibody system，FMS)已被广泛应用于航空航天工程、车辆工程、机器人工程等多个领域.然而，其动力学模型——微分代数方程组(differential algebraic equations，DAEs)在仿真中求解耗时巨大，效率问题尤为突出.为了提高求解效率，已经提出了许多模型降阶方法^[1-2]，可以分为两类：一类与振动模态无关，如Krylov子空间方法^[3]，适当的正交分解^[4-5]等；另一类是模态降阶，如模态综合法(component mode synthesis，CMS)^[6]，自适应模态集成^[7]和全局模态参数化(global modal parameterization，GMP)^[8].因为模态方法简单，易于实现，使柔性多体系统中的柔性分量自由度显著降低，所以是一类应用广泛的模型降阶技术.国内唐嫕璇博士在基于绝对节点坐标法(absolute nodal coordinate formulation，ANCF) 描述的柔性多体系统模型降阶方面开展了较为深入的研究，针对具有大范围运动和大变形的柔性多体系统，提出了一种有效的模态降阶方法——连续局部线性化模态降阶方法^[9]，通过对柔性多体系统进行多次线性化，再对每个线性系统建立各自的降阶模型，最终达到对整个柔体多体系统降阶的目的.但在降阶过程中该方法就模态选取这一问题未进行考虑，因为许多低频模态集并不促进结构和控制系统之间的相互作用^[10].因此，在对柔性多体系统进行模态降阶时，尤其是处于大变形、大旋转状态的柔性多体系统，如果模态选取不当，可能会导致精度问题，甚至求解结果发散.

为了解决上述问题，文献[11]提出了一种基于非线性压缩感知的模态自适应选择方法，它将模态坐标定义为正交模态向量上的稀疏系数；柔性多体系统的动态响应分析定义为以采样后的动力学方程为约束的模态坐标范数优化问题；给出改进的贪婪Gauss-Newton迭代算法(modified greedy Gauss-Newton，MGGN)求解优化问题，可以得到稀疏的模态坐标，选择模态坐标中非零元素对应的模态向量，即完成了模态的自适应选择.但此方法存在局限性，仅应用在浮动坐标法(floating frame of reference formulation，FFRF)^[12]描述的柔性多体系统动态响应分析中.然而，FFRF在对柔性旋转梁仿真时会出现变形随转速无限增大甚至发散这一错误结果，有违柔性旋转梁绕端部旋转时，伴随转速增大而表现出端部横向变形先快速变大、后快速变小的实验现象^[13].

因此，为了满足机械系统的刚性和柔性部件高速运动时，大范围移动和旋转运动导致的柔性部件大的弹性变形需求.本文提出了一种基于ANCF^[14]的模态自适应选择方法.因为ANCF正好适用于发生任意位移的柔性体的大旋转和大变形分析，并且ANCF在数值积分方面具有优势——恒定质量矩阵和没有Coriolis力项的动力学方程，这也保证了刚体结构旋转时刚体惯性的精确建模.

本文主要工作归纳为以下4点：

1) 采用基于连续介质力学的ANCF梁单元对系统进行建模，并将系统的动力学方程划分成两个部分——边界区域相关的部分和内部区域相关的部分.

2) 应用Craig-Bampton法将非线性系统投影到模态空间上^[15]，模态坐标和模态向量被定义为稀疏系数和稀疏基.

3) 根据Latin超立方抽样设计的采样矩阵，减少动力学方程的数量.以采样后的方程作为约束条件来构造模态坐标范数优化问题，用贪婪Gauss-Newton(GGN)算法求解优化问题，得到稀疏的模态坐标(动态响应分析中起主要作用)，选择模态坐标非零元素对应的模态向量，即完成了模态的自适应选择.

4) 最后通过两个实例验证了该模态自适应选择方法应用在柔性多体系统大变形和大旋转运动中的有效性.

1.   基于ANCF的平面梁系统建模

1.1   不含剪切的二维梁单元动力学方程

图 1显示了在全局坐标系中由多个梁单元组成的梁部件k的初始和当前构型，长度为l^k；图 2显示了梁单元i的初始和当前构型，长度为l.梁单元i中轴线上任意一点p^ki在全局坐标系下的位置矢量可以用单元节点坐标和形函数来表示：

图  1  梁部件k的初始构型和当前构型

Figure  1.  Undeformed and deformed configurations of beam component k

下载: 全尺寸图片 幻灯片

图  2  梁单元i的初始构型和当前构型

Figure  2.  Undeformed and deformed configurations of beam element i

下载: 全尺寸图片 幻灯片

其中S^ki为全局坐标系下定义的形状函数，

I∈R^2×2为单位矩阵，L=x^ki/l，x^ki为梁单元i在未变形状态下的局部坐标；r₁^ki和r₂^ki为r^ki向量坐标列阵中的元素；e^ki为单元i的节点坐标向量，

$e_{1}=\left.r_{1}^{k i}\right|_{x^{k i}=0}$和$e_{2}=\left.r_{2}^{k i}\right|_{x^{k i}=0}$是梁单元$i$左侧节点在全局坐标系下的坐标, $e_{5}=\left.r_{1}^{k i}\right|_{x^{k i}=l}$和$e_{6}=\left.r_{2}^{k i}\right|_{x^{k i}=l}$是右侧节点在全局坐标系下的坐标; $e_{3}=\left.\frac{\partial r_{1}^{k i}}{\partial x^{k i}}\right|_{x^{k i}=0}$和$e_{4}=\left.\frac{\partial r_{2}^{k i}}{\partial x^{k i}}\right|_{x^{k i}=0}$是梁单元$i$左侧节点斜率, $e_{7}=\left.\frac{\partial r_{1}^{k i}}{\partial x^{k i}}\right|_{x^{k i}=l}$和$e_{8}=$ $\left.\frac{\partial r_{2}^{k i}}{\partial x^{k i}}\right|_{x^{k i}=l}$是右侧节点斜率.

式(1)对时间求导可以得到单元动能公式，如下:

式中, $\left({ }^{\cdot}\right)$表示对时间求导，ρ^ki，V^ki为梁单元i的密度和体积；$\boldsymbol{M}^{k i}=\int_{0}^{V} \rho^{k i}\left(\boldsymbol{S}^{k i}\right)^{\mathrm{T}} \boldsymbol{S}^{k i} \mathrm{d} V^{k i}$为梁单元i的质量矩阵，是一个常数矩阵且离心力和Coriolis力都为零.

在确立了质量矩阵以后，要考虑广义力的表达式.其中包括外部施加的力和系统柔性部件变形引起的弹性力.由虚功原理可知，作用在梁单元i上任意一点的外力F^ki对单元做的虚功为

式中，δr^ki为广义外力所对应的虚位移；δe^ki为节点虚位移；Q_f^ki为单元节点坐标的广义外力矩阵.对于做大变形、大旋转的柔性梁，应考虑梁单元自重(重力视为分布力)对其变形的影响，虚功可以通过对单元体积求积分得到，进而求得单元广义重力矩阵：

广义弹性力通过虚功原理得到虚功表达式，可以用来定义梁单元i的刚度矩阵；应变能也可以用来定义刚度矩阵.本文考虑通过应变能来获取刚度矩阵表达式，在连续介质力学公式中，单元总应变能U^ki用纵向应变能U_l^ki和横向应变能U_t^ki的线性组合来描述^[16]：

式中, E^ki为弹性模量，A^ki为梁单元截面面积，I^ki为截面惯性矩；根据单元总应变能U^ki的表达式，可以得出刚度矩阵也应该由两部分组成：纵向变形相关的刚度矩阵和横向变形相关的刚度矩阵.先介绍纵向变形相关的刚度矩阵及推导过程，根据几何非线性理论，纵向应变为

为了得到纵向应变更加一般的表达式，将节点坐标向量写成两个向量的和:

其中e_r^ki代表刚体的位移，e_f^ki是e^ki减去e_r^ki的差值,

x和y是任意刚体位移，c和s分别表示cos θ和sin θ，θ表示任意刚体在全局坐标系下旋转的角度；由等式$\left(\boldsymbol{e}_{\mathrm{r}}^{k i}\right)^{\mathrm{T}}\left(\frac{\partial \boldsymbol{S}^{k i}\left(x^{k i}\right)}{\partial x^{k i}}\right)^{\mathrm{T}}\left(\frac{\partial \boldsymbol{S}^{k i}\left(x^{k i}\right)}{\partial x^{k i}}\right) \boldsymbol{e}_{\mathrm{r}}^{k i}=1$，代入到式(9)中，可得到更加一般的纵向应变，下文中$\left(\frac{\partial \boldsymbol{S}^{k i}\left(x^{k i}\right)}{\partial x^{k i}}\right)^{\mathrm{T}}\left(\frac{\partial \boldsymbol{S}^{k i}\left(x^{k i}\right)}{\partial x^{k i}}\right)$用S_l^ki来表示.

根据式(12)和(8)中梁单元i的纵向应变能，与纵向变形有关的刚度矩阵为

由于e_r^ki总是以式(11)的形式存在，所以e_r^ki在式(12)、(13)中是任意的，甚至可以方便的设置为e_r^ki=[0, 0, 1, 0, l, 0, 1, 0]^T.

横向变形有关的刚度矩阵推导如下：通常情况下纵向变形总是比横向变形要小得多，于是可以将纵向变形看作是一个很小的量，在这种情况下纵向变形的变形梯度$\sqrt{\left(\boldsymbol{r}^{k i}\right)^{\prime T}\left(\boldsymbol{r}^{k i}\right)^{\prime}}$约等于常数1，则无限小的弧长$\mathrm{d} s^{k i}=\sqrt{\left(\boldsymbol{r}^{k i}\right)^{\prime \mathrm{T}}\left(\boldsymbol{r}^{k i}\right)^{\prime}} \mathrm{d} x^{k i}$可以写成$\mathrm{d} s^{k i} \simeq \mathrm{d} x^{k i}$，其中$\left(\boldsymbol{r}^{k i}\right)^{\prime}=\frac{\partial \boldsymbol{r}^{k i}}{\partial x^{k i}}$.因此曲率$\kappa^{k i}=\left|\frac{\mathrm{d}^{2} \boldsymbol{r}^{k i}}{\mathrm{d}\left(s^{k i}\right)^{2}}\right| \simeq\left|\frac{\mathrm{d}^{2} \boldsymbol{r}^{k i}}{\mathrm{d}\left(x^{k i}\right)^{2}}\right|$，结合式(1)可以得到

将式(14)代入到U_t^ki的表达式中，得到横向变形有关的刚度矩阵：

它是一个恒定不变的常数矩阵.于是式(8)中与横向变形有关的应变能可以写成如下形式：

已知梁单元i的动能公式、质量矩阵和刚度矩阵，根据Lagrange方程以及虚功原理可以推导出梁单元i的动力学方程如下：

其中$\boldsymbol{C}_{\boldsymbol{e}^{k i}}^{k i}$是单元约束方程C^ki关于单元绝对节点坐标e^ki的Jacobi矩阵，λ^ki是Lagrange乘子向量.

1.2   柔性梁部件的动力学方程

平面梁系统中梁部件k的质量矩阵M^k、广义外力阵Q_f^k以及纵向变形相关的刚度矩阵K_l^k和横向变形相关的刚度矩阵K_t^k，由ANCF梁单元拼接如下：

其中N为单元数，B^ki是与单元相关的Boole矩阵.梁部件k的动力学方程可以通过组装单元的动力学方程(17)得到

其中m_k是约束方程的个数.

在求梁单元纵向变形的刚度矩阵时总是和单元节点坐标相关，因此每次更新耗时巨大，本文将式(13)中K_l^ki的单元绝对节点坐标提到积分外，K_l^ki改写成

其中$\otimes$代表的是矩阵进行Kronecker积运算，n∈R^8×8为单位阵，$\boldsymbol{S}_{\mathrm{L}}=\int_{0}^{l} \boldsymbol{S}_{1}^{k i} \otimes \boldsymbol{S}_{1}^{k i} \mathrm{d} x^{k i}$.那么与纵向变形相关的刚度矩阵K_l^k则为

其中e^k为梁部件k上所有绝对节点坐标，e_r^k为第k个梁部件的刚体位移，N^k为单位阵，它的行、列数与e^k和e_r^k相同，

K_lin^k是固定的，所以每次更新K_l^k特别方便，只需将每个时间步长下的绝对节点坐标代入相乘即可.

1.3   ANCF平面梁系统模态分析方法

为了提高求解效率，需要将动力学模型进行降阶处理.根据文献[15]，通过CMS对ANCF平面梁系统矩阵进行降阶处理，其基本思想是将复杂非线性系统投影到模态空间上，并截取低阶模态构成的系统，忽略高阶模态的影响，从而实现降阶. 因为K_t^k是恒定的，可以假设非线性的纵向弹性力K_l^ke^k是一种伪外力，将其移动到式(19)的右端得到新的动力学方程：

并将系统的绝对节点坐标和动力学方程划分为两个部分——边界区域相关的部分和内部区域相关的部分.动力学方程可以改写成如下形式：

其中，M_bb^k和K_tbb^k代表与边界区域相关的质量矩阵和横向变形刚度矩阵；M_aa^k和K_taa^k代表与内部区域相关的质量矩阵和横向刚度矩阵；M_ba^k，M_ab^k，K_tba^k和K_tab^k代表内部区域坐标和边界区域坐标相互耦合的质量矩阵和刚度矩阵；Q^k=Q_f^k-K_l^ke^k被划分成了边界区域相关的 Q_b^k和内部区域相关的$\boldsymbol{Q}_{\mathrm{a}}^{k}; \boldsymbol{C}_{e_{\mathrm{b}}^{k}}^{k}, \boldsymbol{C}_{e_{\mathrm{a}}^{k}}^{k}$表示约束方程C^k有关边界区域坐标和内部区域坐标的Jacobi矩阵.

因为式(24)等号左侧是线性的，所以可以应用Craig-Bampton方法，将非线性系统投影到模态空间上，内部区域坐标$\boldsymbol{e}_{\mathrm{a}}^{k} \in R^{n_{\mathrm{a}} \times 1}$用一组模态坐标$\boldsymbol{\xi}^{k} \in R^{n_{\mathrm{a}} \times 1}$来表示，T^k为模态变换矩阵，通过$\boldsymbol{p}^{k}=\left\{\left(\boldsymbol{e}_{\mathrm{b}}^{k}\right)^{\mathrm{T}}, \left(\boldsymbol{\xi}^{k}\right)^{\mathrm{T}}\right\}^{\mathrm{T}}$来表示$\boldsymbol{e}^{k}=\left\{\left(\boldsymbol{e}_{\mathrm{b}}^{k}\right)^{\mathrm{T}}, \left(\boldsymbol{e}_{\mathrm{a}}^{k}\right)^{\mathrm{T}}\right\}^{\mathrm{T}}$，具体如下：

其中$\boldsymbol{I}_{\mathrm{b}} \in R^{n_{\mathrm{b}} \times n_{\mathrm{b}}}$是单位阵，n_b是边界区域相关坐标的自由度，$\boldsymbol{\varPhi}_{c}^{k} \in R^{n_{\mathrm{a}} \times n_{\mathrm{b}}}$是约束模态矩阵，见式(27)，n_a是内部区域相关坐标的自由度，$\boldsymbol{\varPhi}_{\mathrm{n}}^{k} \in R^{n_{\mathrm{a}} \times n_{\mathrm{a}}}$是固有模态矩阵，通过求解本征方程(28)得到：

其中$\varOmega_{m}^{k}$和$\boldsymbol{\varPhi}_{n m}^{k} \in R^{n_{\mathrm{c}} \times 1}$分别为梁部件k内部区域的第m个固有频率和它对应的模态向量，n_c是模态向量的维度.

2.   模态自适应选择方法

2.1   全模态稀疏表示内部区域坐标

由于只有少量的模态在模型的动态行为中有很高的参与因素，因此可以将全模态坐标看作是稀疏或近似稀疏.模态坐标定义为稀疏系数，全模态向量定义为稀疏基，将稀疏模态坐标扩展到整个平面梁系统.结合Lagrange乘子向量λ^k，可以得到等式

其中边界区域坐标e_b^k和Lagrange乘子向量λ^k的数量远远少于模态坐标ξ^k.

2.2   设计采样矩阵

为了减少动力学方程的数量，提升数值计算效率，需要对运动方程进行采样.由于模态向量作为稀疏基来表示内部区域坐标，所以采样矩阵A^k作用于内部区域相关的方程

其中Ω^k代表的是需要保留的动力学方程在系统方程中的索引.本文采用Latin超立方体抽样创建Ω^k.采样后的方程(24)以及简化式表示如下：

其中(·)_Ω^k代表的是梁部件k采样以后的矩阵或向量.

2.3   平面梁系统动态响应分析的l₁范数优化问题

本小节提出了基于ANCF的模态自适应选择方法，以ANCF平面梁系统采样后的动力学方程为约束条件，构造的模态坐标l₁范数优化问题，求解每个时间t步长下的优化问题，可以得到稀疏的模态坐标，选择模态坐标非零元素对应的模态向量，即完成了模态的自适应选择.

本小节先介绍优化问题的建立，通过求解动力学方程(19)可以获得梁系统每个时间步长下的动态响应，根据方程(32)和(19)，柔性多体系统的动态响应的非线性映射算子定义如下：

其中η为缩放因子，$\hat{\boldsymbol{\lambda}}^{k}=\eta \boldsymbol{\lambda}^{k}, \boldsymbol{f}\left(\boldsymbol{e}^{k}\right)=\left(\boldsymbol{Q}^{k}\right)_{\varOmega^{k}}-\left(\boldsymbol{K}_{t}^{k}\right)_{\varOmega^{k}} \boldsymbol{e}^{k}$. 采用一阶向后差分公式(backward differentiation formula，BDF)作为时间积分法，

将BDF公式和式(29)代入到式(33)中，平面梁系统的非线性映射算子可以改写成如下形式：

于是构造出ANCF平面梁系统动态响应t时间步长下的l₁范数优化问题：

其中，‖·‖₁代表 1-范数，s^k表示柔性部件k的模态坐标ξ_t^k的稀疏性.

2.4   优化问题的求解

根据文献[17]提出的求解非线性欠定方程组的稀疏解的GGN算法，它可以得到一个非常精确的解，本文采用此算法.此外，由于不同的工况，在模态坐标中会产生不同的幅度.因为将模态坐标的某些元素设置为零的阈值难以确定，所以设置一个非线性映射算子P_A对每个柔性部件保留s^k个绝对值最大的且非零的模态坐标.

为了保证GGN算法的收敛性和提高算法的效率，文献[17]提出了一种求非支持集最大下降的方法.通过求解最大问题，可以确定最大下降方向的指标.

其中‖·‖₂代表 2-范数，Φ_c^k是Φ^k(c)的Jacobi矩阵，L=Φ_c^k(∶, S)由Φ_c^k的列组成；支持集S由已知的支持集 S^k和未知的支持集S^uk组成，S^k由边界区域坐标e_b^k和Lagrange乘子向量$\hat{\boldsymbol{\lambda}}^{k}$在c中的位置索引组成；S^uk是模态坐标ξ^k非零元素在c中的位置索引，Φ_c^k(∶, l)代表的是 Φ_c^k的第l列，I是单位阵，运算符(·)^$\dagger$表示求伪逆阵.GGN算法的执行步骤如图 3所示.

图  3  GGN算法流程图

Figure  3.  The flowchart for the GGN algorithm

下载: 全尺寸图片 幻灯片

图 3中，GGN算法有4个步骤：

步骤1  输入动态响应上一个时间步长的结果$\boldsymbol{e}_{t-1}^{k}, \dot{\boldsymbol{e}}_{t-1}^{k}$，已知的支持集S^k，以及柔性部件k的稀疏系数s^k和采样矩阵A^k.

步骤2  用向前差分公式初始化e_t^k，初始的$\hat{\boldsymbol{\lambda}}_{t}^{k}$设置为零向量.利用模态变换矩阵以及e_t^k和$\hat{\boldsymbol{\lambda}}_{t}^{k}$组成的向量求得c.用非线性映射算子P_A(ξ_t^k)处理模态坐标ξ_t^k，将其代入到c中，更新c.找到模态坐标ξ_t^k非零元素在c中的索引，初始化未知的支持集S^uk.再次通过模态变换矩阵和c，得到绝对节点坐标e_t^k，最后利用BDF求出$\dot{\boldsymbol{e}}_{t}^{k}$和$\ddot{\boldsymbol{e}}_{t}^{k}$.

步骤3  计算Φ^k(c)和定义在式(35)的Jacobi矩阵Φ_c^k，通过求解式(37)获得在非支持集中的最大下降指数l_max.如果f(l_max)小于给定的公差，则将该指数添加到支持集S中.最后迭代值d通过计算求得.

步骤4  更新c，ξ_t^k以及S.计算非线性映射算子Φ^k(c).如果Φ^k(c)超过了给定的误差，则返回步骤3，否则结束迭代.最后得到$\boldsymbol{e}_{t}^{k}, \dot{\boldsymbol{e}}_{t}^{k}, \dot{\boldsymbol{e}}_{t}^{k}$和$\hat{\boldsymbol{\lambda}}_{t}^{k}$，作为t步长下的动态响应结果.

3.   数值结果

本文的算例都是在处理器配置为Intel(R) Core(TM) i7-7500U CPU@2.70 GHz和存储器配置为12 GB的计算机上使用MATLABⒸR2016a执行的.

3.1   自由落体的单摆运动

本文通过两个例子来验证基于ANCF的模态自适应选择方法的有效性.将自由落体的单摆运动作为第一个例子.如图 4，梁通过旋转接头与地面连接，在初始状态下保持水平，初速度为零.

图  4  自由落体单摆

Figure  4.  The free falling pendulum

下载: 全尺寸图片 幻灯片

梁的几何参数和材料参数见表 1.

表  1  单摆的几何与材料参数

Table  1.  Geometry parameters and material parameters of the pendulum

parameter	value
length l/m	1
square sectional area A/m2	4×10-4
Young’s modulus E/Pa	7×105
density ρ/(kg/m3)	7.2×103
moment of inertia I/m4	1.333×10-8

下载: 导出CSV 

| 显示表格

用160个ANCF不含剪切的二维梁单元，总自由度为644来表示自由落体的单摆运动，总模拟时间和时间步长分别设置为1.1 s和0.001 s.单摆的稀疏度设置为100，采样数设置为308.公差ε设为1×10^-11.图 5给出了仿真过程中若干时刻单摆梁中线的构型，用提出的方法和传统的ANCF进行比较，单摆自由端的横向位移如图 6所示.显然，由提出的方法得到的仿真结果与传统的ANCF高度一致.

图  5  不同时刻柔性单摆构型

Figure  5.  Configurations of the free falling pendulum at different moments

下载: 全尺寸图片 幻灯片

图  6  单摆自由端的垂直位置

Figure  6.  Vertical positions of the free end of the pendulum

下载: 全尺寸图片 幻灯片

此外，图 7展示了提出的方法在不同时间步长下对单摆模态坐标的自适应选择，实线为模态坐标的值，黑圈表示被挑选出的100个具有重大贡献的模态，即选择100个绝对值最大的模态坐标.但是在不同的时间步长下选择方式是不同的，换句话说，该方法实现了在每个时间步长下动态选择模态.

图  7  提出的方法在不同时间步长下对模态的自适应选择(单摆)

Figure  7.  Adaptive selection of modal coordinates with the proposed method for different time steps(pendulum)

下载: 全尺寸图片 幻灯片

本文提出的模态自适应选择方法计算时间比传统ANCF要短，前者总的计算时间为531.855 s，后者为787.071 s，效率至少提高了30%，若减少采样数，效率可以提升更多，详情见表 2.

表  2  单摆传统ANCF和所提出方法的计算效率(单位: s)

Table  2.  Computation efficiency of the ANCF and the proposed method for pendulums (unit: s)

model	matrix operation	updated Jacobian matrix, stiffness matrix and residue etc	total time
ANCF	176.411	598.334	787.071
proposed	118.341	412.514	531.855

下载: 导出CSV 

| 显示表格

3.2   平面3-RRR并联机器人

平面3-RRR并联机器人示意图如图 8所示，图中移动平台的局部坐标系固定在平台的中心点上，全局坐标系在图中右侧.驱动杆与机架(底座)相连，从动杆与移动平台相连^[18].

图  8  全局坐标系下的3-RRR并联机器人示意图

Figure  8.  Geometry and global coordinates of the 3-RRR mechanism

下载: 全尺寸图片 幻灯片

机构的几何参数和材料参数见表 3.为了扩大连杆弹性变形的影响，设置了小于铝的弹性模量.在仿真中，3-RRR并联机器人的移动平台所预设的运动轨迹为一个半径为0.1 m的圆:

表  3  机构的几何与材料参数

Table  3.  Geometry parameters and material parameters of the mechanism

material parameter	driving link	passive link		member	length l/m
thickness T/m	0.01	0.005		driving link	0.245
width W/m	0.03	0.01		passive link	0.242
Young’s modulus E/Pa	2.01×1011	7×108		moving stage	0.112
density ρ/(kg/m3)	2.7×103	2.7×103		fixed stage	0.400

下载: 导出CSV 

| 显示表格

ω是一个角速度，值为2π rad/s.在模拟开始时，平台的初始位置为x=0.1 m，y=0 m，z=0 rad，初始速度为$\dot{x}=0 \mathrm{~m} / \mathrm{s}, \dot{y}=0.1 \pi \;\mathrm{m} / \mathrm{s}, \dot{z}=0 \;\mathrm{rad} / \mathrm{s}$.

用自然坐标描述机构驱动杆的运动，驱动杆和移动平台视为刚性部件，移动平台用x, y方向上的位移和x轴正方向上的转角θ表示其运动，用100个ANCF不含剪切的二维梁单元，即自由度为404来表示机构中一个从动杆件(柔性部件)的运动，机构总的自由度为1 227，总模拟时间和步长分别设置为2 s和0.001 s.每一个从动杆件的稀疏度设置为80，采样数设置为200.公差ε设为1×10^-11.

通过对3-RRR并联机器人的运动学计算，得到了开始动态分析所需的初始位置和速度.图 9显示了不同时间的机构构型.可以看出，驱动杆件一直保持未变形，而柔性的从动杆件在平台移动时有明显的变形.

图  9  在不同时刻的机构构型

Figure  9.  Configurations of the mechanism at different moments

下载: 全尺寸图片 幻灯片

用所提出的方法和传统的ANCF进行比较^[19].通过对3-RRR并联机器人的运动学计算，得到了起始时刻动力学分析所需的初始位置和速度.图 10为机构移动平台在x，y方向上的位移和转角θ.显然，仿真结果是高度一致的.

图  10  移动平台在x，y方向上的位移以及转角θ

Figure  10.  Rotation angle θ and displacements in directions x and y of the platform

下载: 全尺寸图片 幻灯片

用提出的方法求解机构左端从动杆件具有重大贡献的模态，实现了模态的自适应选择，如图 11所示.在图中，自适应地选择了80个绝对值最大的模态坐标.

图  11  提出的方法在不同时间步长下对模态的自适应选择

Figure  11.  Adaptive selection of modal coordinates with the proposed method for different time steps

下载: 全尺寸图片 幻灯片

本文所提出的模态自适应选择方法的计算时间比传统ANCF要短，总的计算时间分别为421.281 s和1 390.966 s，所用时间是传统ANCF的30.03%，计算效率提升显著, 详情见表 4.

表  4  传统ANCF和所提出方法的计算效率(单位: s)

Table  4.  Computation efficiency of the ANCF and the proposed method (unit: s)

model	matrix operation	updated Jacobian matrix, stiffness matrix and residue etc	total time
ANCF	566.283 6	766.148	1 390.966
proposed	184.131	223.639	421.281

下载: 导出CSV 

| 显示表格

4.   结论

本文针对柔性大变形系统提出了基于ANCF的模态自适应选择方法.利用两个实例证明了该方法的有效性，与传统的ANCF的数值分析结果高度吻合并且计算效率至少提升了30%，在误差允许的范围内，减少采样数，计算效率更高.另外本文还改进了刚度矩阵的计算方式，提高了每次迭代更新的效率.虽然通过两个实例证实了该方法的有效性，尤其是在描述由离心力引起的大变形模型时，但在分析轴向变形很大的模型上有些许不足.本文为了应用Craig-Bampton方法，将非线性系统映射到模态空间上.使用式(23)对动力学方程等号左端进行线性化处理，这在一定程度上降低了纵向变形与横向弯曲的耦合性，使得方法存在一定的局限性.后续将针对这一问题展开深入研究，改进方法的适用性，通过弹性力高耦合的实例做进一步的验证.