Bennett连杆机构是一种具有单自由度的特殊空间四连杆机构^[1]。CHEN和YOU^[2]提出了一种Bennett机构的衍生形式，其具有完全折叠与完全展开的性质；LI^[3]提出了一种更普遍的Bennett机构衍生形式，能够完全展开为空间四边形。以Bennet连杆机构衍生形式为基本单元，通过一定的连接方式可拓展成覆盖更大平面的网格组合。Bennett连杆机构具有完全收拢状态占用空间小、便于运输储存等优点，在建筑工程^[3]、航天技术^[2]、应急救灾^[4]和医疗器械^[4]等领域有广阔的应用前景。

现有的机构动力响应研究大多集中在具有确定参数的机构中。然而，实际应用中机构参数等具有不确定性，如材料不确定性、几何不确定性和荷载不确定性。因此，参数确定的机构动力学分析无法反映不确定性参数对机构动力行为的影响。目前描述不确定性参数的方法主要有两类，即概率方法^[5]和非概率方法^[6]。概率方法通常用于描述已知参数概率密度函数的不确定性分析，而非概率方法主要用于描述参数有界但概率密度函数未知的不确定性分析。概率方法包括蒙特卡洛模拟法^[7]、响应面法^[8]、非统计方法^[5]等，其中非统计方法包括微分分析方法^[5]、多项式混沌方法等。实际应用中机构参数的完整概率信息往往难以获得，使用概率方法具有一定的局限性。

在非概率分析中存在许多不确定性分析方法。WEI^[9]等提出了一种分析框架，利用蒙特卡罗采样和优化方法量化具有随机不确定性和区间不确定性的人行桥动力响应。SUN^[10]等将Karhunen-Loève展开和Kriging模型结合，提出了一种具有混合不确定性的含平面间隙铰柔性机械臂的动力响应分析方法。除上述方法外，区间方法是一种确定含区间参数的多项式响应上下界高效非概率算法之一，广泛应用于分析含区间参数的不确定性问题。Chebyshev多项式方法通过求解区间多项式函数可直接确定区间响应值的上下界，是一种有效的区间计算方法。XIANG^[11]等基于Chebyshev多项式，提出一种含有区间不确定性的含间隙机构动力响应的不确定性分析方法。随机分析中常用的代理模型，如Kriging模型和混沌多项式，均对采样点数据存在一定的要求，而Chebyshev多项式则是在一定区间内选取采样点，对数据本身无要求。在实际应用中Bennett机构较难获得完整参数概率分布，因此本文将利用Chebyshev多项式量化Bennett机构响应不确定性。

目前Bennett机构的研究大多集中于机构的运动学分析^[12]，对Bennett机构的动力学分析较少。Bennett机构的动力学分析涉及机构运动和结构变形的耦合。有限元法由于集成的整体刚度矩阵通常是奇异的，难以有效处理机构运动问题。非线性有限元法中的几何精确方法^[13]较难更新旋转参数和处理奇异性问题。绝对节点坐标法^[14]计算成本较高且可能难以收敛。多体动力学可以模拟机构运动过程，进行位移、速度等响应分析^[15]；另一方面，多体动力学也可引入考虑摩擦、接触等非线性分析模型，进行结构加速度、内力等动力响应分析。拉格朗日方程在虚位移原理和达朗贝尔原理的基础上引入广义坐标建立不含约束力的质点系动力学方程。牛顿欧拉法应用达朗贝尔原理，将动力学问题转化为牛顿—欧拉形式的静力学平衡问题。凯恩法建立在分析力学的基础上，以广义速度为自变量，引入了偏速度、偏角速度、广义主动力和广义惯性力的定义，建立了代数方程形式的动力学方程。Bennett机构动力响应不确定性分析是建立在动力学分析基础上的，有必要采用有效的动力分析方法保证结果可靠性。

有限质点法(finite particle method, FPM)是基于向量式力学发展的、面向结构行为的分析方法。该方法可用于结构的复杂行为分析^{[16 − 17]}，并已被成功应用于动不定结构^[18]、可展开结构^[19]和含间隙机构^[20]的分析。有限质点法以质点为基本要素，采用显式积分方案，不需要集成刚度矩阵，避免了刚度矩阵奇异。在机构分析中，可以考虑机构位移和连杆变形的耦合，构建各种复杂机构或系统的动力学模型^[18]。因此结合Chebyshev多项式和有限质点法，可有效开展Bennett机构全过程动力响应不确定性分析。

本文基于Chebyshev多项式和有限质点法，对Bennett机构展开全过程进行动力响应不确定性分析。首先介绍Bennett机构的几何特征；其次给出Bennett机构建模方法和相应的有限质点法分析方法；然后介绍基于Chebyshev多项式的区间方法，并提出一种易与有限质点法结合的非侵入式算法，建立机构参数与动力响应之间的关系；最后通过算例验证方法的有效性，并研究不确定性参数对Bennett机构动力响应的影响。

1   Bennett机构及其衍生形式

Bennett连杆机构是一种特殊的单自由度空间过约束四连杆机构^[1]。它由4个连杆和4个转动铰组成，每个转动铰的转轴方向和与其相连的两个连杆垂直，如图1(a)所示。为保证机构为单自由度，需满足：

图  1  Bennett机构及其衍生形式

Figure  1.  Bennett linkage and its alternative form

下载: 全尺寸图片 幻灯片

式中：${a_{ij}}$为连接节点$i$和$j$连杆的长度；$a$和$b$为机构连杆长度；${\alpha _{ij}}$为节点$i$处转动铰转轴方向向量${{\boldsymbol{e}}_{{\text{r}}i}}$到节点$j$处转动铰转轴方向向量${{\boldsymbol{e}}_{{\text{r}}j}}$的旋转角度；${\theta _i}$为与节点$i$相连连杆到转轴方向向量${{\boldsymbol{e}}_{{\text{r}}i}}$的旋转角度。若为等边Bennett机构，则满足$a = b = l$，$\alpha + \beta = \pi$，$\tan \dfrac{{{\theta _1}}}{2}\tan \dfrac{{{\theta _2}}}{2} = \dfrac{1}{{\cos \alpha }}$。

上述等边Bennett机构不具有完全展开和完全折叠的性质。CHEN和YOU^[2]提出了一种Bennett机构衍生形式来实现机构的完全折叠与完全展开。如图1(b)所示，${{{A}}_0}{{{B}}_0}{{{C}}_0}{{{D}}_0}$为Bennett机构，ABCD为其对应的衍生形式。在CHEN和YOU^[2]研究的基础上引入参数${\mu _{\text{B}}}$，可形成一种可完全展开为曲面的Bennett机构衍生形式，其物理模型如图1(c)所示。

2   Bennett机构有限质点法分析

在有限质点法中，Bennett机构被离散成质点和单元。质点代表机构的位置、质量、受力、变形和边界条件，质点间的相互作用关系通过单元表示。本节将给出Bennett机构的有限质点法分析模型、质点运动方程以及各单元计算方法。

2.1   Bennett机构简化分析模型

目前常用的机构建模方法包括精细建模和简化建模。精细建模将构件离散为实体单元，可以考虑机构尺寸和接触约束，但实体单元数量较多导致计算成本较大。简化建模将构件离散为梁单元，具有建模方便、计算成本低等优点。

本文给出一种Bennett机构简化建模方法，如图2所示，其中标识了模型中的质点和各类单元。质点的位置根据Bennett机构物理模型的几何尺寸计算得到。一根连杆由若干个梁单元模拟；连杆截面通过在截面处虚设梁单元(下文简称虚梁单元)模拟；转动铰由转动铰单元模拟；展开驱动由扭簧单元模拟；连杆截面处的接触通过面接触单元模拟。

图  2  Bennett机构简化计算模型

Figure  2.  Simplified computational model of Bennett linkage

下载: 全尺寸图片 幻灯片

2.2   质点运动方程及相关单元

在有限质点法中，质点运动遵循牛顿第二定律。质点有3个平动自由度和3个转动自由度，其运动方程可表示为：

式中：$m$为质点质量；d为质点平动位移；${{\boldsymbol{F}}^{{\mathrm{ext}}}} = {\left[ {{F_x},{F_y},{F_{\textit{z}}}} \right]^{\mathrm{T}}}$为质点外力；$\displaystyle\sum {\boldsymbol{f}} = {\left[ {\displaystyle\sum {f_x},\displaystyle\sum {f_y},\displaystyle\sum {f_{\textit{z}}}} \right]^{\mathrm{T}}}$为梁单元传递至质点的内力；${{\boldsymbol{F}}^{{\mathrm{dmp}}}} = - \zeta m\dot {\boldsymbol{d}}$为阻尼力，$\zeta$为质量比阻尼系数，表示为单位质量下的阻尼系数；$\theta$为质点角位移；${{\boldsymbol{M}}^{{\mathrm{ext}}}} = {\left[ {{M_x},{M_y},{M_{\textit{z}}}} \right]^{\mathrm{T}}}$为质点外力矩；$\displaystyle\sum {\boldsymbol{m}} = {\left[ {\displaystyle\sum {m_x},\displaystyle\sum {m_y},\displaystyle\sum {m_{\textit{z}}}} \right]^{\mathrm{T}}}$为梁单元传递至质点的内力矩；${{\boldsymbol{M}}^{{\text{dmp}}}} = - \zeta {\boldsymbol{I\dot \theta }}$为阻尼力矩。有限质点法采用显式积分方案，可由中心差分法求解。n+1时刻的${{\boldsymbol{d}}^{n + 1}}$和${{\boldsymbol{\theta }}^{n + 1}}$可根据n−1和n时刻质点的${{\boldsymbol{d}}^{n - 1}}$、${{\boldsymbol{\theta }}^{n - 1}}$、${{\boldsymbol{d}}^n}$和${{\boldsymbol{\theta }}^n}$迭代计算。

梁单元^[20]基于欧拉梁原理描述两质点间相互作用关系。通过虚拟逆向运动扣除梁单元刚体位移，从而获得纯变形。纯变形可分解为轴向变形${\it{\Delta} _{\rm{e}}}$、扭转变形$\Delta \varphi _{{\text{B}}x}^{}$、弯曲变形$\Delta \varphi _{{\text{A}}i}^{}$和$\text{Δ}{\varphi }_{\text{B}i}\left(i=y, {\textit{z}}\right)$共6个分量。随后根据纯变形计算梁单元内力增量^[20]。最后将内力增量与初始时刻质点内力叠加得到当前时刻的总内力大小。虚梁单元设置在连杆端部截面处，用于模拟Bennett机构截面尺寸。虚梁单元密度、弹性模量的设置需要考虑其对机构运动和计算临界时间步长的影响。

转动铰单元^[20]用于模拟Bennett机构中构件间的转动铰。转动铰单元由轴承和轴颈组成，仅有绕旋转轴转动一个自由度，轴承和轴颈分别通过一个质点模拟。两质点的平动自由度相互耦合，具有相同的位置、速度和加速度。两质点与转动轴垂直方向的自由度相互耦合，而沿转动轴方向的自由度相互独立。

扭簧单元用于模拟驱动Bennett机构展开的扭簧。扭簧单元的作用力通过评估两个截面之间二面角的变化来计算。首先计算截面法向量${{\boldsymbol{n}}_1}$、${{\boldsymbol{n}}_2}$和转轴方向向量${{\boldsymbol{e}}_{\text{r}}}$，随后利用${\text{atan}}2$函数计算得到截面间的二面角。通过扭簧单位长度刚度、扭簧长度和截面间二面角计算得到弯矩。将扭簧产生的弯矩转化为力作用到对应质点上，进而驱动机构展开。DONG等^[21]推导了具体计算方法，此处不再赘述。

接触单元用于模拟Bennett机构完全展开或完全闭合时机构的端部截面或侧面的接触作用。接触单元中截面间二面角的计算方法与扭簧单元相同，此处不再赘述。其中的接触状态判定说明如下：设定截面角度预判阈值${\theta ^{{\text{tol}}}}$，当二面角角度小于阈值时，添加预判标记，以标识在接下来的运动中可能发生接触；在之后的迭代计算中，若二面角大于$\pi$且存在预判标记，则表明二面角发生了突变，两截面发生接触。发生接触后，可采用罚函数法计算约束力矩。

3   基于Chebyshev多项式的动力响应近似方法

3.1   Chebyshev多项式

Chebyshev多项式是一种常用的近似连续多项式的方法。已经证明，Chebyshev多项式的逼近精度比其他大多数类型的截断级数更高^{[6, 22]}。此外，Chebyshev区间方法不需要计算区间函数的高阶导数，上下界可通过求解区间多项式(也称为Chebyshev代理模型)确定，区间计算的高估现象也可有效控制^[22]。文献[22]给出了Chebyshev多项式方法的细节，本节简要介绍Chebyshev多项式原理。

对于一维问题，Chebyshev多项式的n次多项式可表示为^[22]：

式中，变量$\xi \in [ {\underline {\chi } ,\overline \chi } ]$，$\theta \in [ {0,\pi } ]$。基于Chebyshev多项式，$\xi$在区间$[ {\underline {\chi } ,\overline \chi } ]$内的任一连续函数$f( \xi )$可被近似为：

式中：${f_i}$为常系数；$n$为Chebyshev多项式的最高阶。对于$\xi \in \left[ { - 1,1} \right]$，Chebyshev多项式满足正交性：

基于以上Chebyshev多项式正交性，常系数${f_i}$可由Mehler 积分方法计算^[22]：

式中：$m$为计算常系数的插值点个数，应不小于$n + 1$。为方便计算，在本文中$m = n + 1$；${\xi _j}$为第$j$个插值点(即采样点)，满足：

由式(7)可以看出，常系数是不同插值点处函数值和多项式值的线性组合。

多维Chebyshev多项式可用多个一维Chebyshev多项式的张量积表示。对于变量${\xi _i} \in [ {{{\underline {\chi } }_i},{{\overline \chi }_i}} ]$，$i = 1,2, \cdots ,k$，$k$维Chebyshev多项式为：

式中：下标${a_i} = 0,1,2, \cdots ,n$代表Chebyshev多项式的阶次；${\theta _i} = {\text{arccos}}\left( {\dfrac{{2{\xi _i} - ( {{{\overline \chi }_i} + {{\underline {\chi } }_i}} )}}{{{{\overline \chi }_i} - {{\underline {\chi } }_i}}}} \right) \in [ {0,\pi } ]$；${\boldsymbol{\xi }} = [ {\xi _1},{\xi _2}, \cdots , {\xi _k} ]^{\text{T}}$为$k$维向量。与式(5)类似，一个$k$维连续函数$f( {\boldsymbol{\xi }} )( {{\boldsymbol{\xi }} \in [ {{\underline {\bf\textit{χ}}},\overline {\bf\textit{χ}}} ]} )$可被近似表示为：

式中：$l$为下标${i_1},{i_2}, \cdots ,{i_k}$中为0的个数；系数${f_{{i_1}{i_2} \cdots {i_k}}}$可通过式(11)计算：

其中，上述采样点的计算公式为：

对于含有多个参数的Bennett机构系统，可利用上述Chebyshev多项式近似方法建立不确定参数与Bennett机构系统动力响应之间的关系，得到系统响应的上下边界。

3.2   机构动力响应近似方法

本节将上述Chebyshev多项式近似方法与FPM结合，提出适用于机构动力响应不确定性分析的非侵入式算法。

对于含有区间不确定参数的结构系统，FPM的质点运动方程式(3)可改写成：

其中$\zeta$是系统的确定参数向量，$[ \varsigma ] \in [ {\underline {\varsigma } ,\overline \varsigma } ]$. 是区间参数向量，其定义与3.1节中区间变量相同。在本文中，区间参数向量可包含Bennett机构连杆长度、弹性模量等参数。

根据式(7)和式(10)可将区间不确定参数转化成${\left( {n + 1} \right)^k}$个样本点，其中$n$为Chebyshev多项式阶数，$k$为区间变量个数。式(13)可分解成${\left( {n + 1} \right)^k}$个含不确定参数分量等式分别求解。通过上述方法，区间不确定参数问题转换成多个确定参数问题。利用上述确定参数分别建立机构计算模型，通过FPM计算其系统动力响应；随后引入时间参数$t$，通过Chebyshev多项式拟合不确定参数与系统动力响应的全过程动力响应关系：

式中：${\boldsymbol{\xi }}$为系统区间参数；${\boldsymbol{F}}\left( {{\boldsymbol{\xi }},t} \right)$为系统动力响应；${{\boldsymbol{F}}_{{i_1}{i_2} \cdots {i_k}}}\left( t \right)$为Chebyshev多项式的常系数向量。由于采样点已确定，常系数向量是仅与$t$(样本点动力响应)相关的函数；相应的${T_{{i_1}{i_2} \cdots {i_k}}}\left( {\boldsymbol{\xi }} \right)$仅与采样点有关，与$t$(样本点动力响应)无关。常系数向量${{\boldsymbol{F}}_{{i_1}{i_2} \cdots {i_k}}}\left( t \right)$可通过(11)计算：

根据以上基于Chebyshev多项式和FPM所建立的代理模型，动力响应边界可由区间运算得到。由于${T_{{i_1} \cdots {i_k}}}\left( {{{\boldsymbol{\xi }}_s}} \right)$为多个余弦函数相乘，满足${T_{{i_1} \cdots {i_k}}}\left( {{{\boldsymbol{\xi }}_s}} \right) \in \left[ { - 1,1} \right]$，即Chebyshev多项式可表示为余弦函数，结合式(14)和式(15)可计算得到动力响应边界：

式中，${i_1},{i_2}, \cdots ,{i_k}$不能全为0。

上述机构动力响应近似方法的计算流程如下，图中$s$代表第$s$个样本向量，$N$代表样本向量数，$T\mathrm{_{fin}}$代表模拟时长。

图  3  机构动力响应近似方法计算流程

Figure  3.  Flowchart of approximation method of dynamic response of mechanisms

下载: 全尺寸图片 幻灯片

4   数值算例

4.1   算例验证

本节对两个典型机构(滑块曲柄和双摆机构)进行不确定性分析，并与论文^[6]的数据进行对比，验证本文所提出方法的有效性。

4.1.1   曲柄滑块

本节将曲柄滑块中的曲柄长度作为不确定参数进行不确定分析。图4绘制了曲柄滑块的模型示意图，其中点A、B和C分别为曲柄、连杆和滑块的中心。曲柄滑块的参数与Wu等^[6]的研究相同，参数如表1所示。${\theta }_{1}和{\theta }_{2}$分别表示曲柄、连杆局部坐标系与全局坐标系的夹角。曲柄中心、连杆中心和滑块在全局坐标系下的坐标分别为$\left( {{x_1},{y_1}} \right)$、$\left( {{x_2},{y_2}} \right)$和$\left( {{x_3},0} \right)$。在曲柄滑块初始状态，$\left[ {{\theta _1},{\theta _2}} \right] = \left[ {0,0} \right]$，机构在自重和与滑块相连的弹簧阻尼器驱动下开始运动。当${\theta }_{1}和{\theta }_{2}$均为0时，弹簧阻尼器中的弹簧力与阻尼力均为0。

假设曲柄连杆长度${l_1}$的不确定度为其设定值的1%，则有：

式中，${\xi _1} \in \left[ { - 1,1} \right]$。

分别采用基于3阶和5阶Chebyshev多项式的机构动力响应近似方法对曲柄滑块系统进行不确定分析，并利用区间运算得到代理模型的边界。其中，3阶和5阶Chebyshev多项式的采样点分别为4和6。建立FPM模型，曲柄和和连杆分别由5个和10个梁单元模拟。总模拟时间为2 s，时间步长设为$1 \times {10^{ - 6}}$ s。

图  4  曲柄滑块示意图

Figure  4.  Schematic of slider crank.

下载: 全尺寸图片 幻灯片

表  1  曲柄滑块参数

Table  1.  Parameters of slider crank mechanisms

参数	数值
曲柄长度 ${l}_{1} /{\mathrm{m}}$	0.15
连杆长度 ${l}_{2}/\text{m}$	0.56
曲柄质量 ${m_1}$ /kg	0.37
连杆质量 ${m_2}$ /kg	0.77
滑块质量 ${m_3}$ /kg	0.45
弹簧阻尼器阻尼系数 $c$ /(N/(m/s))	1
弹簧阻尼器弹性模量 $k$ /(N/m)	5
扭矩 $\tau$ /(Nm)	-0.5

下载: 导出CSV 

| 显示表格

图5绘制了滑块的位移时程曲线。可以看出，本文计算得到的响应边界与文献[6]结果基本吻合，证明了该方法求解含有单个不确定参数模型的有效性。此外，基于3阶Chebyshev多项式的模拟结果与5阶Chebyshev结果基本一致。已计算3阶Chebyshev多项式的模拟结果与5阶Chebyshev结果相对文献(Scan method)计算结果的相对误差， 图6和图7分别绘制了下边界和上边界的相对误差时程曲线。3阶Chebyshev多项式与5阶Chebyshev的误差曲线均基本重合，误差保持在5%以内。其中，3阶Chebyshev多项式的上下边界模拟结果误差最大值分别为3.68%和5.06%，5阶Chebyshev多项式的上下边界模拟结果误差最大值分别为3.76%和4.74%。因此综合考虑计算精度和计算成本，选择3阶Chebyshev多项式作为代理模型进行计算。

图  5  滑块位移响应区间

Figure  5.  Interval of displacement response of crank

下载: 全尺寸图片 幻灯片

图  6  不同阶数Chebyshev多项式拟合结果下边界误差时程曲线

Figure  6.  Time history of error of lower boundary for different orders of Chebyshev polynomial fitting results

下载: 全尺寸图片 幻灯片

图  7  不同阶数Chebyshev多项式拟合结果上边界误差时程曲线

Figure  7.  Time history of error of upper boundary for different orders of Chebyshev polynomial fitting results

下载: 全尺寸图片 幻灯片

4.1.2   双摆机构

本节将双摆机构的几何参数连杆长度作为区间参数进行不确定分析。图8绘制了双摆机构的初始形态，该机构从静止状态在重力作用下运动。双摆机构的参数和WU等^[6]的研究相同。两根连杆长度、质量相同，为${l_1} = {l_2} =$1 m，${m_1} = {m_2} =$1 kg。重力加速度$g={\left[ 0\;\; 9.81\;\;0 \right]}^{\text{T}}\text{ m}/{\text{s}}^{2}$。${\theta }_{1}和{\theta }_{2}$分别表示连杆局部坐标系与全局坐标系间的夹角，上下杆中心坐标分别为$\left( {{x_1},{y_1}} \right)$和$\left( {{x_2},{y_2}} \right)$，初始角度和初始角速度设为$[{\theta }_{1},{\theta }_{2},\dot{{\theta }_{1}},\dot{{\theta }_{2}}]=\left[\dfrac{\pi }{4},\dfrac{\pi }{3},0,0\right]$。

假设双摆两连杆长度的不确定性为其设定值的5%，则有

式中，${\xi _1},{\xi _2} \in \left[ { - 1,1} \right]$。

图  8  双摆示意图

Figure  8.  Schematic of double pendulum.

下载: 全尺寸图片 幻灯片

采用基于3阶Chebyshev多项式的机构动力响应近似方法对双摆机构进行不确定分析，并利用区间运算得到代理模型的边界。采样点个数为${6^2} = 36$。建立FPM模型，两根连杆均由10个梁单元模拟。总模拟时间为10 s，时间步长设为$3 \times {10^{ - 6}}$ s。

图9绘制了上杆中心的水平位移时程曲线。可以看出，3阶和5阶Chebyshev多项式分析结果基本吻合，与文献结果差异较小^[6]，证明了该方法能有效求解含有两个不确定性参数的机构动力响应。

图  9  上杆中心的水平位移响应区间

Figure  9.  Interval of horizontal displacement response of upper rod center.

下载: 全尺寸图片 幻灯片

关于Chebyshev多项式阶数的选取，一般而言可以通过将Chebyshev多项式的拟合结果与精确解Scan method的结果进行比较，判断拟合精度是否满足要求，并在保证拟合精度的条件下选取尽可能低的Chebyshev多项式展开的阶数，以提高计算效率。对于机构动力响应不确定性分析，当不确定参数较少时，如1个或2个时，3阶Chebyshev多项式和5阶Chebyshev多项式拟合结果基本相同；当不确定参数较多时，需要选取高阶数Chebyshev多项式进行拟合。在一般机构的简单分析中，选择3阶Chebyshev多项式即可兼顾拟合精度和计算效率的要求。

4.2   Bennett机构动力响应不确定性分析

本节利用上述机构动力响应近似方法对Bennett机构进行动力响应不确定性分析。所讨论的Bennett机构由4根连杆、4个扭簧和4个转动铰组成。Bennett机构在扭簧驱动下展开，模型参数设置如表2所示。

表  2  Bennett机构模型几何与材料参数

Table  2.  Geometrical and material parameters of Bennett linkage

参数	数值
连杆长度L	0.45
截面边长a/m	0.035
$\omega$/($^\circ$)	54.7
${\mu _{\text{B}}}$/($^\circ$)	0
$\lambda$/($^\circ$)	45.0
连杆弹性模量$E$/Pa	$1 \times {10^9}$
连杆泊松比$\nu$	0.3
连杆密度$\rho$/(${\text{kg}}/{{\text{m}}^3}$)	800
虚梁截面边长${a_{{\text{vb}}}}$/m	0.05
虚梁弹性模量${E_{{\text{vb}}}}$/Pa	$1 \times {10^9}$
虚梁泊松比$\nu$	0.3
虚梁密度$\rho$/(${\text{kg}}/{{\text{m}}^3}$)	2
单位长度接触刚度${K^{\rm{C}}}$/(${\text{N}}/{\text{rad}}$)	$1 \times {10^4}$
单位长度扭簧刚度${K^{\rm{F}}}$/(${\text{N}}/{\text{rad}}$)	5

下载: 导出CSV 

| 显示表格

该机构的FPM计算模型如图10所示。单根连杆由8个梁单元模拟，机构共被离散为76个质点(红色)、32个梁单元(蓝色)、40个虚梁单元(绿色)、4个转动铰单元(黄色)、4个扭簧单元和4个接触约束单元，如图10中编号所示。机构总模拟时间为1 s，时间步长为$5 \times {10^{ - 8}}$ s。为了便于讨论不确定参数对Bennett机构动力响应的影响，本文将模型质量比阻尼系数$\zeta$设为0。

图  10  Bennett机构的计算模型

Figure  10.  Computational model of Bennett linkage

下载: 全尺寸图片 幻灯片

4.2.1   连杆长度不确定影响

本节讨论Bennett机构连杆长度不确定时的动力响应。假设所有连杆长度完全相同，不确定性为1%，表示为：

式中，$\xi = \left[ { - 1,}\;\;1 \right]$为区间参数。本节以3阶Chebyshev多项式建立代理模型，并使用区间运算对代理模型进行不确定性分析，所需采样点数量为4。

图11绘制了端点B、D间的距离时程曲线。在初始阶段，机构完全闭合，BD距离为0；从0到约0.28 s，机构在扭簧驱动下展开，BD距离逐渐增大；随后，当机构完全展开时，连杆截面发生碰撞，此时BD距离达到最大；由于连杆截面碰撞为弹性碰撞，截面反弹，机构发生收拢，接近完全闭合状态后，机构重新展开。

图  11  BD距离时程曲线

Figure  11.  Time history of distance BD

下载: 全尺寸图片 幻灯片

图11还绘制了BD距离时程曲线的上下界。通过本文方法计算出的响应值上下边界完全包络连杆长度确定时的Bennett机构BD距离时程曲线。在机构第一次展开过程中，BD距离上下边界较接近，最大差值为0.068 m，连杆长度不确定性影响较小；在第二次展开过程中，BD距离的上下边界区间明显增大，最大值为0.183 m，占对应时刻确定参数模型BD距离的87.1%，连杆长度不确定性影响显著增大，机构位移响应不确定性显著增强。这是由于区间算法的包裹效应(Wrapping effect)^[6]导致。

图12绘制了端点B处$z$方向的速度时程曲线。在第一次连杆截面碰撞前，速度时程曲线没有明显振荡，上下边界曲线差距较小。在第一次完全展开截面碰撞后，速度时程曲线振荡显著加剧，上下边界曲线区间增大。在0.68 s之后，机构开始第二次展开过程，速度时程曲线上下边界差距进一步增大。在[0.85 s, 0.9 s]内，曲线振荡较弱，上下边界曲线最大差值为5.64 m/s。在[0.9 s, 1.0 s]内，曲线振荡显著增强，上下边界曲线最大差值为5.77 m/s。

图  12  B点$z$方向速度时程曲线

Figure  12.  Time history of z-direction velocity of point B

下载: 全尺寸图片 幻灯片

需要说明的是，本节不确定性分析时假定Bennett机构4根连杆长度完全相同，即具有相同的不确定性。在实际应用中，4根连杆长度的不确定性可能不同，而不同的连杆长度可能导致Bennett机构无法正常展开。

4.2.2   连杆弹性模量不确定影响

本节将讨论Bennett机构连杆弹性模量不确定时的动力响应。本节所讨论的Bennett机构含有两组相互独立的不确定参数，即假定AC与BD、BC与DA的弹性模量分别相同。假设两组连杆弹性模量不确定性均为5%，表示为：

式中，区间参数${\xi _1},{\xi _2} \in \left[ { - 1,}\;1 \right]$。利用3阶Chebyshev多项式建立代理模型，所需采样点数量为${4^2} = 16$。

图13绘制了端点B、D间的距离时程曲线。在第一次Bennett机构展开过程中，BD距离上下边界时程曲线基本重合；在第二次Bennett机构展开过程中，上下边界曲线在机构展开部分基本重合；但当机构进入收拢过程中，上下边界差距逐渐增大。这可能是由于发生第二次机构连杆端部截面碰撞后，不同弹性模量的连杆储存的应变能不同，导致机构运动速度不同。总体而言，连杆弹性模量不确定性对BD距离影响相对较小。

图14绘制了端点B处z方向的速度时程曲线。在机构发生连杆截面碰撞前，速度时程曲线上下边界基本重合，此时弹性模量不确定性对机构速度影响较小；在机构连杆截面碰撞后，曲线振荡显著增加，上下边界曲线差异增大。在约0.91 s时刻，速度响应上下边界最大差值为8.42 m/s，占对应时刻确定参数模型速度响应的277.0%。此外，两次碰撞后速度响应上下边界区间都有较大增加。以上结果表明连杆弹性模量的不确定性对速度响应的影响比位移响应更大，而强外部作用会显著增加参数不确定性的影响。

图  13  BD距离时程曲线

Figure  13.  Time history of distance BD

下载: 全尺寸图片 幻灯片

图  14  B点$z$方向速度时程曲线

Figure  14.  Time history of z-direction velocity of point B

下载: 全尺寸图片 幻灯片

图15绘制了梁单元能量时程曲线的上边界和确定参数模型的梁单元能量时程曲线。在机构截面碰撞前，梁单元能量时程曲线及其上边界曲线基本为0；当机构截面碰撞时，梁单元能量突增，随后在一定范围内振荡；在第二次连杆机构截面碰撞后，梁单元能量振幅进一步增大。第一次碰撞前、第一次碰撞后以及第二次碰撞后的梁单元能量上边界最大值与对应时刻梁单元能量比值分别为124.5%、312.7%和258.7%。值得注意的是，在连杆截面第一次碰撞后，梁单元能量上边界曲线与确定参数模型的梁单元曲线间差距显著增加，并且二者差距在第二次截面碰撞后进一步增大。该现象进一步说明碰撞会显著增加机构参数不确定性对动力响应的影响。需要说明的是，代理模型计算出的下边界曲线存在负值，违背梁单元应变能为非负值的原则，这可能是因为不同样本点结果振荡过强，响应值在不同时刻逼近零值附近，使Chebyshev多项式产生错误估计。

图  15  梁单元能量时程曲线

Figure  15.  Time history of beam element

下载: 全尺寸图片 幻灯片

5   结论

本文针对Bennett机构动力响应不确定性分析问题，提出了基于Chebyshev多项式和有限质点法的动力响应不确定性分析方法。给出了该方法的原理和计算流程，通过两个典型机构算例验证了方法的有效性，并开展了Bennett机构动力响应不确定性分析。主要结论如下：

(1) 提出了一种基于Chebyshev多项式和有限质点法的分析区间不确定性方法。对曲柄滑块和双摆机构算例进行不确定性分析，结果表明本文方法计算结果与文献结果基本一致，证明了该方法的有效性。

(2) 利用本文方法开展了Bennett机构动力响应不确定性分析。结果表明连杆长度不确定性对机构位移、速度有较大影响，BD距离响应区间最大差值占对应时刻确定参数模型BD距离的87.1%；连杆弹性模量对位移影响不显著，对速度、梁单元能量影响较大，其中速度响应上下边界最大差值占对应时刻确定参数模型速度的277.0%。

(3) 碰撞等强外部作用会显著增加参数不确定性对机构动力响应的影响。在含构件长度、弹性模量等不确定参数的Bennett机构中，机构截面碰撞前响应上下边界差异均不明显；机构截面碰撞后响应边界差异显著增大。特别是梁单元能量在第一次和第二次碰撞后响应上边界分别增至312.7%和258.7%。

本文所提出的机构动力响应不确定性分析方法不仅适用于Bennett机构，还可推广至其他机构分析中。后续研究将进一步考虑节点间隙大小等不确定性参数对机构动力响应的影响，并将考虑有限质点法与其他不确定性分析方法结合，例如将有限质点法与概率分析方法结合以获得机构动力响应概率分布。