CN111159636A - 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法 - Google Patents

Abstract

本发明涉及动力学系统优化技术领域，提出一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法。首先，基于绝对节点坐标方法，建立柔性多体系统的质量矩阵、刚度矩阵和广义力列阵；其次，建立柔性多体系统的动力学方程和优化目标函数；再次，基于直接微分法或伴随变量法，建立柔性多体系统动力学的半解析灵敏度计算公式；最后，求解柔性多体系统动力学微分代数方程，获得灵敏度计算结果。本发明根据绝对节点坐标方法和多体系统动力学理论，建立柔性多体系统的半解析灵敏度计算公式，以解决柔性多体系统动力学的灵敏度分析问题，目的在于提供一套柔性多体系统灵敏度分析的新策略，为大规模复杂的柔性多体系统的灵敏度计算提供便利。

Description

基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度 分析方法

技术领域

本发明属于多体系统动力学优化技术领域，涉及一种柔性多体系统动力学的灵敏度分析方法，尤其涉及一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法。

背景技术

近年来，多体系统动力学及其优化分析在航空航天、机械、汽车等领域发挥着越来越重要的作用。如果使用基于梯度的优化算法进行多体系统动力学的优化设计，则需要对多体系统进行灵敏度分析。灵敏度分析不仅可以用于确定部分优化算法的迭代方向，也可以表征设计变量对目标函数的影响程度，从而可以通过降低设计变量的个数来提高优化效率。因此，多体系统动力学的灵敏度分析已成为多体系统优化设计的核心问题。

目前，对多体系统灵敏度分析的研究大多是针对刚性多体系统，由于柔性多体系统的高维和强非线性，使其求解相对困难，因此，对柔性多体系统灵敏度分析的研究相对较少。柔性多体系统的建模方法主要有增量有限元法、浮动坐标方法、大转动矢量方法和绝对节点坐标方法。相比于其它三种方法，绝对节点坐标方法有效地避开了柔性体有限转动的参数化问题，在全局坐标系下采用插值函数来描述大变形和大范围转动，能精确地描述柔性多体系统高度的几何非线性。

相比于传统静态优化设计，多体系统动力学动态优化设计的目标函数和约束方程包含状态变量，因此，对其进行优化设计的灵敏度包括状态灵敏度(状态变量对设计变量的灵敏度) 和设计灵敏度(目标函数对设计变量的灵敏度)两部分。并且，优化设计还受多体系统状态方程的约束，其状态方程往往表现为一组高维和强非线性的微分代数方程组，给灵敏度的求解带来了困难。目前，对多体系统动力学灵敏度分析的主要方法有：有限差分法、直接微分法和伴随变量法。

有限差分法是一种近似计算方法。它只需要对设计变量做扰动，然后采用差商的方式计算目标函数对设计变量的灵敏度。但其计算的工作量会随设计变量的数量成比例增加，计算效率和精度相对较低。1991年，Greene等采用有限差分法计算了线性、结构和瞬态响应问题的灵敏度，发现这种方法在实际计算时，对计算程序时间和精度的要求很高。直接微分法是通过直接将多体系统的动力学方程和约束方程对设计变量求导，计算出状态灵敏度，从而求解多体系统的灵敏度。该方法由Krishnaswami和Bhatti在1984年首先提出，随后，Chang 等人研究了采用直接微分法计算约束动力系统设计灵敏度系数矩阵的一般方法，并通过两个算例说明了该方法的有效性；Dias等人采用直接微分法建立了刚-柔多体系统灵敏度分析的灵敏度方程，并通过两个刚-柔多体系统的例子，将计算结果与有限差分法进行了比较；Neto等人将直接微分法应用于具有复合材料的柔性多体系统的灵敏度求解。伴随变量法是通过引入一系列伴随变量，消掉方程中包含状态灵敏度的相关项，从而求解得到的一系列伴随变量方程来计算灵敏度。该方法于1978年由Haug和Arora首次提出，因其具有计算速度较快的特点，在近几年应用较为广泛。Liu采用伴随变量法推导出了约束柔性多体系统的一阶和二阶灵敏度分析方程；Zhang和Chen采用伴随变量法，将以大变形问题梁板建模为重点的绝对节点坐标公式扩展到柔性多体系统的设计灵敏度分析中；Alexander等人采用浮动坐标方法进行柔体系统建模，将伴随变量法应用于具有运动回路的柔性多体系统。

然而，上述提到的直接微分法和伴随变量法是采用解析方法来计算设计变量的灵敏度，当多体动力学系统规模很大、结构较复杂时，采用解析方法进行灵敏度计算所需推导的动力学公式将会十分繁杂，甚至可能会出现无法解析求得多体系统对某些设计变量的灵敏度的情况。采用半解析灵敏度分析方法，可以为大规模复杂的柔性多体系统的灵敏度计算提供便利。

发明内容

本发明提出一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法。该方法基于绝对节点坐标方法，建立柔性多体系统的动力学模型，以解决柔性多体系统的灵敏度分析问题，目的在于提供一种柔性多体系统灵敏度分析的新策略，为大规模复杂的柔性多体系统的灵敏度计算提供便利。

为了达到上述目的，本发明采用的技术方案为：

一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，步骤如下：

一、基于绝对节点坐标方法，建立柔性多体系统单元的质量矩阵、刚度矩阵和广义力列阵

绝对节点坐标方法在全局坐标系下描述柔性多体系统的运动形态，将每个节点的广义坐标描述为位置矢量坐标和该点的斜率矢量坐标。本发明基于三维建模的一维二节点梁单元：如图2所示，每个节点有6个坐标，分别为3个位置矢量坐标和3个斜率矢量坐标，每个单元包含2个节点。假设柔性多体系统节点的个数为n。

所述梁单元坐标表示为e＝[e_u e_v]^T，其中，e_p(p＝u,v)表示节点p的坐标：

式中，r_pX,r_pY,r_pZ(p＝u,v)分别表示节点u,v的位置矢量在全局坐标系中的分量；

分别表示节点u,v的斜率坐标。

因此，梁单元上任意点的位置可以表示为：

式中，

表示定义在总体坐标系上的单元形函数矩阵，写为：

S＝[S₁I,S₂I,S₃I,S₄I] (3)

其中，

为单位矩阵；形函数S₁＝1-3ξ²+2ξ³，S₂＝L(ξ-2ξ²+ξ³)，S₃＝3ξ²-2ξ³，S₄＝L(ξ³-ξ²)，其中，变量ξ＝x/L,x∈[0,L]。L为单元的初始长度，x为单元的局部坐标。

1、建立柔性多体系统单元的质量矩阵

相比于其它的柔性多体系统建模方法，采用绝对节点坐标方法推导的系统质量矩阵是常数矩阵。一般地，通过单元的动能计算公式

可以得到单元的质量矩阵

式中，ρ为单元的密度；A为单元的横截面积；V为单元的体积。

2、建立柔性多体系统单元的广义力列阵

柔性多体系统的广义力列阵包括广义外力列阵和广义弹性力列阵。其中，广义外力列阵包括重力和接触力等，广义弹性力列阵包括纵向拉伸变形和横向弯曲变形。

假设系统只受到重力作用，通过虚功原理，可以得到重力虚功：

式中，G为单元的重力矢量，δr为单元的虚位移。单元的广义重力

表示为：

其中，g为重力加速度。

假设梁单元为各向同性，单元的总应变能U包括与纵向拉伸变形相关的应变能U_l和与横向弯曲变形相关的应变能U_t：

式中，ε为单元的纵向拉伸应变；κ为单元在当前构型下的曲率；E为单元的弹性模量； I为单元的截面惯性矩。

单元广义弹性力

表示为单元应变能对单元坐标的导数：

式中，下标x表示单元形函数对物质坐标x的导数；双下标xx表示单元形函数对物质坐标x的二阶导数。单元的纵向拉伸应变ε的具体表达式为：

其中，l为单元的初始纵向长度；l_s为变形后单元的纵向长度，可以通过对微元弧长积分得到：

3、建立柔性多体系统单元的刚度矩阵

由柔性多体系统单元的广义力列阵，可以得到系统单元的刚度矩阵

二、建立柔性多体系统的动力学方程和优化目标函数

本发明研究的问题是基于微分代数方程形式的柔性多体系统动力学数学模型，系统的动力学方程表示为：

式中，

为系统的广义质量矩阵；Φ＝[Φ₁,Φ₂,…,Φ_m]^T为系统的位置约束列阵；

为位置约束的Jacobian矩阵；

为动力学方程中的Lagrange乘子矢量；q、

和

分别为系统的广义坐标、广义速度和广义加速度；

为系统的广义外力列阵；

为系统的广义弹性力列阵，其依赖于系统的设计变量。

目标函数是评价优化设计方案的标准，在多体系统动力学优化问题中，其一般表示为如下的积分形式：

式中，前两部分G⁰和G^f分别与系统的初态和终态有关，第三部分H为被积分项，其与系统的中间过程有关。其中，上标0和f分别表示相关参量的初始和终止时刻值； q＝[q₁,q₂,…,q_6n]^T为状态变量，描述系统的动态响应；b＝[b₁,b₂,…,b_c]^T为设计变量，是优化设计过程中计算选择并需要最终确定的待求参数，其下标c表示设计变量的数目；t⁰和t^f分别表示初始时刻和终止时刻，也可以表示某些状态变量或其速度达到某一定值的特定时刻，一般地，其可以由下式确定：

式中，Ω⁰和Ω^f分别表示初始和终止时刻条件。

系统的初始状态依赖于设计变量，其应满足以下相容附加条件：

式中，

和

分别表示初始状态位置和速度的相容附加条件，其需要使得

和

满秩。

三、基于直接微分法或伴随变量法，建立柔性多体系统动力学的半解析灵敏度计算公式

本发明提出的柔性多体系统动力学的半解析灵敏度分析方法，相较于解析方法，其不需要解析地推导优化设计目标函数和柔性多体系统动力学方程对设计变量的导数，而是将解析求导公式中对设计变量求导的部分用局部有限差分代替。传统半解析灵敏度分析方法是基于单元级别进行灵敏度计算，而改进半解析灵敏度分析方法是基于总体级别，其无需提取扰动前后柔性多体系统单元的相关矩阵信息，使程序实现过程更为简便，具有较高的计算效率。

(a)基于直接微分法建立柔性多体系统动力学的半解析灵敏度计算公式

1、建立柔性多体系统的半解析灵敏度计算公式

灵敏度是目标函数对设计变量的偏导数，一般形式的柔性多体系统灵敏度计算的解析公式可以表示为：

式中，下标表示对相应变量的导数。其中，q_b、

和λ_b为状态灵敏度，分别表示为状态变量相关量q、

qⁱ和λ_b(i＝0,f)对设计变量的导数；

和

分别为系统起始和终止时刻对设计变量的导数，由系统的起始和终止时刻条件确定，将约束条件式(14)对设计变量求导，可得到：

本发明采用改进半解析灵敏度分析方法，将上述解析方法中对设计变量的求导项分别进行总体级别的局部有限差分，可得到改进半解析灵敏度分析方法的灵敏度计算公式：

其中，方程右边的

和

为总体级别的局部有限差分项，其分别表示为：

式中，Δb表示对设计变量的扰动。

2、建立柔性多体系统的状态灵敏度计算公式

采用直接微分法求解柔性多体系统的灵敏度，需要求解系统的状态灵敏度。将柔性多体系统的动力学方程对设计变量求导，可得：

为保证计算结果的准确性，采用约束违约自动稳定方法求解，此时，系统的约束方程同时包含位置约束、速度级约束及加速度级约束：

式中，

为系统的速度级约束；

为系统的加速度级约束；η₁和η₂为自定义参数。

将式(25)对设计变量求导，并与式(24)中的第一个式子联立新的矩阵微分代数方程组：

其中

式中，为了与位置约束Φ进行区分，将

定义为Baumgarte约束。

采用解析方法求解上述矩阵微分代数方程组，可解析地求得柔性多体系统的状态灵敏度

q_b、

和λ_b。

本发明采用改进半解析灵敏度分析方法，将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分，可得到改进半解析灵敏度分析方法的状态灵敏度计算公式：

其中，方程右边的

和

为总体级别的局部有限差分项，其分别表示为：

式中，M、Φ_q、Q_s、Q_k、

和λ都是基于总体级别。

传统半解析灵敏度分析方法，是将上述解析公式中对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

式中，M_e为单元级别的质量矩阵；Φ_qe为单元级别的位置约束Jacobian矩阵；Q_se为单元级别的广义外力列阵；Q_ke为单元级别的广义弹性力列阵；

为单元级别的Baumgarte约束；

和λ_e分别为单元加速度矢量和Lagrange乘子矢量。

3、建立柔性多体系统的状态灵敏度初值计算公式

求解上述状态灵敏度计算公式，必须要先给定状态灵敏度的初值

和

将系统初始状态的相容附加条件对设计变量求导，可得：

将

入上式，并与系统初始状态的位置约束Φ⁰和速度级约束

关于设计变量的导数方程联立，可得：

采用解析方法求解上述矩阵微分方程组，可解析地求得柔性多体系统的状态灵敏度初值

和

本发明采用改进半解析灵敏度分析方法，将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分，可得到改进半解析灵敏度分析方法的状态灵敏度初值计算公式：

其中，方程右边的

和

为总体级别的局部有限差分项，其分别表示为

式中，Φ⁰、

和

都是基于总体级别。

传统半解析灵敏度分析方法，是将上述解析公式中对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

式中，

和

分别为系统初始状态单元级别的位置约束和速度级约束；

和

为系统单元级别的初始相容附加条件。

(b)基于伴随变量法建立柔性多体系统动力学的半解析灵敏度计算公式

1、引进伴随变量

首先，引进伴随变量

和

将其分别转置左乘柔性多体系统的动力学方程，并在t⁰到t^f上积分,可得：

将式(51)分别关于

和

依次进行分部积分，可以得到：

接着，引进伴随变量

和

将其分别转置左乘柔性多体系统的位置约束、速度级约束及起始条件对设计变量的导数方程，可以得到：

2、建立伴随变量方程

将柔性多体系统的灵敏度计算公式分别关于

和

依次进行分部积分，可以得到：

将所有左乘伴随变量的方程式(52)-(58)分别累减到上述两次分部积分后的柔性多体系统的灵敏度计算公式中，令状态灵敏度

q_b、

和λ_b的同类项系数为零，即可得到伴随变量方程：

通过求解上述一系列伴随变量方程，即可求得相应的伴随变量。

3、建立柔性多体系统的半解析灵敏度计算公式

消掉状态灵敏度

q_b、

和λ_b后的柔性多体系统的灵敏度计算公式，与伴随变量μ、ν、σ⁰、ρ⁰、ξ⁰、σ^f、ρ^f、ξ^f、π和β有关，表示为：

采用解析方法求解上述微分方程，可解析地求得柔性多体系统的灵敏度。

本发明采用改进半解析灵敏度分析方法，将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分，可得到改进半解析灵敏度分析方法的灵敏度计算公式：

其中，方程右边的

和

为总体级别的局部有限差分项，其分别表示为

式中，M、Φ_q、Q_s、Q_k、Φ、Φⁱ、

和λ(i＝0,f)都是基于总体级别；μ、ν、π、β、σⁱ、ρⁱ和ξⁱ(i＝0,f)是伴随变量。

传统半解析灵敏度分析方法，是将上述解析公式中对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

式中，Φ_e为单元级别的位置约束；M_e、Φ_qe、Q_se、Q_ke、

和λ_e(i＝0,f) 都是基于单元级别。

四、求解柔性多体系统动力学微分代数方程组，获得灵敏度计算结果

基于绝对节点坐标描述的含约束柔性多体系统的灵敏度计算公式，是指标3的微分代数方程组(DAEs)，其计算方法分为显示方法和隐式方法。显示方法不需要平衡迭代，计算速度较快，但该方法是条件稳定的，其计算结果受积分步长影响较大，收敛性差；隐式方法需要在积分步内迭代求解，计算速度较慢，但该方法是无条件稳定的，可以选取较大的积分步长进行计算。针对本发明的具有大变形、大位移的柔性多体系统动力学问题，通常采用隐式方法进行求解。

本发明采用广义α算法进行求解。考虑公式的简洁性，令Q＝Q_s-Q_k，则离散DAEs的时间区间后，可得到如下非线性代数方程组：

式中，h表示积分时间步长；上标“(_d)”表示第d个迭代步；ζ和γ为算法参数；α为引入的算法辅助基矢量，其满足：

其中，a_m和a_f为算法参数。

为了保证计算结果的精度和数值稳定性，算法参数a_m、a_f、ζ和γ分别取为：

式中，ρ_r∈[0,1]是谱半径，其取值大小决定算法能量耗散的频率范围，本发明取ρ_r＝0.8。令辅助基矢量和Lagrange乘子矢量组成的列向量为未知变量，采用Newton–Raphson法迭代求解上述非线性代数方程组。随时间步进，即可得到仿真时间内的柔性多体系统的灵敏度计算结果。其中，采用Newton–Raphson迭代法求解时，涉及到的雅克比矩阵J为：

其中

至此，即完成了本发明提出的基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法。

本发明的有益积极效果：

1.本发明提出的柔性多体系统动力学的半解析灵敏度分析方法，相较于有限差分法，其不需要反复地求解柔性多体系统动力学的微分代数方程组，具有更高的计算效率和计算精度；

2.本发明提出的柔性多体系统动力学的半解析灵敏度分析方法，相较于解析方法，其不需要解析地推导优化设计目标函数和柔性多体系统动力学方程对设计变量的导数，而是将解析求导公式中对设计变量求导的部分进行总体级别的局部有限差分。因此，它可以灵活地处理各种类型的设计变量，并且克服了部分多体系统动力学方程对设计变量解析导数求解困难的问题，具有更强的通用性；

3.本发明提出的柔性多体系统动力学的半解析灵敏度分析方法，将基于单元级别的传统半解析灵敏度计算方法转换为基于总体级别的改进半解析灵敏度计算方法，无需提取扰动前后柔性多体系统单元的相关矩阵信息，避免了相近数值加减造成有效数字损失的精度问题，并且使程序实现过程更为简便，具有较高的计算效率。

附图说明

图1是本发明实施流程图。

图2是基于三维建模的一维二节点梁单元。

图3是柔性绳索系统动力学模型。

图4是柔性绳索系统在1s时间范围内的仿真结果。

图5是解析法、改进半解析法和有限差分法的灵敏度计算结果随扰动步长变化曲线；其中，图(a)为目标函数对柔性绳索长度的灵敏度；图(b)为目标函数对柔性绳索横截面积的灵敏度；图(c)为目标函数对柔性绳索密度的灵敏度；图(d)为目标函数对柔性绳索弹性模量的灵敏度。

图6是解析法、改进半解析法和有限差分法的灵敏度计算时间历程曲线；其中，图(a)为目标函数对柔性绳索长度的灵敏度；图(b)为目标函数对柔性绳索横截面积的灵敏度；图(c)为目标函数对柔性绳索密度的灵敏度；图(d)为目标函数对柔性绳索弹性模量的灵敏度。

图7是解析法和改进半解析法的部分状态灵敏度时间历程曲线；其中，图(a)为绳索末端点x方向的位置坐标对柔性绳索长度的灵敏度；图(b)为绳索末端点x方向的位置坐标对柔性绳索横截面积的灵敏度；图(c)为绳索末端点x方向的位置坐标对柔性绳索密度的灵敏度；图 (d)为绳索末端点x方向的位置坐标对柔性绳索弹性模量的灵敏度。

具体实施方式

结合图1，本发明的具体实施方式如下：

实施方式基于描述柔性多体系统的绝对节点坐标建模方法和柔性多体系统动力学理论，可用于柔性多体系统动力学灵敏度问题的分析，具体按以下步骤实现：

一、基于绝对节点坐标方法，推导柔性多体系统单元的质量矩阵M_e、刚度矩阵K_e及广义力列阵Q_ge和Q_ke。

二、基于柔性多体系统动力学理论，在步骤1的基础上，建立柔性多体系统的动力学方程。具体步骤为：

(1)给出该柔性多体系统的尺寸参数、划分单元的数目和单元材料信息；

(2)给出该柔性多体系统的初始状态信息，包括初始位置、初始速度等；

(3)给出该柔性多体受到的重力等外力信息；

(4)建立该柔性多体系统总体的质量矩阵、刚度矩阵和广义力列阵；

(5)给出该柔性多体系统的约束方程；

(6)在步骤(1)-(5)的基础上，便可建立柔性多体系统的动力学方程，即微分代数方程组 (DAEs)。

三、给出柔性多体系统优化的设计变量，以及优化的初态和终态约束条件，建立系统的优化目标函数。

四、基于直接微分法或伴随变量法，采用改进半解析方法，求解柔性多体系统的灵敏度。

(a)基于直接微分法求解柔性多体系统的灵敏度。具体步骤为：

(1)将系统的目标函数对设计变量求导，并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分，得到灵敏度计算公式Ψ_b1；

(2)将系统的起始和终止时刻条件对设计变量求导，并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分，得到

(3)将系统初始状态的相容附加条件对设计变量求导，将

代入，并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分；

(4)将系统的位置约束和速度级约束对设计变量求导，并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分；

(5)联立步骤(4)和步骤(5)的矩阵微分方程组，解得系统的状态灵敏度初值

和

(6)采用约束违约自动稳定方法，重新建立系统的约束方程；

(7)将系统的动力学方程对设计变量求导，并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分；

(8)利用步骤(5)的状态灵敏度初值，采用广义α算法求解步骤(7)的矩阵微分方程组，得到系统的状态灵敏度

q_b、

和λ_b；

(9)将步骤(8)得到的系统的状态灵敏度，代入步骤(1)中的灵敏度计算公式，即可解得柔性多体系统的灵敏度。

(b)基于伴随变量法求解柔性多体系统的灵敏度。具体步骤为：

(1)将系统的目标函数对设计变量求导，并对

和

依次进行分部积分，得到系统的灵敏度计算公式Ψ_b2；

(2)引入伴随变量μ和ν，将它们分别转置左乘系统的动力学方程，并在t⁰到t^f上积分；

(3)将步骤(2)得到的第一个方程对

和

依次进行分部积分；

(4)引进伴随变量σⁱ、ρⁱ、ξⁱ、π和β(i＝0,f)将它们分别转置左乘系统的位置约束、速度级约束和起始条件对设计变量的导数方程；

(5)将步骤(2)得到的第二个方程、步骤(3)和步骤(4)得到的方程，分别累减到步骤(1) 的灵敏度计算公式中；

(6)令步骤(5)公式中的状态灵敏度

q_b、

和λ_b的同类项的系数为零，得到一系列伴随变量方程，以及消掉状态灵敏度后的灵敏度计算公式Ψ_b3；

(7)采用广义α算法求解步骤(6)的伴随变量方程，得到伴随变量μ、ν、σ⁰、ρ⁰、ξ⁰、σ^f、ρ^f、ξ^f、π和β；

(8)将步骤(6)的灵敏度计算公式中的关于设计变量的求导项，分别进行总体级别的局部有限差分，得到灵敏度计算公式Ψ_b4；

(9)将步骤(7)得到的伴随变量，代入步骤(8)的灵敏度计算公式中，即可解得柔性多体系统的灵敏度。

仿真实例：利用本发明方法，针对柔性绳索系统算例，展开数值仿真。

图3为柔性绳索系统动力学模型。采用绝对节点坐标方法建模，系统的状态变量为单元节点的位置坐标和斜率坐标。绳索系统的原长L为1.2m，横截面积A为0.0018m²，杨氏模量 E为7MPa,截面惯性矩I为1.215m⁴，密度ρ为5540kg/m³。在分析过程中，将绳索系统等分为10个单元，单元的长度l为0.12m。假设柔性绳索系统的初始状态为水平放置，其在A处被固定，只受重力作用，则系统在1s时间范围内的仿真结果如图4所示。图中曲线分别表示柔性绳索系统在不同时刻下的当前构型。取设计变量b＝[L,A,ρ,E]^T，分别表示为柔性绳索的长度、横截面积、密度和弹性模量。当系统的起止时刻t⁰＝0,t^f＝1时，给定目标函数

其中，x为绳索末端节点在水平方向的位置坐标，L为绳索原长。接下来，分别针对直接微分法和伴随变量法，采用本发明方法，进行柔性绳索系统的灵敏度计算。

表1当δ＝1E-05时，不同灵敏度计算方法的灵敏度计算结果和计算时间

表1为半解析方法的扰动值δ取1E-05时，针对直接微分法(Direct)和伴随变量法(Adjoint)，分别采用解析法(AM)、传统半解析方法(TSAM)及本发明提出的改进半解析方法(PSAM)，求解柔性绳索系统的目标函数对设计变量的灵敏度计算结果和计算时间。从中可以得到：分别针对直接微分法和伴随变量法进行计算，三种灵敏度分析方法对设计变量的灵敏度计算结果基本一致，证明了本发明的柔性多体系统动力学半解析灵敏度分析方法的正确性；考虑计算时间，半解析灵敏度分析方法的计算时间比解析方法长，但相比于传统半解析灵敏度分析方法，本发明提出的改进半解析灵敏度分析方法计算用时较少，并且其计算用时与解析法相差很小。

图5为取不同的扰动值δ＝1E-01,1E-02,…,1E-12时，分别采用解析方法(AM)、本发明提出的改进半解析灵敏度分析方法(PSAM)和有限差分法(FDM)计算柔性绳索系统的目标函数对4个设计变量的灵敏度计算结果。从中可以得到：有限差分法的计算结果受扰动量大小影响明显，就绳索长度而言，有限差分法的灵敏度计算结果与解析解相差很大，就绳索的横截面积、密度和弹性模量而言，有限差分法的灵敏度计算结果与解析解相差不大，但只能在较小的扰动范围内保持数值稳定性；而本发明提出的改进半解析灵敏度分析方法的计算结果受扰动量大小影响较小，其计算结果与解析解基本一致，并且可以在较大的扰动范围内保持数值稳定性。

图6为半解析方法的扰动值δ取1E-05时，分别采用解析方法(AM)、本发明提出的改进半解析灵敏度分析方法(PSAM)和有限差分法(FDM)计算柔性绳索系统在1s时间范围内，目标函数对4个设计变量的灵敏度计算结果随时间变化曲线。从中可以得到：对于4个设计变量，有限差分法的计算结果始终与解析解相差很大；而本发明提出的改进半解析灵敏度分析方法的计算结果始终与解析解保持一致，证明了本发明的柔性多体系统动力学半解析灵敏度分析方法具有较高的计算精度。

图7为半解析方法的扰动值δ取1E-05时，分别采用解析方法(AM)和本发明提出的改进半解析灵敏度分析方法(PSAM)计算柔性绳索系统在1s时间范围内，绳索末端点x方向的位置坐标对4个设计变量的状态灵敏度计算结果随时间变化曲线。从中可以得到：对于4个设计变量，本发明提出的改进半解析灵敏度分析方法的计算结果始终与解析解保持一致，也证明了本发明的柔性多体系统动力学半解析灵敏度分析方法具有较高的计算精度。

表2不同规模柔性绳索系统下，不同灵敏度计算方法的计算时间(s)

表2为半解析方法的扰动值δ取1E-05时，分别针对直接微分法(Direct)和伴随变量法 (Adjoint)，采用本发明提出的改进半解析灵敏度分析方法(PSAM)和有限差分法(FDM) 计算不同规模柔性绳索系统在1s时间范围内的灵敏度的计算用时。其中，最左边一列为柔性绳索系统划分单元的数目。从中可以得到：分别针对直接微分法和伴随变量法，本发明提出的改进半解析灵敏度分析方法计算用时较少，其中，伴随变量法用时最少；而有限差分法的灵敏度计算时间相对较长，是其它两种方法的1-6倍，证明了本发明的柔性多体系统动力学半解析灵敏度分析方法具有较高的计算精度。

Claims (7)

1.一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，首先，基于绝对节点坐标方法，建立柔性多体系统的质量矩阵、刚度矩阵和广义力列阵；其次，建立柔性多体系统的动力学方程和优化目标函数；再次，基于直接微分法或伴随变量法，建立柔性多体系统动力学的半解析灵敏度计算公式；最后，求解柔性多体系统动力学微分代数方程，获得灵敏度计算结果；包括下述步骤：

一、基于绝对节点坐标方法，建立柔性多体系统单元的质量矩阵、刚度矩阵和广义力列阵；

绝对节点坐标方法在全局坐标系下描述柔性多体系统的运动形态，将每个节点的广义坐标描述为位置矢量坐标和该点的斜率矢量坐标；基于三维建模的一维二节点梁单元，假设柔性多体系统节点的个数为n，每个节点有6个坐标，分别为3个位置矢量坐标和3个斜率矢量坐标，每个单元包含2个节点；

所述梁单元坐标表示为e＝[e_u e_v]^T，其中，e_p(p＝u,v)表示节点p的坐标：

式中，r_pX,r_pY,r_pZ(p＝u,v)分别表示节点u,v的位置矢量在全局坐标系中的分量；

分别表示节点u,v的斜率坐标；

因此，梁单元上任意点的位置表示为：

式中，

表示定义在总体坐标系上的单元形函数矩阵，写为：

S＝[S₁I,S₂I,S₃I,S₄I] (3)

其中，

为单位矩阵；形函数S₁＝1-3ξ²+2ξ³，S₂＝L(ξ-2ξ²+ξ³)，S₃＝3ξ²-2ξ³，S₄＝L(ξ³-ξ²)，其中，变量ξ＝x/L,x∈[0,L]；L为单元的初始长度，x为单元的局部坐标；

(1)柔性多体系统单元的质量矩阵

采用绝对节点坐标方法推导的柔性多体系统的质量矩阵是常数矩阵；根据单元的动能计算公式，能够得到单元的质量矩阵

式中，ρ为单元的密度；A为单元的横截面积；

(2)柔性多体系统单元的广义力列阵

柔性多体系统的广义力列阵包括广义外力列阵和广义弹性力列阵；其中，广义外力为重力，广义弹性力包括纵向拉伸变形和横向弯曲变形；

假设系统只受到重力作用，通过求解重力对系统单元产生的虚功，得到单元的广义重力

式中，g为重力加速度；

假设梁单元为各向同性，单元的应变能包括与纵向拉伸变形相关的应变能U_l和与横向弯曲变形相关的应变能U_t；通过将单元的应变能对单元坐标求导，得到单元的广义弹性力

式中，E为单元的弹性模量；ε为单元的纵向拉伸应变；I为单元的截面惯性矩；下标x表示单元形函数对物质坐标x的导数；双下标xx表示单元形函数对物质坐标x的二阶导数；

(3)柔性多体系统单元的刚度矩阵

单元的刚度矩阵

表示为：

二、建立柔性多体系统的动力学方程和优化目标函数；

柔性多体系统的动力学方程表示为：

式中，

为系统的广义质量矩阵；Φ＝[Φ₁,Φ₂,…,Φ_m]^T为系统的位置约束列阵；

为位置约束的Jacobian矩阵；

为动力学方程中的Lagrange乘子矢量；q、

和

分别为系统的广义坐标、广义速度和广义加速度；

为系统的广义外力列阵；

为系统的广义弹性力列阵，其依赖于系统的设计变量；

在多体系统动力学优化问题中，目标函数表示为如下的积分形式：

式中，前两部分G⁰和G^f分别与系统的初态和终态有关，第三部分H为被积分项，其与系统的中间过程有关；其中，上标0和f分别表示相关参量的初始和终止时刻值；q＝[q₁,q₂,…,q_6n]^T为状态变量，描述系统的动态响应；b＝[b₁,b₂,…,b_c]^T为设计变量，是优化设计过程中计算选择并需要最终确定的待求参数，其下标c表示设计变量的数目；t⁰和t^f分别表示初始时刻和终止时刻，也可以表示某些状态变量或其速度达到某一定值的特定时刻，其由下式确定：

式中，Ω⁰和Ω^f分别表示初始和终止时刻条件；

系统的初始状态依赖于设计变量，其应满足以下相容附加条件：

式中，

和

分别表示初始状态位置和速度的相容附加条件，其需要使得

和

满秩；

三、基于直接微分法或伴随变量法，建立柔性多体系统动力学的半解析灵敏度计算公式；

(a)基于直接微分法建立柔性多体系统动力学的半解析灵敏度计算公式

(1)建立柔性多体系统的半解析灵敏度计算公式

改进半解析灵敏度分析方法，是将灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分；其灵敏度计算公式表示为：

其中，

式中，下标表示对相应变量的导数；其中，q_b、

和λ_b为状态灵敏度，分别表示为状态变量相关量q、

qⁱ和λ_b(i＝0,f)对设计变量的导数；

和

分别为系统起始和终止时刻对设计变量的导数，由系统的起始和终止时刻条件确定；

和

为总体级别的局部有限差分项；

(2)建立柔性多体系统的状态灵敏度计算公式

采用直接微分法求解柔性多体系统的灵敏度，需要求解系统的状态灵敏度

q_b、

和λ_b；为保证计算结果的准确性，采用约束违约自动稳定方法求解；

改进半解析灵敏度分析方法，是将状态灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分；其状态灵敏度计算公式表示为：

其中，

式中，

和

