CN111159636B - 一种柔性多体***动力学半解析灵敏度分析方法 - Google Patents
一种柔性多体***动力学半解析灵敏度分析方法 Download PDFInfo
- Publication number
- CN111159636B CN111159636B CN201911224272.9A CN201911224272A CN111159636B CN 111159636 B CN111159636 B CN 111159636B CN 201911224272 A CN201911224272 A CN 201911224272A CN 111159636 B CN111159636 B CN 111159636B
- Authority
- CN
- China
- Prior art keywords
- body system
- flexible multi
- sensitivity
- analytic
- semi
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 174
- 238000010206 sensitivity analysis Methods 0.000 title claims abstract description 72
- 230000035945 sensitivity Effects 0.000 claims abstract description 141
- 238000004364 calculation method Methods 0.000 claims abstract description 117
- 239000011159 matrix material Substances 0.000 claims abstract description 44
- 238000005457 optimization Methods 0.000 claims abstract description 15
- 238000013461 design Methods 0.000 claims description 86
- 239000013598 vector Substances 0.000 claims description 27
- 238000004422 calculation algorithm Methods 0.000 claims description 16
- 230000005484 gravity Effects 0.000 claims description 13
- 230000010354 integration Effects 0.000 claims description 11
- 238000009795 derivation Methods 0.000 claims description 10
- 230000001133 acceleration Effects 0.000 claims description 8
- 238000004088 simulation Methods 0.000 claims description 6
- 230000008569 process Effects 0.000 claims description 5
- 238000005452 bending Methods 0.000 claims description 4
- 230000004044 response Effects 0.000 claims description 3
- 238000001228 spectrum Methods 0.000 claims description 3
- 239000000126 substance Substances 0.000 claims description 3
- 238000002940 Newton-Raphson method Methods 0.000 claims description 2
- 230000009471 action Effects 0.000 claims description 2
- 238000012938 design process Methods 0.000 claims description 2
- 238000006073 displacement reaction Methods 0.000 claims description 2
- 239000000463 material Substances 0.000 claims description 2
- 230000021715 photosynthesis, light harvesting Effects 0.000 claims description 2
- 230000000087 stabilizing effect Effects 0.000 claims 1
- 230000009466 transformation Effects 0.000 claims 1
- 238000004458 analytical method Methods 0.000 description 29
- 230000004069 differentiation Effects 0.000 description 6
- 238000011160 research Methods 0.000 description 3
- 238000003491 array Methods 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 230000006641 stabilisation Effects 0.000 description 2
- 238000011105 stabilization Methods 0.000 description 2
- IOYNQIMAUDJVEI-BMVIKAAMSA-N Tepraloxydim Chemical group C1C(=O)C(C(=N/OC\C=C\Cl)/CC)=C(O)CC1C1CCOCC1 IOYNQIMAUDJVEI-BMVIKAAMSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000008033 biological extinction Effects 0.000 description 1
- 239000002131 composite material Substances 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000013433 optimization analysis Methods 0.000 description 1
- 230000008092 positive effect Effects 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000001052 transient effect Effects 0.000 description 1
- 238000004800 variational method Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Complex Calculations (AREA)
Abstract
本发明涉及动力学***优化技术领域,提出一种基于绝对节点坐标描述的柔性多体***动力学半解析灵敏度分析方法。首先,基于绝对节点坐标方法,建立柔性多体***的质量矩阵、刚度矩阵和广义力列阵;其次,建立柔性多体***的动力学方程和优化目标函数;再次,基于直接微分法或伴随变量法,建立柔性多体***动力学的半解析灵敏度计算公式;最后,求解柔性多体***动力学微分代数方程,获得灵敏度计算结果。本发明根据绝对节点坐标方法和多体***动力学理论,建立柔性多体***的半解析灵敏度计算公式,以解决柔性多体***动力学的灵敏度分析问题,目的在于提供一套柔性多体***灵敏度分析的新策略,为大规模复杂的柔性多体***的灵敏度计算提供便利。
Description
技术领域
本发明属于多体***动力学优化技术领域,涉及一种柔性多体***动力学的灵敏度分析方法,尤其涉及一种基于绝对节点坐标描述的柔性多体***动力学半解析灵敏度分析方法。
背景技术
近年来,多体***动力学及其优化分析在航空航天、机械、汽车等领域发挥着越来越重要的作用。如果使用基于梯度的优化算法进行多体***动力学的优化设计,则需要对多体***进行灵敏度分析。灵敏度分析不仅可以用于确定部分优化算法的迭代方向,也可以表征设计变量对目标函数的影响程度,从而可以通过降低设计变量的个数来提高优化效率。因此,多体***动力学的灵敏度分析已成为多体***优化设计的核心问题。
目前,对多体***灵敏度分析的研究大多是针对刚性多体***,由于柔性多体***的高维和强非线性,使其求解相对困难,因此,对柔性多体***灵敏度分析的研究相对较少。柔性多体***的建模方法主要有增量有限元法、浮动坐标方法、大转动矢量方法和绝对节点坐标方法。相比于其它三种方法,绝对节点坐标方法有效地避开了柔性体有限转动的参数化问题,在全局坐标系下采用插值函数来描述大变形和大范围转动,能精确地描述柔性多体***高度的几何非线性。
相比于传统静态优化设计,多体***动力学动态优化设计的目标函数和约束方程包含状态变量,因此,对其进行优化设计的灵敏度包括状态灵敏度(状态变量对设计变量的灵敏度)和设计灵敏度(目标函数对设计变量的灵敏度)两部分。并且,优化设计还受多体***状态方程的约束,其状态方程往往表现为一组高维和强非线性的微分代数方程组,给灵敏度的求解带来了困难。目前,对多体***动力学灵敏度分析的主要方法有:有限差分法、直接微分法和伴随变量法。
有限差分法是一种近似计算方法。它只需要对设计变量做扰动,然后采用差商的方式计算目标函数对设计变量的灵敏度。但其计算的工作量会随设计变量的数量成比例增加,计算效率和精度相对较低。1991年,Greene等采用有限差分法计算了线性、结构和瞬态响应问题的灵敏度,发现这种方法在实际计算时,对计算程序时间和精度的要求很高。直接微分法是通过直接将多体***的动力学方程和约束方程对设计变量求导,计算出状态灵敏度,从而求解多体***的灵敏度。该方法由Krishnaswami和Bhatti在1984年首先提出,随后,Chang等人研究了采用直接微分法计算约束动力***设计灵敏度系数矩阵的一般方法,并通过两个算例说明了该方法的有效性;Dias等人采用直接微分法建立了刚-柔多体***灵敏度分析的灵敏度方程,并通过两个刚-柔多体***的例子,将计算结果与有限差分法进行了比较;Neto等人将直接微分法应用于具有复合材料的柔性多体***的灵敏度求解。伴随变量法是通过引入一系列伴随变量,消掉方程中包含状态灵敏度的相关项,从而求解得到的一系列伴随变量方程来计算灵敏度。该方法于1978年由Haug和Arora首次提出,因其具有计算速度较快的特点,在近几年应用较为广泛。Liu采用伴随变量法推导出了约束柔性多体***的一阶和二阶灵敏度分析方程;Zhang和Chen采用伴随变量法,将以大变形问题梁板建模为重点的绝对节点坐标公式扩展到柔性多体***的设计灵敏度分析中;Alexander等人采用浮动坐标方法进行柔体***建模,将伴随变量法应用于具有运动回路的柔性多体***。
然而,上述提到的直接微分法和伴随变量法是采用解析方法来计算设计变量的灵敏度,当多体动力学***规模很大、结构较复杂时,采用解析方法进行灵敏度计算所需推导的动力学公式将会十分繁杂,甚至可能会出现无法解析求得多体***对某些设计变量的灵敏度的情况。采用半解析灵敏度分析方法,可以为大规模复杂的柔性多体***的灵敏度计算提供便利。
发明内容
本发明提出一种基于绝对节点坐标描述的柔性多体***动力学半解析灵敏度分析方法。该方法基于绝对节点坐标方法,建立柔性多体***的动力学模型,以解决柔性多体***的灵敏度分析问题,目的在于提供一种柔性多体***灵敏度分析的新策略,为大规模复杂的柔性多体***的灵敏度计算提供便利。
为了达到上述目的,本发明采用的技术方案为:
一种柔性多体***动力学半解析灵敏度分析方法,步骤如下:
一、基于绝对节点坐标方法,建立柔性多体***单元的质量矩阵、刚度矩阵和广义力列阵
绝对节点坐标方法在全局坐标系下描述柔性多体***的运动形态,将每个节点的广义坐标描述为位置矢量坐标和该点的斜率矢量坐标。本发明基于三维建模的一维二节点梁单元:如图2所示,每个节点有6个坐标,分别为3个位置矢量坐标和3个斜率矢量坐标,每个单元包含2个节点。假设柔性多体***节点的个数为n。
所述梁单元坐标表示为e=[eu ev]T,其中,ep(p=u,v)表示节点p的坐标:
因此,梁单元上任意点的位置可以表示为:
S=[S1I,S2I,S3I,S4I] (3)
其中,为单位矩阵;形函数S1=1-3ξ2+2ξ3,S2=L(ξ-2ξ2+ξ3),S3=3ξ2-2ξ3,S4=L(ξ3-ξ2),其中,变量ξ=x/L,x∈[0,L]。L为单元的初始长度,x为单元的局部坐标。
1、建立柔性多体***单元的质量矩阵
式中,ρ为单元的密度;A为单元的横截面积;V为单元的体积。
2、建立柔性多体***单元的广义力列阵
柔性多体***的广义力列阵包括广义外力列阵和广义弹性力列阵。其中,广义外力列阵包括重力和接触力等,广义弹性力列阵包括纵向拉伸变形和横向弯曲变形。
假设***只受到重力作用,通过虚功原理,可以得到重力虚功:
其中,g为重力加速度。
假设梁单元为各向同性,单元的总应变能U包括与纵向拉伸变形相关的应变能Ul和与横向弯曲变形相关的应变能Ut:
式中,ε为单元的纵向拉伸应变;κ为单元在当前构型下的曲率;E为单元的弹性模量;I为单元的截面惯性矩。
式中,下标x表示单元形函数对物质坐标x的导数;双下标xx表示单元形函数对物质坐标x的二阶导数。单元的纵向拉伸应变ε的具体表达式为:
其中,l为单元的初始纵向长度;ls为变形后单元的纵向长度,可以通过对微元弧长积分得到:
3、建立柔性多体***单元的刚度矩阵
二、建立柔性多体***的动力学方程和优化目标函数
本发明研究的问题是基于微分代数方程形式的柔性多体***动力学数学模型,***的动力学方程表示为:
式中,为***的广义质量矩阵;Φ=[Φ1,Φ2,…,Φm]T为***的位置约束列阵;为位置约束的Jacobian矩阵;为动力学方程中的Lagrange乘子矢量;q、和分别为***的广义坐标、广义速度和广义加速度;为***的广义外力列阵;为***的广义弹性力列阵,其依赖于***的设计变量。
目标函数是评价优化设计方案的标准,在多体***动力学优化问题中,其一般表示为如下的积分形式:
式中,前两部分G0和Gf分别与***的初态和终态有关,第三部分H为被积分项,其与***的中间过程有关。其中,上标0和f分别表示相关参量的初始和终止时刻值;q=[q1,q2,…,q6n]T为状态变量,描述***的动态响应;b=[b1,b2,…,bc]T为设计变量,是优化设计过程中计算选择并需要最终确定的待求参数,其下标c表示设计变量的数目;t0和tf分别表示初始时刻和终止时刻,也可以表示某些状态变量或其速度达到某一定值的特定时刻,一般地,其可以由下式确定:
式中,Ω0和Ωf分别表示初始和终止时刻条件。
***的初始状态依赖于设计变量,其应满足以下相容附加条件:
三、基于直接微分法或伴随变量法,建立柔性多体***动力学的半解析灵敏度计算公式
本发明提出的柔性多体***动力学的半解析灵敏度分析方法,相较于解析方法,其不需要解析地推导优化设计目标函数和柔性多体***动力学方程对设计变量的导数,而是将解析求导公式中对设计变量求导的部分用局部有限差分代替。传统半解析灵敏度分析方法是基于单元级别进行灵敏度计算,而改进半解析灵敏度分析方法是基于总体级别,其无需提取扰动前后柔性多体***单元的相关矩阵信息,使程序实现过程更为简便,具有较高的计算效率。
(a)基于直接微分法建立柔性多体***动力学的半解析灵敏度计算公式
1、建立柔性多体***的半解析灵敏度计算公式
灵敏度是目标函数对设计变量的偏导数,一般形式的柔性多体***灵敏度计算的解析公式可以表示为:
式中,下标表示对相应变量的导数。其中,qb、和λb为状态灵敏度,分别表示为状态变量相关量q、qi和λb(i=0,f)对设计变量的导数;和分别为***起始和终止时刻对设计变量的导数,由***的起始和终止时刻条件确定,将约束条件式(14)对设计变量求导,可得到:
本发明采用改进半解析灵敏度分析方法,将上述解析方法中对设计变量的求导项分别进行总体级别的局部有限差分,可得到改进半解析灵敏度分析方法的灵敏度计算公式:
式中,Δb表示对设计变量的扰动。
2、建立柔性多体***的状态灵敏度计算公式
采用直接微分法求解柔性多体***的灵敏度,需要求解***的状态灵敏度。将柔性多体***的动力学方程对设计变量求导,可得:
为保证计算结果的准确性,采用约束违约自动稳定方法求解,此时,***的约束方程同时包含位置约束、速度级约束及加速度级约束:
将式(25)对设计变量求导,并与式(24)中的第一个式子联立新的矩阵微分代数方程组:
其中
本发明采用改进半解析灵敏度分析方法,将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分,可得到改进半解析灵敏度分析方法的状态灵敏度计算公式:
式中,Me为单元级别的质量矩阵;Φqe为单元级别的位置约束Jacobian矩阵;Qse为单元级别的广义外力列阵;Qke为单元级别的广义弹性力列阵;为单元级别的Baumgarte约束;和λe分别为单元加速度矢量和Lagrange乘子矢量。
3、建立柔性多体***的状态灵敏度初值计算公式
本发明采用改进半解析灵敏度分析方法,将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分,可得到改进半解析灵敏度分析方法的状态灵敏度初值计算公式:
(b)基于伴随变量法建立柔性多体***动力学的半解析灵敏度计算公式
1、引进伴随变量
2、建立伴随变量方程
通过求解上述一系列伴随变量方程,即可求得相应的伴随变量。
3、建立柔性多体***的半解析灵敏度计算公式
采用解析方法求解上述微分方程,可解析地求得柔性多体***的灵敏度。
本发明采用改进半解析灵敏度分析方法,将上述解析公式中对设计变量的求导项分别进行总体级别的局部有限差分,可得到改进半解析灵敏度分析方法的灵敏度计算公式:
四、求解柔性多体***动力学微分代数方程组,获得灵敏度计算结果
基于绝对节点坐标描述的含约束柔性多体***的灵敏度计算公式,是指标3的微分代数方程组(DAEs),其计算方法分为显示方法和隐式方法。显示方法不需要平衡迭代,计算速度较快,但该方法是条件稳定的,其计算结果受积分步长影响较大,收敛性差;隐式方法需要在积分步内迭代求解,计算速度较慢,但该方法是无条件稳定的,可以选取较大的积分步长进行计算。针对本发明的具有大变形、大位移的柔性多体***动力学问题,通常采用隐式方法进行求解。
本发明采用广义α算法进行求解。考虑公式的简洁性,令Q=Qs-Qk,则离散DAEs的时间区间后,可得到如下非线性代数方程组:
式中,h表示积分时间步长;上标“(d)”表示第d个迭代步;ζ和γ为算法参数;α为引入的算法辅助基矢量,其满足:
其中,am和af为算法参数。
为了保证计算结果的精度和数值稳定性,算法参数am、af、ζ和γ分别取为:
式中,ρ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,本发明的具体实施方式如下:
实施方式基于描述柔性多体***的绝对节点坐标建模方法和柔性多体***动力学理论,可用于柔性多体***动力学灵敏度问题的分析,具体按以下步骤实现:
一、基于绝对节点坐标方法,推导柔性多体***单元的质量矩阵Me、刚度矩阵Ke及广义力列阵Qge和Qke。
二、基于柔性多体***动力学理论,在步骤1的基础上,建立柔性多体***的动力学方程。具体步骤为:
(1)给出该柔性多体***的尺寸参数、划分单元的数目和单元材料信息;
(2)给出该柔性多体***的初始状态信息,包括初始位置、初始速度等;
(3)给出该柔性多体受到的重力等外力信息;
(4)建立该柔性多体***总体的质量矩阵、刚度矩阵和广义力列阵;
(5)给出该柔性多体***的约束方程;
(6)在步骤(1)-(5)的基础上,便可建立柔性多体***的动力学方程,即微分代数方程组(DAEs)。
三、给出柔性多体***优化的设计变量,以及优化的初态和终态约束条件,建立***的优化目标函数。
四、基于直接微分法或伴随变量法,采用改进半解析方法,求解柔性多体***的灵敏度。
(a)基于直接微分法求解柔性多体***的灵敏度。具体步骤为:
(1)将***的目标函数对设计变量求导,并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分,得到灵敏度计算公式Ψb1;
(4)将***的位置约束和速度级约束对设计变量求导,并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分;
(6)采用约束违约自动稳定方法,重新建立***的约束方程;
(7)将***的动力学方程对设计变量求导,并将公式中关于设计变量的求导项分别进行总体级别的局部有限差分;
(9)将步骤(8)得到的***的状态灵敏度,代入步骤(1)中的灵敏度计算公式,即可解得柔性多体***的灵敏度。
(b)基于伴随变量法求解柔性多体***的灵敏度。具体步骤为:
(2)引入伴随变量μ和ν,将它们分别转置左乘***的动力学方程,并在t0到tf上积分;
(4)引进伴随变量σi、ρi、ξi、π和β(i=0,f)将它们分别转置左乘***的位置约束、速度级约束和起始条件对设计变量的导数方程;
(5)将步骤(2)得到的第二个方程、步骤(3)和步骤(4)得到的方程,分别累减到步骤(1)的灵敏度计算公式中;
(7)采用广义α算法求解步骤(6)的伴随变量方程,得到伴随变量μ、ν、σ0、ρ0、ξ0、σf、ρf、ξf、π和β;
(8)将步骤(6)的灵敏度计算公式中的关于设计变量的求导项,分别进行总体级别的局部有限差分,得到灵敏度计算公式Ψb4;
(9)将步骤(7)得到的伴随变量,代入步骤(8)的灵敏度计算公式中,即可解得柔性多体***的灵敏度。
仿真实例:利用本发明方法,针对柔性绳索***算例,展开数值仿真。
图3为柔性绳索***动力学模型。采用绝对节点坐标方法建模,***的状态变量为单元节点的位置坐标和斜率坐标。绳索***的原长L为1.2m,横截面积A为0.0018m2,杨氏模量E为7MPa,截面惯性矩I为1.215m4,密度ρ为5540kg/m3。在分析过程中,将绳索***等分为10个单元,单元的长度l为0.12m。假设柔性绳索***的初始状态为水平放置,其在A处被固定,只受重力作用,则***在1s时间范围内的仿真结果如图4所示。图中曲线分别表示柔性绳索***在不同时刻下的当前构型。取设计变量b=[L,A,ρ,E]T,分别表示为柔性绳索的长度、横截面积、密度和弹性模量。当***的起止时刻t0=0,tf=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=[eu ev]T,其中,ep表示节点p的坐标,其中p=u,v:
因此,梁单元上任意点的位置表示为:
S=[S1I,S2I,S3I,S4I] (3)
其中,为单位矩阵;形函数S1=1-3ξ2+2ξ3,S2=L(ξ-2ξ2+ξ3),S3=3ξ2-2ξ3,S4=L(ξ3-ξ2),其中,变量ξ=x/L,x∈[0,L];L为单元的初始长度,x为单元的局部坐标;
(1)柔性多体***单元的质量矩阵
式中,ρ为单元的密度;A为单元的横截面积;
(2)柔性多体***单元的广义力列阵
柔性多体***的广义力列阵包括广义外力列阵和广义弹性力列阵;其中,广义外力为重力,广义弹性力包括纵向拉伸变形和横向弯曲变形;
式中,g为重力加速度;
式中,E为单元的弹性模量;ε为单元的纵向拉伸应变;I为单元的截面惯性矩;下标x表示单元形函数对物质坐标x的导数;双下标xx表示单元形函数对物质坐标x的二阶导数;
(3)柔性多体***单元的刚度矩阵
二、建立柔性多体***的动力学方程和优化目标函数;
柔性多体***的动力学方程表示为:
式中,为***的广义质量矩阵;Φ=[Φ1,Φ2,…,Φm]T为***的位置约束列阵;为位置约束的Jacobian矩阵;为动力学方程中的Lagrange乘子矢量;q、和分别为***的广义坐标、广义速度和广义加速度;为***的广义外力列阵;为***的广义弹性力列阵,其依赖于***的设计变量;
在多体***动力学优化问题中,目标函数表示为如下的积分形式:
式中,前两部分G0和Gf分别与***的初态和终态有关,第三部分H为被积分项,其与***的中间过程有关;其中,上标0和f分别表示相关参量的初始和终止时刻值;q=[q1,q2,…,q6n]T为状态变量,描述***的动态响应;b=[b1,b2,…,bc]T为设计变量,是优化设计过程中计算选择并需要最终确定的待求参数,其下标c表示设计变量的数目;t0和tf分别表示初始时刻和终止时刻,也可以表示某些状态变量或其速度达到某一定值的特定时刻,其由下式确定:
式中,Ω0和Ωf分别表示初始和终止时刻条件;
***的初始状态依赖于设计变量,其应满足以下相容附加条件:
三、基于直接微分法或伴随变量法,建立柔性多体***动力学的半解析灵敏度计算公式;
(a)基于直接微分法建立柔性多体***动力学的半解析灵敏度计算公式
(1)建立柔性多体***的半解析灵敏度计算公式
改进半解析灵敏度分析方法,是将灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分;其灵敏度计算公式表示为:
其中,
式(13)中,下标表示对相应变量的导数;其中,qb、和λb为状态灵敏度,分别表示为状态变量相关量q、qi和λb对设计变量的导数,其中i=0,f;和分别为***起始和终止时刻对设计变量的导数,由***的起始和终止时刻条件确定;和为总体级别的局部有限差分项;
(2)建立柔性多体***的状态灵敏度计算公式
改进半解析灵敏度分析方法,是将状态灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分;其状态灵敏度计算公式表示为:
其中,
其中,Δb表示对设计变量的扰动;Me为单元级别的质量矩阵;Φqe为单元级别的约束Jacobian矩阵;Qse为单元级别的广义外力列阵;Qke为单元级别的广义弹性力列阵;为单元级别的Baumgarte约束;和λe分别为单元加速度矢量和Lagrange乘子矢量;
(3)建立柔性多体***的状态灵敏度初值计算公式
改进半解析灵敏度分析方法,是将状态灵敏度的初值计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分;其状态灵敏度的初值计算公式表示为:
(b)基于伴随变量法建立柔性多体***动力学的半解析灵敏度计算公式
伴随变量法是通过引进一系列伴随变量 和将其分别转置左乘柔性多体***的动力学方程、约束方程和起始条件对设计变量的导数方程;然后,将所有左乘伴随变量的方程分别累减到两次分部积分后的柔性多体***的灵敏度计算公式中;最后,消掉式中的状态灵敏度qb、和λb,得到一系列伴随变量方程和新的灵敏度计算公式;该方法得到的灵敏度计算公式包含伴随变量,不包含状态灵敏度;
(1)求解伴随变量
推导得到的伴随变量方程分别为:
通过求解上述一系列伴随变量方程,得到相应的伴随变量;
(2)建立柔性多体***的半解析灵敏度计算公式
改进半解析灵敏度分析方法,是将灵敏度计算解析公式中的对设计变量的求导项分别进行总体级别的局部有限差分;其灵敏度计算公式表示为:
四、求解柔性多体***动力学微分代数方程组,获得灵敏度计算结果;
针对具有大变形、大位移的柔性多体***动力学问题,采用隐式方法进行求解;
采用广义α算法进行求解,令Q=Qs-Qk,则离散DAEs的时间区间后,得到如下非线性代数方程组:
式中,h表示积分时间步长;上标“(d)”表示第d个迭代步;ζ和γ为算法参数;α为引入的算法辅助基矢量,其满足:
其中,am和af为算法参数;
为保证计算结果的精度和数值稳定性,算法参数am、af、ζ和γ分别取为:
式中,ρr∈[0,1]是谱半径,其取值大小决定算法能量耗散的频率范围;令辅助基矢量和Lagrange乘子矢量组成的列向量为未知变量,采用Newton–Raphson法迭代求解上述非线性代数方程组;随时间步进,得到仿真时间内的柔性多体***的灵敏度计算结果。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911224272.9A CN111159636B (zh) | 2019-12-04 | 2019-12-04 | 一种柔性多体***动力学半解析灵敏度分析方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911224272.9A CN111159636B (zh) | 2019-12-04 | 2019-12-04 | 一种柔性多体***动力学半解析灵敏度分析方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111159636A CN111159636A (zh) | 2020-05-15 |
CN111159636B true CN111159636B (zh) | 2021-09-24 |
Family
ID=70556397
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911224272.9A Active CN111159636B (zh) | 2019-12-04 | 2019-12-04 | 一种柔性多体***动力学半解析灵敏度分析方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111159636B (zh) |
Families Citing this family (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN112084592A (zh) * | 2020-09-04 | 2020-12-15 | 上海交通大学 | 折叠式桁架动力学分析***、方法、装置和存储介质 |
CN112906212B (zh) * | 2021-02-05 | 2024-05-24 | 南京航空航天大学 | 一种基于绝对节点坐标法的裸电动力绳***建模方法 |
CN113076677B (zh) * | 2021-04-15 | 2023-03-17 | 朱礼云 | 基于五次埃尔米特形函数的柔性体结构高阶非线性有限元数值模拟方法 |
CN113297730B (zh) * | 2021-05-13 | 2023-04-11 | 广东工业大学 | 基于自适应模态的柔性多体***动态响应计算方法和*** |
CN113688479B (zh) * | 2021-08-27 | 2023-05-02 | 天津大学 | 电磁场有限元快速频率分析的电磁灵敏度分析方法 |
CN114357651A (zh) * | 2022-01-04 | 2022-04-15 | 北京理工大学 | 一种高效的不确定性柔性多体***动力学响应预测方法 |
Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101814103A (zh) * | 2010-04-01 | 2010-08-25 | 西北工业大学 | 基于超单元的多组件布局建模与结构优化设计方法 |
CN103034166A (zh) * | 2012-11-26 | 2013-04-10 | 北京工业大学 | 一种机床关键性几何误差源识别方法 |
CN104375460A (zh) * | 2014-11-17 | 2015-02-25 | 北京工业大学 | 一种数控机床加工精度可靠性敏感度分析方法 |
CN107122515A (zh) * | 2017-03-17 | 2017-09-01 | 北京航空航天大学 | 基于绝对节点坐标法的绳系运输***的动力学分析方法 |
CN107220421A (zh) * | 2017-05-18 | 2017-09-29 | 北京理工大学 | 一种空间复杂柔性结构多体***动力学建模与计算方法 |
CN107545126A (zh) * | 2017-09-28 | 2018-01-05 | 大连理工大学 | 一种基于多体***滑移绳索单元的聚合式张拉整体结构动力响应分析方法 |
CN107943748A (zh) * | 2017-11-22 | 2018-04-20 | 南京理工大学 | 基于Bathe积分策略的多体动力学方程求解方法 |
Family Cites Families (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US10685154B2 (en) * | 2016-12-29 | 2020-06-16 | Dassault Systemes Simulia Corp. | Analytical consistent sensitivities for nonlinear equilibriums, where the only source of nonlinearities is small sliding contact constraints |
-
2019
- 2019-12-04 CN CN201911224272.9A patent/CN111159636B/zh active Active
Patent Citations (7)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101814103A (zh) * | 2010-04-01 | 2010-08-25 | 西北工业大学 | 基于超单元的多组件布局建模与结构优化设计方法 |
CN103034166A (zh) * | 2012-11-26 | 2013-04-10 | 北京工业大学 | 一种机床关键性几何误差源识别方法 |
CN104375460A (zh) * | 2014-11-17 | 2015-02-25 | 北京工业大学 | 一种数控机床加工精度可靠性敏感度分析方法 |
CN107122515A (zh) * | 2017-03-17 | 2017-09-01 | 北京航空航天大学 | 基于绝对节点坐标法的绳系运输***的动力学分析方法 |
CN107220421A (zh) * | 2017-05-18 | 2017-09-29 | 北京理工大学 | 一种空间复杂柔性结构多体***动力学建模与计算方法 |
CN107545126A (zh) * | 2017-09-28 | 2018-01-05 | 大连理工大学 | 一种基于多体***滑移绳索单元的聚合式张拉整体结构动力响应分析方法 |
CN107943748A (zh) * | 2017-11-22 | 2018-04-20 | 南京理工大学 | 基于Bathe积分策略的多体动力学方程求解方法 |
Non-Patent Citations (3)
Title |
---|
A semi-analytical sensitivity analysis for multibody systems described using Level Sets;E.Tromme et al.;《International Conference on Engineering and Applied Sciences Optimization》;20140606;第1-27页 * |
基于绝对节点坐标法的柔性多体***灵敏度分析;王铁成 等;《振动与冲击》;20151231;第34卷(第24期);第89-92、103页 * |
结构优化半解析灵敏度分析的改进算法;陈飙松 等;《计算力学学报》;20181031;第35卷(第5期);第533-539页 * |
Also Published As
Publication number | Publication date |
---|---|
CN111159636A (zh) | 2020-05-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111159636B (zh) | 一种柔性多体***动力学半解析灵敏度分析方法 | |
CN113111430B (zh) | 基于非线性气动力降阶的弹性飞机飞行动力学建模方法 | |
CN109902404B (zh) | 不同阻尼形式的结构时程响应积分的统一递推计算方法 | |
CN107944137B (zh) | 高超声速飞行器弹道状态多场耦合的热气动弹性计算技术 | |
CN113191040A (zh) | 一种考虑结构稳定性的单材料结构拓扑优化方法和*** | |
Ruiz et al. | Eigenvector sensitivity when tracking modes with repeated eigenvalues | |
CN113204906B (zh) | 一种考虑结构稳定性的多相材料拓扑优化设计方法和*** | |
Albers et al. | Integrated structural and controller optimization in dynamic mechatronic systems | |
CN109325284B (zh) | 一种考虑不确定性的蜂窝结构舵面等效参数强稳健性辨识方法 | |
Lu et al. | A hybrid numerical method for vibration analysis of linear multibody systems with flexible components | |
CN110008635B (zh) | 利用Newmark精细积分结合法对弹塑性结构地震反应分析的方法 | |
CN104091003B (zh) | 一种基础运动时柔性壳结构大变形响应的有限元建模方法 | |
CN110795790B (zh) | 一种复杂建筑结构非线性动力时程分析方法 | |
CN112580241A (zh) | 一种基于结构降阶模型的非线性气动弹性动响应分析方法 | |
Hui et al. | A data-driven CUF-based beam model based on the tree-search algorithm | |
CN113505405B (zh) | 等效荷载获取方法、基于等效荷载的拓扑优化方法及*** | |
Awrejcewicz et al. | Nonlinear coupled problems in dynamics of shells | |
CN107562991B (zh) | 完全基于降阶模型的结构非线性屈曲平衡路径跟踪方法 | |
CN116468104A (zh) | 基于变分原理的面向神经算子训练与偏微分方程组求解的统一方法、介质及产品 | |
CN114818550B (zh) | 一种飞机振动试验中时变气动载荷地面等效模拟方法 | |
CN114211478B (zh) | 一种模块化机械臂协调操作最优控制方法及*** | |
Vidal et al. | Modeling of composite and sandwich beams with a generic cross-section using a variable separation method | |
CN112818583B (zh) | 等效静荷载获取方法、拓扑优化方法及*** | |
CN112966421B (zh) | 使用p型有限元法计算薄板结构屈曲载荷因子和相应屈曲形状的计算方法 | |
CN114692446A (zh) | 一种矢量转动球坐标参数化的非线性壳有限元建模方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |