区间过程激励下刚柔耦合系统动态不确定性分析的序列模拟方法

刘延浩 倪冰雨 田万一 姜潮

湖南大學机械与运载工程学院，长沙，410006

摘要：

针对动态不确定性激励作用下的刚柔耦合系统动力学问题，提出了一种基于区间过程模型的动力学不确定性序列模拟方法，旨在通过区间过程序列抽样及刚柔耦合动力学序列模拟计算结构振动与机构运动等系统动态响应的上下边界。介绍了中心刚体柔性梁刚柔耦合系统动力学方程的构建与数值求解方法。针对动态不确定性激励作用下的刚柔耦合系统动力学分析，引入区间过程模型及其区间K-L展开对动态不确定性进行了度量和高效表征，提出了一种求解系统机构运动与结构振动等动力学响应上下边界的序列模拟方法。该方法利用序列模拟策略识别当前模拟序列中对动态响应上边界或下边界具有贡献的区间过程参数样本集，作为下一模拟序列中的局部加密抽样中心，可有效避免直接蒙特卡罗模拟在计算动力学响应上下边界时因过多无效抽样模拟而导致的低效收敛问题。最后，通过三个算例对所提方法的有效性进行了验证。研究结果表明，对于刚柔耦合系统大范围运动及振动响应上下边界的求解，序列模拟方法相比于直接蒙特卡罗模拟方法具有更好的计算效率和计算精度。

关键词：动态不确定性；区间过程模型；区间K-L展开；刚柔耦合动力学；序列模拟方法

中图分类号：TH113

DOI：10.3969/j.issn.1004132X.2024.05.002

开放科学（资源服务）标识码（OSID）：

A Sequential Simulation Method for Dynamic Uncertainty Analysis of

Rigid-flexible Coupling Systems under Interval Process Excitations

LIU Yanhao  NI Bingyu  TIAN Wanyi  JIANG Chao

School of Mechanical and Vehicle Engineering，Hunan University，Changsha，410006

Abstract： For the dynamic problem of rigid-flexible coupling systems under dynamic uncertain excitations， an interval process model-based sequential simulation method was proposed for uncertainty analysis， which aimed to obtain the upper and lower bounds of the system dynamic responses such as structural vibrations and mechanism kinematics， by sequential sampling of the interval process and the rigid-flexible coupling dynamics simulations. The construction and numerical solution of the dynamic equation of the rigid-flexible coupling systems with central rigid body and flexible beam were introduced. Aiming at the dynamic analysis of rigid-flexible coupling systems under uncertain dynamic excitations， the interval process model and the interval K-L expansion were introduced to quantify and represent the dynamic uncertainty efficiently， and a sequential simulation method was proposed to solve the upper and lower bounds of the dynamic responses of the system mechanism motions and structural vibrations. The method used a sequential simulation strategy to identify the interval process parameter sample sets that contributed to the upper or lower bounds of dynamic responses in the cur rent simulation sequence， and served as the local encrypted sampling center in the next simulation sequence， which might effectively avoid the inefficient convergence problem caused by excessive invalid sampling simulations when calculating the upper and lower bounds of dynamic response in direct Monte Carlo simulation. Finally， three examples were given to verify the effectiveness of the proposed method. The results show that the sequential simulation method has better computational efficiency and accuracy than that of the direct Monte Carlo simulation method for solving the upper and lower bounds of the rigid-flexible coupling systems large overall motions and vibration responses.

Key words： dynamic uncertainty； interval process model； interval K-L expansion； rigid-flexible coupling dynamics； sequential simulation method

收稿日期：20240415

基金项目：国家自然科学基金（52175224，52235005，52075156）

0  引言

柔性多体系统动力学研究由刚体和柔性体组成的复杂机械系統在经历大范围空间运动时的动力学行为，在航天器、机器人、兵器与车辆等工程领域得到了广泛关注。刚柔耦合系统的动力学问题是柔性多体系统动力学的本质，主要研究柔性体变形与其大范围空间运动之间的相互耦合作用，以及这种耦合所导致的动力学效应［1-3］。然而，由于服役环境的复杂性和运行工况的多样性，系统通常受到不确定性因素影响导致机构运动、振动等动力学响应产生波动。例如，高速与高精度装配机器人的定位精度必须考虑机械臂在高速运动中与其自身弹性变形之间的刚柔耦合效应，然而基体振动、外部扰动等动态不确定性的影响将降低机器人的装配精度。因此，合理量化上述动态不确定性参数对刚柔耦合动力系统的影响，精确预测系统动力学响应的波动特性，并在设计中加以充分考虑，对提高系统性能及可靠性具有重要意义。

区间过程激励下刚柔耦合系统动态不确定性分析的序列模拟方法——刘延浩  倪冰雨  田万一等

中国机械工程 第35卷 第5期 2024年5月