为总体级别的局部有限差分项；为了与位置约束Φ进行区分，将

定义为Baumgarte约束；η₁和η₂为自定义参数；

传统半解析灵敏度分析方法，是将状态灵敏度计算解析公式中的对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

其中，Δb表示对设计变量的扰动；M_e为单元级别的质量矩阵；Φ_qe为单元级别的约束Jacobian矩阵；Q_se为单元级别的广义外力列阵；Q_ke为单元级别的广义弹性力列阵；

为单元级别的Baumgarte约束；

和λ_e分别为单元加速度矢量和Lagrange乘子矢量；

(3)建立柔性多体系统的状态灵敏度初值计算公式

求解上述状态灵敏度计算公式，必须要先给定状态灵敏度的初值

和

改进半解析灵敏度分析方法，是将状态灵敏度的初值计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分；其状态灵敏度的初值计算公式表示为：

式中，

为系统的速度级约束；

和

为总体级别的局部有限差分项；

传统半解析灵敏度分析方法，是将状态灵敏度的初值计算解析公式中的对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

其中，

和

分别为系统初始状态单元级别的位置约束和速度级约束；

和

为系统单元级别的初始相容附加条件；

(b)基于伴随变量法建立柔性多体系统动力学的半解析灵敏度计算公式

伴随变量法是通过引进一系列伴随变量

和

将其分别转置左乘柔性多体系统的动力学方程、约束方程和起始条件对设计变量的导数方程；然后，将所有左乘伴随变量的方程分别累减到两次分部积分后的柔性多体系统的灵敏度计算公式中；最后，消掉式中的状态灵敏度

q_b、

和λ_b，得到一系列伴随变量方程和新的灵敏度计算公式；该方法得到的灵敏度计算公式包含伴随变量，不包含状态灵敏度；

(1)求解伴随变量

推导得到的伴随变量方程分别为：

通过求解上述一系列伴随变量方程，得到相应的伴随变量；

(2)建立柔性多体系统的半解析灵敏度计算公式

改进半解析灵敏度分析方法，是将灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分；其灵敏度计算公式表示为：

式中，

和

为总体级别的局部有限差分项；

传统半解析灵敏度分析方法，是将灵敏度计算解析公式中的对设计变量的求导项分别进行单元级别的局部有限差分，

和

分别表示为

和

其中，Φ_e为单元级别的位置约束；M_e、Φ_qe、Q_se、Q_ke、

和λ_e(i＝0,f)都是基于单元级别；μ、ν、π、β、σⁱ和ρⁱ(i＝0,f)是伴随变量；

四、求解柔性多体系统动力学微分代数方程组，获得灵敏度计算结果；

针对具有大变形、大位移的柔性多体系统动力学问题，采用隐式方法进行求解；

采用广义α算法进行求解，令Q＝Q_s-Q_k，则离散DAEs的时间区间后，得到如下非线性代数方程组：

式中，h表示积分时间步长；上标“(d)”表示第d个迭代步；ζ和γ为算法参数；α为引入的算法辅助基矢量，其满足：

其中，a_m和a_f为算法参数；

为保证计算结果的精度和数值稳定性，算法参数a_m、a_f、ζ和γ分别取为：

式中，ρ_r∈[0,1]是谱半径，其取值大小决定算法能量耗散的频率范围；令辅助基矢量和Lagrange乘子矢量组成的列向量为未知变量，采用Newton–Raphson法迭代求解上述非线性代数方程组；随时间步进，得到仿真时间内的柔性多体系统的灵敏度计算结果。

2.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤一中，单元的纵向拉伸应变ε的具体表达式为：

其中，l为单元的初始纵向长度；l_s为变形后单元的纵向长度，由对微元弧长积分得到：

3.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤三(a)中，采用改进半解析灵敏度分析方法，柔性多体系统的半解析灵敏度计算公式的局部有限差分项

和

分别表示为：

4.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤三(a)中，采用改进半解析灵敏度分析方法，柔性多体系统的状态灵敏度计算公式的局部有限差分项

和

分别表示为：

其中，M、Φ_q、Q_s、Q_k、

和λ都是基于总体级别。

5.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤三(a)中，采用改进半解析灵敏度分析方法，柔性多体系统的状态灵敏度初值计算公式的局部有限差分项

和

分别表示为：

其中，Φ⁰、

和

都是基于总体级别。

6.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤三(b)中，需要将柔性多体系统的灵敏度计算公式分别关于

和

依次进行分部积分：

7.根据权利要求1所述的一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，其特征在于，在步骤三(b)中，采用改进半解析灵敏度分析方法，柔性多体系统的半解析灵敏度计算公式的局部有限差分项

和

分别表示为：

其中，M、Φ_q、Q_s、Q_k、Φ、Φⁱ、

和λ(i＝0,f)都是基于总体级别；μ、ν、π、β、σⁱ、ρⁱ和ξⁱ(i＝0,f)是伴随变量。

Images