本发明公开一种基于自适应模态的柔性多体系统动态响应计算方法，包括以下步骤：S1：将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；S2：设计柔性多体系统的运动方程的抽样矩阵；S3：求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；S4：利用贪婪高斯‑牛顿算法求解l1范数优化问题，得到柔性多体系统动态响应。本发明定义了可以用模态形状稀疏表示的弹性坐标，导出了具有少量运动方程约束的宽范数坐标下的l1范数优化问题，作为每个时间步长的动态响应分析。提出了GGN算法，有效地解决了这种优化问题，使模态坐标能够自适应地选择各种工作模式。
1.一种基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，包括以下步骤：S1：将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
S3：求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
S4：利用贪婪高斯‑牛顿算法求解l1范数优化问题，得到柔性多体系统动态响应；
S4.1：读取每个柔性体的有限元信息，计算每个柔性体的惯性形状积分模态，设置初始参数，所述惯性形状积分模态包括mff和Kff；
S4.2：求解 得到初始时刻的加速度 和拉格朗日乘子λ0，式中，M表示柔性多体系统的质量矩阵，G表示柔性多体系统的约束雅可比矩阵，表示柔性多体系统的加速度，λ表示柔性多体系统的拉格朗日乘子， 表示作用于柔性多体系统的广义外力， 表示与加速度有关的向量；
S4.3：对每个柔性体进行测量，得到其待定运动方程Φ(c)；
S4.4：利用贪婪高斯‑牛顿算法求解 约束，得到pt, 分别表示在t时刻的坐标、速度、加速度和拉格朗日乘子；
S4.5：判断t是否大于Time，若是，则结束计算，若不是，令t＝t+h，返回步骤S4.3。
2.根据权利要求1所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S1中将柔性多体系统的模态坐标定义为正交振型上的稀疏系数，具体为：设模态坐标是稀疏的，柔性多体系统中有N个柔性体，则广义坐标q表示如下：iT
式中，i表示第i个柔性体， 为柔性体i的弹性坐标的向量，R 为柔性iT i
体i相对于惯性系的位置向量，θ 为柔性体i相对于惯性系的角位移，I是单位矩阵，B是由柔性体i的全模态形状组成， 其中 表示柔性多体系统的弹性坐标， 表示柔性多体系统的模态坐标；
3.根据权利要求2所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S2中柔性多体系统的运动方程的抽样矩阵，具体为：式中， 表示相对于柔性体的位置和旋转情况的质量矩阵， 是耦合刚性运动和变形的质量矩阵， 表示单个柔性体与弹性坐标有关的质量矩阵， 表示与单个柔性体的弹性坐标有关的对称刚度矩阵， 和 分别是相对于刚性和弹性坐标的外力， 和分别是关于刚性和弹性坐标的二次速度矢量，其中Gr和Gf表示与刚性和弹性坐标相关i i i的g(q)的雅可比矩阵，Ω表示剩余方程的指数，I(Ω ,:)表示提取矩阵的Ω行。
4.根据权利要求3所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S2中柔性多体系统的运动方程，简化写为：式中，(·)Ω表示柔性体ith采样后的矩阵或向量， 表示柔性多体系统的质量矩阵， 表示柔性多体系统的刚度矩阵， 表示柔性多体系统的约束雅可比矩阵转置形式， 表示柔性多体系统的广义外力， 表示柔性多体系统的二次速度矢量。
5.根据权利要求4所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S3中每个时间步长下柔性多体系统的动态响应采用一阶后向欧拉法。
6.根据权利要求5所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S3中每个时间步长下柔性多体系统的动态响应，非线性映射算子定义为：m×n m×l
式中，Φ∈R 是一个测量矩阵，Φ(·)∈R 表示非线性映射算子， 表示t时刻的弹性坐标和模态坐标，h表示步长，M(pt)表示t时刻质量矩阵， 表示t时刻柔性多体T系统采样后的矩阵或向量，f(pt)表示t时刻作用于柔性体的广义外力，G (pt)表示t时刻约束雅可比矩阵的转置， 表示t时刻拉格朗日乘子，g(pt)表示t时刻约束方程，且a表示系数向量。
7.根据权利要求6所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S3中每个时间步长动态响应的l1范数优化问题，定义为：约束
式中， 表示t时刻柔性多体系统的模态坐标，||·||1表示1范数，s表示柔性多体系统模态坐标的稀疏度，则 的稀疏性可表示为:i
8.根据权利要求7所述的基于自适应模态的柔性多体系统动态响应计算方法，其特征在于，步骤S4.1中所述初始参数包括p0, Time,h，其中p0表示初始时刻的坐标， 表示初始时刻的速度，Time表示模拟时间，h表示时间步长。
9.一种基于自适应模态的柔性多体系统动态响应计算系统，其特征在于，所述系统基于权利要求1至8任一项所述的基于自适应模态的柔性多体系统动态响应计算方法，包括：模态模块，所述模态模块将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
抽样模块，所述抽样模块设计柔性多体系统的运动方程的抽样矩阵；
范数优化模块，所述范数优化模块求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
贪婪高斯‑牛顿求解模块，所述贪婪高斯‑牛顿求解模块利用贪婪高斯‑牛顿算法求解l1范数优化问题，得到柔性多体系统动态响应，具体为：S4.1：读取每个柔性体的有限元信息，计算每个柔性体的惯性形状积分模态，设置初始参数，所述惯性形状积分模态包括mff和Kff；
S4.2：求解 得到初始时刻的加速度 和拉格朗日乘子λ0，式中，M表示柔性多体系统的质量矩阵，G表示柔性多体系统的约束雅可比矩阵，表示柔性多体系统的加速度，λ表示柔性多体系统的拉格朗日乘子， 表示作用于柔性多体系统的广义外力， 表示与加速度有关的向量；
S4.3：对每个柔性体进行测量，得到其待定运动方程Φ(c)；
S4.4：利用贪婪高斯‑牛顿算法求解 约束 得到pt, 分别表示在t时刻的坐标，速度，加速度和拉格朗日乘子；
S4.5：判断t是否大于Time，若是，则结束计算，若不是，令t＝t+h，返回步骤S4.3。
基于自适应模态的柔性多体系统动态响应计算方法和系统
技术领域
[0001] 本发明涉及柔性多体系统技术领域，更具体地，涉及一种基于自适应模态的柔性多体系统动态响应计算方法和系统。
背景技术
[0002] 柔性多体系统(FMS)是由刚性和变形构件以某种方式连接而成的系统。利用有限元法将柔性体离散为有限自由度。对于复杂的组件，利用有限元法将产生数千个自由度。这将给柔性多体动态响应分析的计算效率带来负担。因此，需要采用模型降阶技术来降低自由度。模态截断是一种常用的约简技术。它通过选取低频模态形状将弹性坐标转换为模态坐标。应用此方法后，柔性部件的自由度将显著降低。然而，模态截断主要缺点是不能直接考虑执行器的控制系统或位置，对许多低频模态集对结构和控制系统之间的相互作用没有贡献。
[0003] 公开日为2020年05月15日，公开号为CN111159636A的中国专利公开了一种基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法，首先，基于绝对节点坐标方法，建立柔性多体系统的质量矩阵、刚度矩阵和广义力列阵；其次，建立柔性多体系统的动力学方程和优化目标函数；再次，基于直接微分法或伴随变量法，建立柔性多体系统动力学的半解析灵敏度计算公式；最后，求解柔性多体系统动力学微分代数方程，获得灵敏度计算结果。该专利的计算方法效率较低。
发明内容
[0004] 本发明的首要目的是提供一种基于自适应模态的柔性多体系统动态响应计算方法，提高柔性多体动态响应分析的计算效率且突出适应不同的工况。
[0005] 本发明的次要目的是提供一种基于自适应模态的柔性多体系统动态响应计算系统。
[0006] 为解决上述技术问题，本发明的技术方案如下：
[0007] 一种基于自适应模态的柔性多体系统动态响应计算方法，包括以下步骤：
[0008] S1：将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
[0009] S2：设计柔性多体系统的运动方程的抽样矩阵；
[0010] S3：求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
[0011] S4：利用贪婪高斯‑牛顿算法求解l1范数优化问题，得到柔性多体系统动态响应。
[0012] 优选地，步骤S1中将柔性多体系统的模态坐标定义为正交振型上的稀疏系数，具体为：
[0013] 设模态坐标是稀疏的或近似稀疏的，柔性多体系统中有N个柔性体，则广义坐标q表示如下：
[0014]
[0015] 式中，i表示第i个柔性体柔性体， 为柔性体i的弹性坐标的向iT iT
量，R 为柔性体i相对于惯性系的位置向量，θ 为柔性体i相对于惯性系的角位移，I是单位i
矩阵，B 是由柔性体i的全模态形状组成， 其中 表示柔性
多体系统的弹性坐标， 表示柔性多体系统的模态坐标；
[0016] 结合拉格朗日乘子向量λ，可得：
[0017]
[0018] 式中， 表示信号， 表示基， 表示系数向量。
[0019] 优选地，步骤S2中柔性多体系统的运动方程的抽样矩阵，具体为：
[0020]
[0021] 式中， 表示相对于柔性体的位置和旋转的质量矩阵， 是耦合刚性运动和变形的质量矩阵， 表示单个柔性体与弹性坐标有关的质量矩阵， 表示与单个柔性体的弹性坐标有关的对称刚度矩阵， 和 分别是相对于刚性和弹性坐标的外力， 和分别是关于刚性和弹性坐标的二次速度矢量，其中Gr和Gf表示与刚性和弹性坐标相关i i i的g(q)的Jacobin矩阵，Ω表示剩余方程的指数，I(Ω ,:)表示提取单元矩阵的Ω行。
[0022] 优选地，步骤S2中柔性多体系统的运动方程，简化写为：
[0023]
[0024] 式中，(·)Ω表示柔性体ith采样后的矩阵或向量， 表示柔性多体系统的质量矩阵， 表示柔性多体系统的刚度矩阵，
表示柔性多体系统的约束雅可比矩阵转置形式， 表示柔
性多体系统的广义外力， 表示柔性多体系统的二次速度矢量。
[0025] 优选地，步骤S3中每个时间步长下柔性多体系统的动态响应采用一阶后向欧拉法。
[0026] 优选地，步骤S3中每个时间步长下柔性多体系统的动态响应，非线性映射算子定义为：
[0027]
[0028] 式中，Φ∈Rm×n是一个测量矩阵，Φ(·)∈Rm×l表示非线性映射算子， 表示t时刻的弹性坐标和模态坐标，h表示步长，M(pt)表示t时刻质量矩阵，(·)Ωt表示t时刻柔T性多体系统采样后的矩阵或向量，f(pt)表示t时刻作用于柔性体的广义外力，G (pt)表示t时刻约束雅可比矩阵的转置， 表示t时刻拉格朗日乘子，g(pt)表示t时刻约束方程，且a表示系数向量。
[0029] 优选地，步骤S3中每个时间步长动态响应的l1范数优化问题，定义为：
[0030] 约束
[0031] 式中， 表示t时刻柔性多体系统的模态坐标，||·||1表示1范数，s表示柔性多体系统模态坐标的稀疏度，则 的稀疏性可表示为:
[0032]
[0033] 式中，si表示柔性体i模态坐标的稀疏度。
[0034] 优选地，步骤S4中利用GGN算法求解l1范数优化问题，具体为：
[0035] S4.1：读取每个柔性体的有限元信息，计算每个柔性体的惯性形状积分模态，设置初始参数，所述惯性形状积分模态包括mff和Kff；
[0036] S4.2：求解 得到初始时刻的加速度 和拉格朗日乘子λ0，式中，M表示柔性多体系统的质量矩阵，G表示柔性多体系统的约束雅可比矩阵， 表示柔性多体系统的加速度，λ表示柔性多体系统的拉格朗日乘子， 表示作用于柔性多体系统的广义外力， 表示与加速度有关的向量；
[0037] S4.3：对每个柔性体进行测量，得到其待定运动方程Φ(c)；
[0038] S4.4：利用GGN算法求解 得到p1, 表示在t时刻的坐标，速度，加速度和拉格朗日乘子；
[0039] S4.5：判断t是否大于Time，若是，则结束计算，若不是，令t＝t+h，返回步骤S4.3。
[0040] 优选地，步骤S4.1中所述初始参数包括p0, Time,h，其中p0表示初始时刻的坐标， 表示初始时刻的速度，Time表示模拟时间，h表示时间步长。
[0041] 一种基于自适应模态的柔性多体系统动态响应计算系统，所述系统基于上述所述的基于自适应模态的柔性多体系统动态响应计算方法，包括：
[0042] 模态模块，所述模态模块将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
[0043] 抽样模块，所述抽样模块设计柔性多体系统的运动方程的抽样矩阵；
[0044] 范数优化模块，所述范数优化模块求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
[0045] GGN求解模块，所述GGN求解模块利用GGN算法求解l1范数优化问题，得到柔性多体系统动态响应。
[0046] 与现有技术相比，本发明技术方案的有益效果是：
[0047] 本发明定义了可以用模态形状稀疏表示的弹性坐标，导出了具有少量运动方程约束的宽范数坐标下的l1范数优化问题，作为每个时间步长的动态响应分析。提出了GGN算法，有效地解决了这种优化问题，使模态坐标能够自适应地选择各种工作模式。
附图说明
[0048] 图1为本发明的方法流程示意图。
[0049] 图2为实施例1中GGN算法流程图。
[0050] 图3为实施例1中曲柄滑块机构示意图。
[0051] 图4为实施例1中曲柄滑块机构点C的x位移图。
[0052] 图5为使用实施例1所述的方法与FOM之间的相对误差。
[0053] 图6为使用实施例1所述的方法计算模态坐标的值和选择示意图。
[0054] 图7为实施例2的系统示意图。
具体实施方式
[0055] 附图仅用于示例性说明，不能理解为对本专利的限制；
[0056] 为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；
[0057] 对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。
[0058] 下面结合附图和实施例对本发明的技术方案做进一步的说明。
[0059] 实施例1
[0060] 本实施例提供一种基于自适应模态的柔性多体系统动态响应计算方法，如图1，包括以下步骤：
[0061] S1：将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
[0062] S2：设计柔性多体系统的运动方程的抽样矩阵；
[0063] S3：求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
[0064] S4：利用贪婪高斯‑牛顿算法求解l1范数优化问题，得到柔性多体系统动态响应。
[0065] 步骤S1中将柔性多体系统的模态坐标定义为正交振型上的稀疏系数，具体为：
[0066] 设模态坐标是稀疏的或近似稀疏的，柔性多体系统中有N个柔性体，则广义坐标q表示如下：
[0067]
[0068] 式中，i表示第i个柔性体柔性体， 为柔性体i的弹性坐标的向iT iT
量，R 为柔性体i相对于惯性系的位置向量，θ 为柔性体i相对于惯性系的角位移，I是单位i
矩阵，B是由柔性体i的全模态形状组成， 其中 表示柔性多
体系统的弹性坐标， 表示柔性多体系统的模态坐标；
[0069] 结合拉格朗日乘子向量λ，可得：
[0070]
[0071] 式中， 表示信号， 表示基， 表示系数向量。
[0072] 步骤S2中柔性多体系统的运动方程的抽样矩阵，具体为：
[0073]
[0074] 式中， 表示相对于柔性体的位置和旋转的质量矩阵， 是耦合刚性运动和变形的质量矩阵， 表示单个柔性体与弹性坐标有关的质量矩阵， 表示与单个柔性体的弹性坐标有关的对称刚度矩阵， 和 分别是相对于刚性和弹性坐标的外力， 和分别是关于刚性和弹性坐标的二次速度矢量，其中Gr和Gf表示与刚性和弹性坐标相关i i i的g(q)的Jacobin矩阵，Ω表示剩余方程的指数，I(Ω ,:)表示提取单元矩阵的Ω行。
[0075] 步骤S2中柔性多体系统的运动方程，简化写为：
[0076]
[0077] 式中，(·)Ω表示柔性体ith采样后的矩阵或向量， 表示柔性多体系统的质量矩阵， 表示柔性多体系统的刚度矩阵，
表示柔性多体系统的约束雅可比矩阵转置形式， 表示柔
性多体系统的广义外力， 表示柔性多体系统的二次速度矢量。
[0078] 步骤S3中每个时间步长下柔性多体系统的动态响应采用一阶后向欧拉法。
[0079] 步骤S3中每个时间步长下柔性多体系统的动态响应，非线性映射算子定义为：
[0080]m×n m×l
[0081] 式中，Φ∈R 是一个测量矩阵，Φ(·)∈R 表示非线性映射算子， 表示t时刻的弹性坐标和模态坐标，h表示步长，M(pt)表示t时刻质量矩阵，(·)Ωt表示t时刻柔T性多体系统采样后的矩阵或向量，f(pt)表示t时刻作用于柔性体的广义外力，G (pt)表示t时刻约束雅可比矩阵的转置， 表示t时刻拉格朗日乘子，g(pt)表示t时刻约束方程，且a表示系数向量。
[0082] 步骤S3中每个时间步长动态响应的l1范数优化问题，定义为：
[0083] 约束
[0084] 式中， 表示t时刻柔性多体系统的模态坐标，||·||1表示1范数，s表示柔性多体系统模态坐标的稀疏度，则 的稀疏性可表示为:
[0085]
[0086] 式中，si表示柔性体i模态坐标的稀疏度。
[0087] 步骤S4中利用GGN算法求解l1范数优化问题，具体为：
[0088] S4.1：读取每个柔性体的有限元信息，计算每个柔性体的惯性形状积分模态，设置初始参数，所述惯性形状积分模态包括mff和Kff；
[0089] S4.2：求解 得到初始时刻的加速度 和拉格朗日乘子λ0，式中，M表示柔性多体系统的质量矩阵，G表示柔性多体系统的约束雅可比矩阵， 表示柔性多体系统的加速度，λ表示柔性多体系统的拉格朗日乘子， 表示作用于柔性多体系统的广义外力， 表示与加速度有关的向量；
[0090] S4.3：对每个柔性体进行测量，得到其待定运动方程Φ(c)；
[0091] S4.4：利用GGN算法求解 得到p1, 表示在t时刻的坐标，速度，加速度和拉格朗日乘子，所述利用GGN算法求解的具体过程如图2所示，在算法中 ，符号 Sk表示支持集 ，包括刚性坐 标和拉格朗日 乘子的索 引，和表示提取 的tth列， 由 列组
+
成，对应于支持集Sn。tol和ε表示公差。运算符(·) 表示伪逆；
[0092] S4.5：判断t是否大于Time，若是，则结束计算，若不是，令t＝t+h，返回步骤S4.3。
[0093] 步骤S4.1中所述初始参数包括p0, Time,h，其中p0表示初始时刻的坐标， 表示初始时刻的速度，Time表示模拟时间，h表示时间步长。
[0094] 为了进一步说明本发明的可行性和有效性，选用了一个曲柄滑块系统，具体如图3所示。在本例中，曲柄和滑块被设置为刚性，连杆被视为柔性。连杆机构由3088个低阶三角形单元组成，自由度为6060。连杆的厚度、材料密度、杨氏模量和泊松比分别设置为0.01m、2700kg/m3、70Gpa和0.3。链接的参考条件设置为固定。曲柄和滑块的z轴质量和惯性矩分别设置为0.475kg，2.8169E‑4kg.m2和0.290kg，2.31889E‑4kg.m2。总仿真时间和时间步长分别设置为0.2s和1E‑4s，联动稀疏度设置为120，采样次数设置为400。并将扭矩M＝15Nm施加到曲柄上。在初始配置中，曲柄和连杆沿x方向放置。
[0095] 分别计算FOM算法和GGN算法。C点的x位移如图4所示。FOM算法的结果和GGN的算法高度接近。如图5所示，FOM算法与GGN算法之间的相对误差。其相对误差的大小为1E‑4，这也意味着该方法的结果是准确的。如图6所示，所GGN算法自适应地选择了120个模态坐标的最大绝对值。由计算结果可知，FOM算法的计算时间为37725秒，而GGN算法的计算时间仅为18786秒，使用该方法计算时间几乎减少了50％。
[0096] 实施例2
[0097] 本实施例提供一种基于自适应模态的柔性多体系统动态响应计算系统，如图7，所述系统基于实施例1所述的基于自适应模态的柔性多体系统动态响应计算方法，包括：
[0098] 模态模块，所述模态模块将柔性多体系统的模态坐标定义为正交振型上的稀疏系数；
[0099] 抽样模块，所述抽样模块设计柔性多体系统的运动方程的抽样矩阵；
[0100] 范数优化模块，所述范数优化模块求解每个时间步长下柔性多体系统的动态响应，得到每个时间步长动态响应的l1范数优化问题；
[0101] GGN求解模块，所述GGN求解模块利用GGN算法求解l1范数优化问题，得到柔性多体系统动态响应。
[0102] 相同或相似的标号对应相同或相似的部件；
[0103] 附图中描述位置关系的用语仅用于示例性说明，不能理解为对本专利的限制；
[0104] 显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。