目前，在刚柔耦合系统动力学不确定性分析领域，国内外学者已经进行了相关研究。吴胜宝等［4］对自由大范围运动情况下刚体柔性梁系统的刚柔耦合动力学问题进行了研究，分析了系统在不同外界激励作用下的耦合动力学响应。为了减少工业机器人的制造和各种环境等随机不确定性因素造成结构系统性能下降，胡启国等［5］采用一次二阶矩法得到臂部可靠性指标，并用可靠性指标反映不确定性因素的影响，并提出了一种考虑刚柔耦合的工业机器人多目标可靠性拓扑优化方法。辛鹏飞等［6］以柔性空间机械臂为例，运用基于Chebyshev多项式的区间扩张函数，将含区间参数的微分代数方程转化为Chebyshev多项式插值点处的确定参数的动力学方程，研究得到了一维区间参数和多维区间参数影响下机械臂系统的动力学响应区间边界，形成了预测机械臂末端轨迹区间的方法。WANG等［7］从不确定性和系统工程的角度研究了火炮系统结构动力响应的不确定优化问题，针对火炮发射过程中内部弹道、发射载荷等参数波动对弹丸运动的影响，提出了一种考虑鲁棒性和经济性的火炮结构动力响应区间不确定优化方法。李勇［8］基于某火炮供药机构刚柔耦合动力学模型，对模块药放入、贮存和回位等多个机构动作的不确定性进行了仿真分析。宋清林［9］考虑运动副间隙和协调臂柔性，研究了多个参数波动对火炮协调器动态特性的不同影响，并进一步分析了协调运动的可靠性。YANG等［10］将刚柔耦合卫星中的多源不确定性量化为区间参数，基于标称刚柔耦合卫星动力学模型建立了相应的耦合状态方程，同时，考虑代理模型中的多项式混沌展开，提出了一种区间分析方法预测不确定姿态振动控制系统的不确定状态响应。根据文献调研，目前针对刚柔耦合动力学系统的不确定性分析主要关注加工误差、材料属性、条件约束等参量不确定性，尚缺乏对动态不确定性因素的考虑。然而，对于工程实际中的诸多刚柔耦合动力学系统，振动激励、交变温度载荷等动态不确定性对此类系统尤其是高速高精度执行机构等系统的动力学性能具有重要影响，相关研究对后续系统高性能保障和高可靠性设计具有重要的理论和应用价值。但不同于刚体动力学系统和结构振动系统，刚柔耦合动力学系统通常同时伴随机构运动和结构振动的复杂耦合作用，使其动力学方程呈现出高度非线性、强耦合和时变性等特点；在随机动态载荷等动态不确定性作用下，对系统机构运动及结构振动的耦合动力学响应进行求解属于强非线性系统的动力学不确定性分析，这对分析方法在计算效率和精度等方面提出了挑战。

本文综合考虑刚柔耦合系统动态不确定性问题的工程特性，提出了一种基于区间过程模型的动力学不确定性序列模拟方法，以获得刚柔耦合系统机构运动状态参数及结构振动响应随时间变化的上下边界。该方法在刚柔耦合系统动态不确定性分析方面的优势主要包括：①对于许多实际工程问题，由于测试条件或经济成本的限制等，通常外部载荷、路面激励等动态不确定参数的测试数据量难以满足随机过程建模的要求，而区间过程由于在任意时刻使用区间而非精确概率分布来描述参数的变化范围，在整个时间历程上通过两条边界曲线描述时变参数的不确定性，可很大程度上降低对参数样本量的依赖；②对于诸多柔性多体系统，所受到的动态不确定性（如摩擦磨损、冲击振动等）通常来源于上游/前端间隙、磨损量、温度等物理属性的综合影响，而此类前端参数的不确定性通常认为属于有界不确定性，因此在复杂系统中演化影响形成的动态不确定性也属于有界动态不确定性，适用于描述为区间过程；③基于区间过程的刚柔耦合动力学不确定性求解可获得结构动态响应边界或机构运动轨迹包络，这种边界或包络形式对机构运动精度分析和机构运动干涉判断十分有利和方便；④刚柔耦合动力学问题通常具有强非线性，并且由于涉及大范围机构运动和结构振动耦合计算，往往计算复杂，因此，为便于工程分析，适宜采用非嵌入式不确定性求解方法。现有考虑区间过程的不确定性分析主要针对结构振动问题，且基本上都是处理线弹性振动系统和少数弱非线性振动系统。本文方法首先引入区间过程模型及区间K-L展开方法对刚柔耦合动力学系统中的动态不确定性进行量化和高效表征，进而利用序列模拟策略识别当前刚柔耦合动力学分析模拟序列中对系统机构运动与结构振动等动力学响应上边界或下边界具有贡献的区间过程参数样本集，并作为下一模拟序列中的局部加密抽样中心，有效避免了直接蒙特卡罗模拟在计算动力学响应上下边界时因过多无效抽样模拟而导致的低效收敛问题，从而有效提高了不确定性求解效率。此外，本文简要介绍了刚柔耦合系统动力学建模与数值求解方法，给出了基于区间K-L展开的刚柔耦合系统动力学不确定性分析序列模拟求解方法，并通过三个算例验证了本文方法的有效性。

1  刚柔耦合系统动力学分析简介

1.1  刚柔耦合系统物理模型

图1所示的中心刚体柔性梁刚柔耦合系统属于一类典型的刚柔耦合问题［4］，其中，中心刚体做平面运动，其上以悬臂方式连接柔性梁。坐标系OXY为系统的惯性坐标系，在柔性梁上建立浮动坐标系oxy，图1中o为刚体与柔性梁的连接处；θ为ox轴与OX轴的夹角；x为未变形时柔性梁的轴线；刚体质量为m，刚体绕质心c的转动惯量为Jch，半径为R；柔性梁的质量为m1，长度为L，密度为ρ，横截面积为S，弹性模量为E，截面惯性矩为I；τ为作用在转轴上的力矩；r为柔性梁变形后一点z到O点的矢径；ω1、ω2分别为梁的轴线纵向伸长量和为横向弯曲变形量；ωc为耦合变形量。

1.2  刚柔耦合系统动力学方程构建与求解

多体系统动力学方程的推导方法有多种，一般常用的有牛顿欧拉法、拉格朗日法和高斯法等［11-12］，在构建多体系统动力学方程中，多采用拉格朗日法，该类方程可直接表示为系统控制输入的函数。对于完整系统用广义坐标表示的动力学方程，通常系指第二类拉格朗日方程。

如图1所示，从惯性坐标系OXY原点引往柔性梁变形后一点z（x，y）的矢径r为［4］

r=rc+ro+ρo+u（1）

式中，rc为刚体质心关于惯性坐标系的矢径；ro为刚体质心至浮动坐标系原点的矢径；ρo为梁变形前z点关于浮动坐标系原点的矢径；u为梁变形前后末端位移矢径。

设λ为浮动坐标系相对于惯性坐标系的方向余弦矩阵，则矢径r在惯性系下的坐标阵为

r=rc+λ（ro+ρo+u）（2）

rc=（X，Y）T  ro=（R，0）T

ρo=（x，y）T  u=（ux，uy）

λ=cos θ－sin θsin θcos θ

将式（2）对时间t求导可得到z点速度为

r·=r·c+λ·（ro+ρo+u）+λu·（3）

系统的动能和势能分别为

Ek=12∫Vρr·Tr·dV+12mr·2c+12Jchθ2（4）

Ep=12∫L0ES（ω1x）2dx+12∫L0EI（2ω2x2）2dx（5）

式中，V为体积。

求解单自由度柔性臂的振动模型可采用假设模态法对模态坐标进行解耦处理，再对结果进行线性变换得到实际的坐标响应。采用假设模态法描述柔性梁的变形，纵向变形ω1和横向变形ω2分别表示为［4］

ω1（x，t）=Φx（x）A（t）

ω2（x，t）=Φy（x）B（t）（6）

Φx（x）=（Φ（1）x（x），Φ（2）x（x），…，Φ（i）x（x），…，Φ（N）x（x））

A（t）=（A1（t），A2（t），…，Ai（t），…，AN（t））T

Φy（x）=）（Φ（1）y（x），Φ（2）y（x），…，Φ（i）y（x），…，Φ（N）y（x））

B（t）=（B1（t），B2（t），…，Bi（t），…，BN（t））T

Φ（i）x（x）=sin（2i－12πxL）  i=1，2，…，N

Φ（i）y（x）=cosh（βix）－cos（βix）－γi（sinh（βix）－

sin（βix））  i=1，2，…，N

γi=cosh（βiL）+cos（βiL）sinh（βiL）+sin（βiL）

式中，Φx（x）、Φy（x）分别为梁的纵向振动和横向振动的模态函数行矢量;A（t）、B（t）分别为纵向振动和横向振动的模态坐标列矢量；N为模态函数的阶数；βi为模态函数行矢量特征值，可由特征方程1+cosh（βiL）cos（βiL）=0求出。

采取广义坐标q=（X，Y，θ，A，B），运用第二类拉格朗日方程可得

ddt（Tq·）－Tq=－Uq+Fq（7）

Fq=（FX，FY，τ，0，0）

式中，Fq为外主动力对应的广义力；FX、FY分别为刚体上合力主矢在X、Y方向上投影；τ为刚体上合力关于质心的矢矩。

将式（4）和式（5）代入式（7），得到刚柔耦合系统动力学方程为

M（q）q¨=Q（q，q·，t） （8）

其中，M、Q分别为广义质量阵和广义力阵，可分别表示为

M=M110M13M14M150M22M23M24M25M31M32M33M34M35M41M42M43M440M51M52M530M55

Q=［QX  QY  Qθ  QA  QB］T（9）

进一步地，可以得到如下方程：

M110M13M14M15

0M22M23M24M25

M31M32M33M34M35

M41M42M43M440

M51M52M530M55X¨Y¨θ¨A¨B¨=QXQYQθQAQB（10）

式（9）、式（10）中部分变量为矩阵，各变量的具体计算表达式可以参考文献［4］。

目前，对上述刚柔耦合系统动力学方程的求解已经可以使用较为成熟的数值方法［13］，主要包括：四阶龙格库塔法（fourth-order Runge-Kutta method，RK4）［14］、自适应变步长的龙格库塔费尔博格法（Runge-Kutta-Fehlberg method，RKF45）；三阶吉尔预估校正法（third-order Gear predictor-

corrector algorithm，3 Gear）；基于纽马克（Newmark）法［15］的希爾伯特修斯泰勒法（Hilber-Hughes-Taylor method，HHT）［16-17］等。上述方法中仅最后2种为隐式积分法，其他均为显式积分法。其中，对于可取较大步长值的仿真计算建议优先选择HHT和3 Gear法等隐式方法，而对于只能取小时间步长的仿真可考虑使用RK4和RKF45法等显式方法。

此外，对于实际工程中的复杂刚柔耦合系统，上述刚柔耦合系统动力学方程的构建方法往往难以应用，因此，对于实际工程中复杂的刚柔耦合问题，通常利用多刚体动力学和有限元法建立系统的刚柔耦合模型，进而获得系统的动力学响应。

2  基于区间K-L展开的刚柔耦合系统动态不确定性序列模拟方法

2.1  区间过程模型及区间K-L展开简介

在传统的随机过程理论中，不确定过程X（t）在任意时刻t被视为一随机变量，任意时刻变量之间的相关性由自相关函数来描述。然而随机过程模型的构建需要大量的样本，在工程实际中很多时候难以获取。而相对于精确的概率分布函数等统计信息，基于有限的样本和工程经验，参数的边界信息往往相对容易获得。如图2所示，在区间过程模型中，任意时刻的不确定参数被描述成一个有界闭区间，不同时间点变量之间的相关性通过相应的相关系数函数来描述［18-20］。

通常将区间过程{X（t）∈XI（t），t∈T}简记为X（t）∈XI（t）或XI（t），其中T为参数t的参数集。关于参数t的函数XL（t）、XU（t）分别为区间过程X（t）∈XI（t）的下边界（函数）和上边界（函数），将区间过程XI（t）的样本函数记为x（t）。区间过程X（t）∈XI（t）的中值函数和半径函数分别定义为

Xm（t）=XU（t）+XL（t）2（11）

Xr（t）=XU（t）－XL（t）2（12）

对于区间过程XI（t），其中值函数Xm（t）反映该不确定过程随时间t的整体变化趋势。半径函数Xr（t）非负，反映了XI（t）的波动程度。半径函数越小，表明该区间过程的不确定性越弱，其中当半径函数Xr（t）=0时，则区间过程XI（t）退化为一确定性函数。XI（t）的协方差函数CovXIXI（ti，tj）和相关系数函数ρXIXI（ti，tj）=CovXIXI（ti，tj）Xr（ti）Xr（tj），ti，tj∈T用于表征区间过程在时间上的相关性。关于区间过程的详细介绍可参考笔者前期研究工作［18-19］。

上述区间过程模型一般不能直接应用于结构不确定分析中，因而应先将其通过数量有限的一组区间变量表示（即区间过程的离散）。基于时间离散的方法为保证精度往往需要大量区间变量来表示时间的连续，但后续计算成本过高。除基于时间离散的方法以外，对于区间过程还可以通过级数展开的形式，将区间过程变换为无穷多个不相关的标准区间变量的确定性时间函数的组合，再保留有限重要项实现对区间过程的近似。

区间过程的区间K-L展开可表征为

X（t）=Xm（t）+∑∞j=1Xr（t）λjφj（t）ζj（13）

式中，ζj∈ζI=［－1，1］（j=1，2，…）为标准区间变量，且满足∑∞i=1ζ2j≤1;λj∈［0，∞）、φj（t）分别为自相关系数函数ρXIXI（t，t′）的特征值与特征函数。

根据Mercer定理，对自相关系数函数ρXIXI（t，t′）有如下谱分解：

ρXIXI（t，t′）=∑∞j=1λjφj（t）φj（t′）（14）

其中，t，t′为任意的两个时刻，特征值λj与特征函数φj（·）为如下第二类齐次方程Fredholm积分方程的解，即

∫TρXIXI（t，t′）φj（t）dt=λjφj（t′）（15）

正交特征函数满足∫Tφi（t）φj（t）dt=δij，其中φi（t）为区别于φj（t）的特征函数，δij为Kronecker-delta函数。在式（15）中，令t=t′，此时有ρXIXI（t，t′）=1，等号两边在T上同时积分可得

∑∞j=1λj=|T|（16）

式中，|T|为时间参数集T的长度。

此外，区间过程的抽样通过区间K-L展开实现。在实际工程中，利用无穷多个区间变量对随时间变化的不确定参数进行度量既不可能也无必要。在实际应用时，可将式（13）所表示的级数项根据特征值的大小进行降序排列，并对前M项进行截断。

记ζ（s）j（j=1，2，…，M）为标准区间变量ζj∈ζIj（j=1，2，…，M）的第s组样本，满足∑Mj=1ζ2j≤1，则区间过程X（t）∈XI（t）的样本函数x（s）（t）可近似表示为

x（s）（t）=Xm（t）+∑Mj=1Xr（t）λjφj（t）ζ（s）j（17）

其中，特征值λj与特征函数φj（·）可根据式（15）确定。由于截断近似带来的误差无法避免，故可将区间过程截断近似的近似度κ定义为

κ=1|T|∑Mj=1λj（18）

其中，截断项数M可根据近似度的要求进行确定；为保证抽样函数的准确性，通常要求κ≥0.95。区间过程X（t）∈XI（t）的样本函数x（s）（t）的具体抽样步骤可参考文献［18］。

2.2  刚柔耦合系统动态响应上下边界的序列模拟求解

在区间过程动态激励作用下，刚柔耦合系统的动力学响应包括结构振动和机构运动特性参数等，同样具有上下边界的动态不确定性［21］，如图3所示。本文所研究的问题主要是对系统动态响应上下边界的求解。根据1.2节中所建立的刚柔耦合动力学方程（式（8）），下面结合区间K-L分解对输入不确定激励进行表征，并给出该方程的非确定性描述形式。

首先，考虑由式（8）所确定的刚柔耦合系统假设受X方向上區间过程激励FX（t）。输入激励FX（t）由区间过程描述，即FX（t）∈FIX（t）=［FLX（t），FUX（t）］，FLX（t）、FUX（t）分别为不确定激励FX（t）的下边界函数和上边界函数，FX（t）的自相关系数函数为ρXIXI（t，t′）。利用区间K-L展开对激励FX（t）进行表示，则有

FX（t）=FmX（t）+∑∞j=1FrX（t）λjφj（t）ζj（19）

将由式（19）所描述的激励代入至刚柔耦合动力学方程（式（8））即可得到完整的非确定性描述形式。同理可得系统受到不同方向不确定激励时的动力学方程，本文不再展开描述。

图4展示出了用于寻找上边界的序列模拟（sequential simulation，SS）方法的流程图。序列模拟是通过一个序列过程来实现的，下边界寻找的过程可通过类似的方式实现。在每个序列中，执行直接蒙特卡罗（Monte Carlo，MC）模拟过程（即下述步骤（1）～步骤（2）），而每个序列中的采样策略需要根据从先前序列获得的结果响应来更新［22］。序列模拟方法的详细步骤如下：

（1）输入区间过程X（t）∈XI（t）及其上下边界函数、自相关系数函数，根据近似度的要求从M维单位球Ωζ={ζ|∑Mj=1ζ2j≤1}内抽样获得标准区间变量ζj∈ζIj=［－1，1］的样本ζ（s）j（j=1，2，…，M）；进而通过区间K-L展开将标准区间变量ζj∈ζIj（j=1，2，…，M）的样本ζ（s）j（j=1，2，…，M）映射为区间过程X（t）∈XI（t）的初始样本函数x（s）（t）（s=1，2，…，ninitial），其中ninitial为初始样本的个数。

（2）将生成的样本输入到所构建的刚柔耦合系统动力学方程或通过多刚体动力学和有限元法所建立的刚柔耦合动力学仿真模型中，求解得到输出响应为{y（x（s）（t），t），s=1，2，…，n}，n=nglobal+nlocal。nglobal、nlocal分别为上一个序列中新生成的全局样本和局部样本的数量。对于初始序列，有nglobal=ninitial，nlocal=0。

（3）從时间参数t的每个时刻{tj，tj∈［0，ts］}中求出所有样本响应{y（x（s）（t），tj），s=1，2，…，n}对应的最大值{y^U（tj），tj∈［0，ts］}，并获得输出响应的上边界y^U（t），其中ts为样本响应的总时间。

（4）识别对上边界y^U（t）的贡献因子，即识别在连续参数时间t上对整个时间历程［0，ts］或部分时间间隔Δt具有贡献最大值的响应实现{y（x（）（t），t）}，并对输入的不确定参数寻找相应的贡献样本{x（*）（t）}。

（5）如果从最近两个序列获得的上界之间的偏差足够小，则认为计算达到收敛并转到步骤（6）；否则生成新的样本（全局样本和局部样本），并转到步骤（2）；全局样本从整个有界空间XI（t）中产生，局部样本在贡献样本的邻域内产生。

（6）将最近一次模拟序列中得到的上边界作为近似解输出。

上述方法为区间过程激励下刚柔耦合系统机构运动与结构振动等动力学响应边界分析提供了一种通用的序列模拟策略。本节的后续内容中，将讨论序列模拟过程中的一些重要问题，如贡献样本的贡献评估、贡献样本邻域的定义、每个序列中局部样本的数量、收敛标准的设置等。

在每个序列模拟中都需要识别出有贡献的样本。如前所述，通常大多数输入样本对输出响应的上边界或下边界没有贡献。例如，在每个序列模拟中获得的上边界必须由已有响应中某些部分组成，因此，这些能够构成响应边界分量的输入样本就称之为贡献样本。如图5所示，将每一个贡献样本x（*）（t）的响应实现记为y（x（）（t），t），其贡献可以通过响应实现覆盖上边界（或下边界）的持续时间来评估，由此每个贡献样本对应的贡献指数定义为

γ［x（）（t）］=Δt［x（）（t）］ts（20）

其中，Δt［x（）（t）］为贡献样本x（）（t）的响应实现覆盖上边界（或下边界）的持续时间，图5中x（1）（t）的响应所覆盖上边界的时间非连续，即由两个时间区域组成，“+”号表示两个区域共同组成了该贡献样本响应所覆盖上边界的总时间Δt［x（1）（t）］。

由序列模拟方法的过程可知，每个序列模拟中都将生成全局样本和局部样本（初始序列中没有局部样本）。全局样本是在输入区间不确定参数的整个有界空间内生成的，而局部样本则来自贡献样本的邻域。全局样本的数目可以根据要解决的数值问题的规模、非线性程度或计算花费来给出。此外，在有贡献的样本附近的局部样本被认为更有可能产生响应边界值，因此，在下一个序列模拟中，优选在贡献样本的邻域内采样。

由上文可知，对于区间过程X（t）的抽样即可等价于对有限个标准区间变量ζj∈ζIj=［－1，1］的抽样；获得区间变量的样本ζ（s）j，借助于式（20）中的转换即可近似获得X（t）的样本函数。同理可知，在贡献样本x（）（t）的邻域内进行抽样亦等价于在对应的区间变量样本ζ（s）j的邻域内进行抽样，进而借助于式（17）中的转换获得局部样本。图6以二维标准区间变量空间为例对不同形状的贡献样本邻域进行介绍。图6a所示的圆形邻域是一个圆形区域，对于多个标准随机变量，邻域变为以贡献样本为中心、以r*为半径的多维球形区域。图6b所示的正方形邻域，它为以贡献样本为中心、边长为2r*的正方形或多维立方体形区域。为了便于描述，本文引入比率β来表示邻域相对于整个有界空间的大小。对于圆形邻域和方形邻域，定义为 r*=βMV（Ωζ），其中V（Ωζ）为有界空间的体积或面积；w（ζI）为区间变量的长度，对于标准区间变量，有w（ζI）=2。通过测试一些区间不确定性问题，

本文建议将比率β设置在1/12到1/4之间。

各个贡献样本邻域内的局部样本数nlocal［x（）（t）］可根据贡献指数γ［x（）（t）］给出，且有nlocal［x（）（t）］∝γ［x（）（t）］。

定义y-（k）（t）为第k次序列模拟所获得的上边界，如图7所示，第k次模拟序列的绝对误差可表示为ε（k）（t）=|y-（k）（t）－y-（k－1）（t）|。则收敛准则可以描述为：

maxt ε（k）（t）≤ε0（21）

式中，ε0为用于控制误差的预设值。

3  算例分析

3.1  振动基柔性臂

考虑图8所示的振动基柔性臂系统，坐标系OXY表示系统的惯性坐标系，在柔性臂上建立浮动坐标系oxy，其中柔性臂可以o点为旋转中心翻转至指定角度，柔性臂的末端与电机轴固接，属于悬臂梁模型约束条件，θ为翻转部分相对于X轴的夹角。x为未变形时柔性臂的轴线。X轴所在平面表示系统未发生垂直振动时所在水平面，O点处为基座质心；Y轴为基座质心的垂直振动的方向。基座的质量取mz=0.5 kg，柔性臂的参数取长度L=8 m，截面积S=7.3×10－3  m2，惯性矩I=8.22×10－9 m4，

密度ρ=2.7667×103 kg/m3，弹性模量E=68.952 GPa。

本算例同时考虑柔性臂结构柔性和基座垂直振动，其中，基座底部受到一动态载荷f（t），且被度量为一中值为0、半径为2 N的平稳区间过程f（t）∈fI（t），其相关性系数函数为ρXIXI（t，t′）=exp（－|t－t′|/lt），其中lt为相关长度，本文取lt=0.3 s。对于指数型自相关系数函数，相关长度lt越大，区间过程f（t）的时间相关性越强，样本函数相对也越平滑［15］。与此同时柔性臂受一确定性力矩驱动下绕o点旋转，且驱动力矩为

τ（t）=τ0sin2πtTm  0≤t≤Tm

0t>Tm（22）

式中，Tm为力矩作用时间，取Tm=10 s；τ0为力矩幅值，取τ0=5 N·m。

在实际中，模态坐标一般由前几阶模态主导，本算例为单自由度柔性臂的振动模型，取前两阶模态即可得到精准的结果［9］。且假定系统只在Y方向做平动，采取广义坐标q=（Y，θ，A，B），其中，A（t）=（A1（t），A2（t））T，B（t）=（B1（t），B2（t））T。

使用第1节所述刚柔耦合系统动力学建模方法，可得到如下方程：

M22M23M24M25

M32M33M34M35

M42M43M440

M52M530M55

Y¨θ¨A¨B¨=QYQθQAQB（23）

本算例中的柔性臂模态取前两阶，则其横向变形ω2（x，t）可描述为如下模态方程：

ω2（x，t）=∑2i=1Φ（i）y（x）Bi（t）=

Φ（1）y（x）B1（t）+Φ（2）y（x）B2（t）（24）

基于上述刚柔耦合方程组，本算例分别通过直接蒙特卡罗模拟方法（以下简称MC法）和序列模拟方法（以下简称SS法）分别计算了柔性臂末端横向位移ω2与柔性臂翻转角θ随时间变化的上下边界，将仿真时间设为t=15 s，且仿真的初始条件均为0。为了保证结果的精度，本文将MC法的模拟次数设置为20 000。在调用SS法对该问题进行求解时，响应上下边界的收敛的最大误差均设置为maxt ε（k）（t）≤ε0=1×10－7，且应用比率β=0.2的多维球形作为邻域。应用MC法与SS法对上述问题进行求解，所得结果如图9和图10所示，结果表明MC法所得边界与SS法所得上下边界高度重合，甚至比MC法所得边界更宽。与基于区间运算的区间不确定性分析方法不同，SS法和MC法均为非嵌入式方法，且给出近似内部解，得到的响应边界将严格地被包络于实际结果之内，因此，计算的响应边界越宽，表明越接近精确响应边界。表1、表2所示为不同时刻下SS法所求边界相比于MC法所求边界的误差，结果表明不同时刻下SS法相对于MC法的误差最大为7.69%。与此同时，相比于MC法中对原模型进行了20 000次调用，SS法在求解上下边界过程中分别对原刚柔耦合方程进行了380次和222次的调用，分別仅有MC法调用次数的1/50和1/100左右。由上述讨论可知，在该刚柔耦合问

题中，SS法能够以较小的计算量获得较高精度的解。

3.2  刚柔耦合运动机构

依托某工程实际刚柔耦合运动机构，建立图11所示刚柔耦合运动机构模型，其中，杆1～杆3和质量块A、B、C、D均为刚体，杆4为柔性体。坐标系OXYZ表示系统的惯性坐标系，在杆4上建立浮动坐标系oxy。杆1和杆2、杆3和杆4分别在B、D处固连；杆2在C处与杆3铰接且可绕质量块C的质心C点旋转，转角θ为杆3转动后与OY轴的夹角。杆1～杆4的长度依次为40 mm、60 mm、70 mm、120 mm，且直径均为6 mm；质量块A的尺寸为10 mm×10 mm×20 mm；其余质量块尺寸均为10 mm×10 mm×10 mm。假定质量块A底部中心受到一不确定性激励f（t），且在其作用下机构产生的时变加速度a可通过底部加速度传感器进行测量，并以区间过程模型进行度量，即a（t）∈aI（t）=［－2，2］m/s2。且杆3绕C点处转动的运动规律为θ=π6sin（20πt）rad。本算例仿真分析选用的材料参数列入表3。

使用UG和ABAQUS软件建立刚柔耦合模型时，杆4是重点关注的部分，作柔性化处理生成柔性体。其余部分作为刚性体建模，然后将刚性和柔性部分装配好，形成刚柔耦合模型，如图12所示。

分别使用蒙特卡罗直接模拟程序和前文所述的序列模拟程序并联合ADAMS仿真软件对系统动力学进行分析，仿真时间t=5 s。为了保证结果的精度，本文将MC法的模拟次数设置为10 000。调用SS方法对该问题进行求解时，响应上下边界的收敛的最大误差设置为ε0=1×10－3 mm，且采用比率β=0.15的多维立方体形邻域；SS法的模拟总数为201次。比较了SS法和MC法评估杆4的末端横向位移ω2响应图的上边界和下边界，仿真结果如图13所示。其中，不同时刻下SS法所求边界相比于MC法所求边界的误差如表4所示。

the end of a flexible rod

如图13和表4所示，虽然SS法模拟次数要少得多，但是SS法仍给出了更好的解。在这种情况下，SS法给出的上限甚至比MC法给出的上限更高，下限比MC法给出的下限更低。如前所述，MC法和SS法均可获得响应近似边界，所得响应界限将严格包含在实际结果内，且表明响应界限越宽，越接近准确结果。由此可以看出，SS法比MC法使用的模拟次数要少得多，且结果更好，故SS法在该数值示例中显示出了更高的计算效率。

3.3  某型码垛机器人执行机构

ABB的某型号码垛机器人的最大载重量为450 kg，作业半径为3.2 m，该本体结构见图14。其结构为典型的双平行四边形结构，主要由基座、腰部、腕部和臂部机构组成［23］，臂部机构按照一定比例组成平行四边形，具有4个相互独立的自

由度，分别为腰部相对于底座的回转运动、大臂的前后仰俯运动、小臂的上下仰俯运动、腕部的回转运动。

码垛机器人臂部结构为悬臂的细长杆件，在极端振动工况下的运动过程中会出现相对较大的柔性变形。由于码垛机器人臂部柔性变形对机器人的运行精度及稳定性影响较大，以刚性件对其进行分析难以满足使用要求，因此采用刚柔耦合方法可更加准确地分析码垛机器人的运行状态。

使用UG和ABAQUS软件建立刚柔耦合模型时，小臂是重点关注的部分，作柔性化处理生成柔性体。其余部分作为刚性体建模，然后将刚性和柔性部分装配好，形成刚柔耦合模型。如图15所示，忽略螺钉、垫片、轴套等对系统影响较小的零件，简化各轴、电机、减速机等，然后设置单位、重力加速度和材料属性等参数。

考虑机构工作时的极端振动工况，假定底座中心受到一垂直不确定性激励f（t）。为模拟此振动工况，本算例以加速度a（t） 的形式对模型进行驱动，且被度量为一中值为0、半径为20 m/s2的平稳区间过程，其相关性系数函数为ρXIXI（t，t′）=exp（－|t－t′|/lt），取lt=0.4 s。机构一个工作周期的工作轨迹如图16所示：从0～0.5 s，机器人末端执行器向下移动700 mm，然后停顿0.5 s抓取400 kg的载荷；从1.0～1.5 s，末端执行器向上抬高700 mm；然后腰部旋转90°同时伴随手臂运动达到码放位置，这些动作在1.5～3.0 s之间完成；之后用0.5 s码放负载；最后用0.5 s返回到初始位置。完成一个完整的作业总共用时4 s。

分别使用MC法和上文所述的SS法并联合ADAMS仿真软件对上述仿真条件下机器人末端执行器的速度、加速度曲线进行分析，仿真时间t=4 s。

基于上述刚柔耦合系统仿真分析，本算例分别通过MC法和SS法依次计算了机器人末端执行器的速度、加速度随时间变化的上下边界。为了保证结果的精度，本文将MC法的模拟次数设置为10 000。调用SS法对该问题进行求解时，响应上下边界的收敛的最大误差均设置为ε0=1×10－4，且应用比率β=0.15的多维球形作为邻域。应用MC法和SS法对上述问题进行求解，所得结果如图17和图18所示，结果表明SS法所得边界与MC法所得上下边界高度重合。表5和表6所示为不同时刻下SS法所求边界相比于MC法所求边界的误差，结果表明两者最大误差为9.76%。与此同时，相比于MC法中对系统仿真模型进行了10 000次调用，SS法在求解上下边界过程中分别对原刚柔耦合方程进行了333次和376次的调用，分别仅有MC法调用次数的3%与5%左右。

如图17、图18和表5、表6所示：在整个工作周期内，加速度波动较大，说明机器人在运动过程中存在较大振动，尤其是在运动状态发生改变时，运动极不平稳，从而导致机器人运动精度降低。且在不同时刻下SS法所得边界与MC法模拟

10 000次后的边界高度重合甚至更宽。如前所述，这两种方法计算所得均為近似内部解，得到的响应边界将严格地被包络于实际结果之内，计算所得响应边界越宽，表明越接近精确响应边界，因此在较为复杂的刚柔耦合问题中，序列模拟方法仍能够以较小的计算量获得较高精度的解。

4  结论

本文针对动态不确定性激励作用下的刚柔耦合系统动力学问题提出了一种基于区间过程模型的动力学不确定性序列模拟方法，可高效计算结构振动与机构运动耦合动力学响应的上下边界。考虑到直接蒙特卡罗模拟方法在进行系统动力学响应上下边界时，因绝大多数区间过程抽样及基于此的刚柔耦合动力学分析模拟对机构运动状态参数和结构振动响应的边界求解不具有任何作用，导致计算效率极低的问题，本文方法利用序列抽样模拟的方式逐步逼近系统在整个时域上的动态响应上下边界。该方法通过识别当前刚柔耦合

动力学计算模拟序列中对动态响应上边界或下边界具有贡献的区间过程参数样本集，作为下一模拟序列中的局部加密抽样中心，从而有效提高了响应边界的计算收敛性。算例验证结果表明，对于刚柔耦合系统大范围运动及振动响应上下边界的求解，相比于直接蒙特卡罗模拟方法，序列模拟方法能以较小的计算量获得较高精度的解。

参考文献：

［1］  陈思佳，刚柔耦合问题与空间多杆柔性机械臂的动力学建模理论研究［D］. 南京：南京理工大学，2012.

CHEN Sijia.Research on Dynamic Modeling Theory of Rigid-flexible Coupling Problem and Spatial Multi-rod Flexible Manipulator［D］. Nanjing：Nanjing University of Science and Technology， 2012.

［2］  王斌锐， 方水光， 严冬明. 机器人手臂的刚柔耦合建模及摆动模态对比［J］.中国机械工程，2012，23（17）：2092-2097.

WANG Binrui， FANG Shuiguang， YAN Dongming. Rigid-flex Coupling Modeling and Swing Mode Comparison of Robot Arm［J］.China Mechanical Engineering，2012，23（17）：2092-2097.

［3］  庄未， 方琛玮， 刘晓平. 1P5R柔性关节机械臂动力学分析与控制［J］.中国机械工程，2009，20（24）：2907-2911.

ZHUANG Wei， FANG Chenwei， LIU Xiaoping. Dynamic Analysis and Control of 1P5R Flexible Joint Manipulator［J］.China Mechanical Engineering， 2009， 20（24）：2907-2911.

［4］  吴胜宝，章定国. 大范围运动刚体-柔性梁刚柔耦合动力学分析［J］. 振动工程学报，2011，24（1）：1-7.

WU Shengbao， ZHANG Dingguo. Dynamic Analysis of Rigid-flexible Coupling of Large-scale Motion Rigid Body-flexible Beam［J］. Journal of Vibration Engineering， 2011， 24（1）：1-7.

［5］  胡启国，周松，考虑刚柔耦合的工业机器人多目标可靠性拓扑优化［J］. 计算机集成制造系统，2020，26（3）：623-631.

HU Qiguo， ZHOU Song. Multi-objective Reliability Topology Optimization of Industrial Robots Considering Rigid-flexible Coupling［J］. Computer Integrated Manufacturing Systems， 2020，26（3）：623-631.

［6］  辛鹏飞，荣吉利，项阳，等. 柔性空间机械臂区间参数不确定性分析［J］. 北京理工大学学报，2017，37（10）：1056-1060.

XIN Pengfei， RONG Jili， XIANG Yang， et al. Uncertainty Analysis of Interval Parameters of Flexible Space Manipulator［J］. Journal of Beijing Institute of Technology， 2017，37（10）：1056-1060.

［7］  WANG L， YANG G， XIAO H， et al. Interval Optimization for Structural Dynamic Responses of an Artillery System under Uncertainty［J］. Engineering Optimization， 2020， 52（2）：343-366.

［8］  李勇. 某供輸弹机供药机构的可靠性研究［D］.南京：南京理工大学，2017.

LI Yong. Research on the Reliability of the Medicine Supply Mechanism of a Certain Bomb Feeding Machine［D］. Nanjing：Nanjing University of Science and Technology， 2017.

［9］  宋清林.考虑间隙和柔性体的火炮协调器动力学与可靠性分析［D］.沈阳：东北大学，2018.

SONG Qinglin. Dynamics and Reliability Analysis of Artillery Coordinator Considering Gaps and Flexible Bodies［D］. Shenyang：Northeastern University， 2018.

［10］  YANG C， LU W， XIA Y. Reliability-constrained Optimal Attitude-vibration Control for Rigid-flexible Coupling Satellite Using Interval Dimension-wise Analysis［J］. Reliability Engineering & System Safety， 2023， 237：109382.

［11］  郝秀清，胡福生，陈建涛.基于牛顿欧拉法的3PTT并联机构动力学分析及仿真［J］.中国机械工程，2006（增刊2）：32-36.

HAO Xiuqing， HU Fusheng， CHEN Jiantao. Dynamic Analysis and Simulation of 3PTT Parallel Mechanism Based on Newton-Euler Method［J］.China Mechanical Engineering，2006（S2）：32-36.

［12］  朱畅. 基于奇异摄动法的振动基柔性机械臂的非线性鲁棒控制研究［D］. 武汉：武汉科技大学，2018.

ZHU Chang. Research on Nonlinear Robust Control of Vibration-based Flexible Manipulator Based on Singular Perturbation Method［D］. Wuhan：Wuhan University of Science and Technology， 2018.

［13］  郭晛. 受约束柔性多体系统动力学数值计算方法若干问题研究［D］. 南京：南京理工大学，2019.

GUO Xian. Research on Several Issues in Numerical Calculation Methods for Constrained Flexible Multi-body System Dynamics［D］. Nanjing：Nanjing University of Science and Technology， 2019.

［14］  袁兆鼎，费景高，刘德贵. 刚性微分方程初值问题数值解法［M］. 北京：科学出版社，1987.

YUAN Zhaoding， FEI Jinggao， LIU Degui. Numerical Solution to Initial Value Problem of Rigid Differential Equations［M］. Beijing：Science Press， 1987.

［15］  NEWMARK N M. A Method of Computation for Structural Dynamics［J］. Journal of the Engineering Mechanics Division，1959，85（3）：67-94.

［16］  HILBER H， HUGHES T， TAYLOR R. Improved Numerical Dissipation for Time Integration Algorithms in Structural Dynamics［ J ］. Earthquake Engineering & Structural Dynamics，1977，5（2）：283-292.

［17］  NEGRUT D， RAMPALLI R， OTTARSSON G， et al. On an Implementation of the Hilber-Hughes-Taylor Method in the Context of Index 3 Differential-algebraic Equations of Multibody Dynamics［J］. Journal of Computational and Nonlinear Dynamics， 2007，2（1）：73-85.

［18］  倪冰雨. 区间过程与区间场模型及在结构不确定性分析中的应用［D］.长沙：湖南大学，2018.

NI Bingyu. Interval Process and Interval Field Models and Their Application in Structural Uncertainty Analysis［D］. Changsha：Hunan University， 2018.

［19］  JIANG C， NI B， HAN X， et al．Non-probabilistic Convex Model Process：a New Method of Time-variant Uncertainty Analysis and Its Application to Structural Dynamic Reliability Problems， Computer Methods in Applied Mechanics and Engineering， 268 （2014） 656-676.

［20］  JIANG C， HAN X， LU G Y， et al． Correlation Analysis of Non-probabilistic Convex Model and Corresponding Structural Reliability Technique， Computer Methods in Applied Mechanics and Engineering， 2011，200：2528-2546.

［21］  曾光， 段民封， 倪冰雨， 等. 机械臂架系统的非随机振动分析［J］. 中国机械工程， 2019， 30（10）：1142.

ZENG Guang， DUAN Minfeng， NI Bingyu， et al. Non-random Vibration Analysis of Robotic Boom System［J］. China Mechanical Engineering， 2019， 30（10）：1142.

［22］  NI B Y， JIANG C， WU P G， et al. A Sequential Simulation Strategy for Response Bounds Analysis of Structures with Interval Uncertainties［J］. Computers & Structures， 2022， 266：106785.

［23］  陳亚梅. 基于ADAMS的码垛机器人动力学刚柔耦合分析［J］. 包装工程，2021，42（17）：261-265.

CHEN Yamei. Rigid-flexible Coupling Analysis of Palletizing Robot Dynamics Based on ADAMS［J］. Packaging Engineering， 2021， 42（17）：261-265.

（编辑  胡佳慧）

作者简介：

刘延浩，男，1999年生，硕士研究生。研究方向为机械结构动力学分析。E-mail：lyh1@hnu.edu.cn。

倪冰雨（通信作者），男，1989年生，副教授。研究方向为结构不确定性分析。E-mail：nby@hnu.edu.cn。