CN114347029A - 一种用于气动软体机器人快速模拟的模型降阶方法 - Google Patents

Abstract

本发明属于模型降阶领域，提供一种用于气动软体机器人快速模拟的模型降阶方法。首先，根据已有的气动机器人模型，通过剖分网格、建立应变能约束构建其动力学方程；其次，基于动力学方程，根据初始状态计算线性模态以及模态导数，并由此构建非线性模态；再次，根据非线性模态建立转换矩阵，对动力学方程进行降阶；最后，通过数值积分方法，求解气动机器人变形。本发明采用基于位置动力学方法建立气动机器人仿真框架，通过位移对广义坐标求二阶导获得模态导数，利用线性模态和模态导数形成降阶矩阵，目的在于建立一种快速、稳定的气动机器人动力学仿真系统，为气动机器人的设计提供参考。

Description

一种用于气动软体机器人快速模拟的模型降阶方法

技术领域

本发明属于模型降阶领域，涉及一种气动机器人动力学仿真的模型降阶方法。

背景技术

近年来，随着人们对机器人的研究越来越趋于轻量化、简易化和快速化，充气软体机器人作为一种新型软体机器人应用而生并且逐渐成为研究热点。这种机器人利用柔软材料制作，通过充气实现驱动，能够适应各种非结构化环境，与人类的交互也更安全。由于软体机器人在实际测试前需要大量物理仿真验证所涉及的构型是否可以完成预期规定的动作，然而每一次物理仿真都十分耗时，这使得气动软体机器人的设计和仿真仍面临挑战。

模型降阶是解决气动软体机器人快速仿真问题的关键环节。目前常用的模型降阶方法诸如特征正交分解法(POD)、动力模态分解法(DMD)等，这些方法大多都是针对线性系统，特别是线性定常系统。然而，气动软体机器人在充气力的作用下大多产生非线性变形。目前对于非线性模型降阶常用的方法是利用全阶模型的信息通过统计技术提取具有代表性的低维空间，即进行一个未降阶的完整仿真，从总选取某些时刻系统的变形u₁，u₂，…，u_n，由这些变形构成一个空间A，通过对A＝[u₁,u₂,…u_n]进行奇异值分解(SVD)得到子空间U，根据选取的降阶最大阶数r选取U的前r列作为转换矩阵。实际中，对全局模型进行一次高精度的仿真往往巨有挑战性并且耗时。

发明内容

本发明提出一种气动机器人动力学仿真的模型降阶方法。该方法采用基于位置动力学方法建立气动机器人模型的动力学方程，并且利用模态导数构成非线性模态，以解决气动机器人在运动过程中存在的非线性模型降阶问题，目的在于提供一套快速、稳定的气动机器人动力学仿真系统，以为气动机器人的设计和测试提供一个快速仿真的参考。

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

一种气动机器人模型降阶的方法，该方法首先，根据已有的气动机器人模型，通过剖分网格、建立应变能约束构建其动力学方程；其次，基于动力学方程，根据初始状态计算线性模态以及模态导数，并由此构建非线性模态；再次，根据非线性模态建立转换矩阵，对动力学方程进行降阶；最后，通过数值积分方法，求解气动机器人变形。包括以下步骤：

第一步：建立气动机器人基于位置动力学方法的动力学方程

基于位置动力学方法通过建立各种约束建立节点之间的关系，动力学方程由n个节点和Q个约束组成：

其中，

代表系统内力，是内能E对位移u的梯度。能量函数E(u)＝1/2C(u)^Tα^-1C(u)则由约束函数C＝[C₁(u)C₂(u)...C_Q(u)]^T构成。α是对角柔度阵，即刚度阵逆矩阵。u＝[u₁u₂...u_n]^T是节点位移列向量,

是节点加速度，M是对角矩阵,f_ext是外力的列向量。

基于位置动力学方法本质是根据约束函数将每个点投影到合适的位置，从当前位置到下一个位置的变化可以表示为约束函数的梯度。传统基于位置动力学方法的约束通常采用距离约束、二面角约束、体积约束等几何约束，这些约束构造简单，计算结果满足视觉合理性，但是不包含物理量，无法与实际材料相结合，计算结果无法反映变形后物体的各种物理特性。本发明在此基础上采用连续介质力学中的StVD超弹性材料构建内能从而建立应变约束，引入实际工程物理量，实现对模型定性和定量的分析。例如针对哈佛微型机器人实验室开发了一种软体抓手，这种软体机器人与传统的机械型机器人相比更适用于多变的复杂环境，以及抓取易碎易损坏的物体，如拿纸杯、茶壶等等。本发明利用这种软体材料的弹性模量和泊松比建立软体抓手的应变约束，从而构建该模型的动力学方程，通过仿真，不仅可以模拟此软体机器人在不同充气力作用下的变形，还可以定量分析各种充气力作用下的物理特性，例如应力、应变分布，以及受力情况，这是传统基于位置动力学无法实现的。引入的应变约束表示如下：

其中，V是初始构型下四面体单元的体积；λ和μ是拉梅常数，它们均由材料的弹性模量和泊松比决定；E是系统内能；刚度阵K_c与基于位置动力学方程中的α^-1等效，属于单元刚度阵；约束C＝[ε_xxε_yyε_zzε_xyε_xzε_xy]^T由应变ε组成。虽然基于位置动力学方法计算快且稳定性强，但是对于充气结构，为了能够精确模拟它的变形，需要对其进行大量网格划分，这导致约束增多，计算成本增加。因此，本发明对基于位置动力学方法构建的动力学方程进行降阶，在保证一定精度的同时，大幅提升计算效率，可以实现实时计算。

第二步：建立气动机器人模型降阶动力学方程

在降阶模型中，位移向量u可以表示为u＝Uq，其中

是全局位移空间下r维子空间的基矩阵，不随时间发生变化，

为广义坐标向量。对于U，一般选择一组正交基，使得U^TMU＝I_r,其中

为单位矩阵。将u＝Uq带入公式(1),并左乘U^T，就会得到降阶动力学方程，它的变量是q，但也同时描述了全局变形u(t)＝Uq(t)。

其中，

分别是降阶后的内力和外力：

f_ext(q)＝U^Tf_ext(Uq) (5)

同样的，整个降阶系统的刚度阵可以表示为：

对于StVK材料，它的应变能函数由格林应变构成，因此应变能可以表达为位移变量u的四阶多元多项式函数。内力f_int(u)作为应变能对于位移变量u的导数，每一个节点的内力都可以表示为一个三阶多元多项式函数。因此，降阶后的内力

中的每一项都是关于q的三次多项式：

其中，Pⁱ，Q^ij，S^ijk属于常向量系数。降阶后的刚度阵

作为

的雅可比矩阵，它的每一项都是关于q的二次多项式，

的第l列可以表示为：

在材料参数和降阶矩阵U确定后多项式系数Pⁱ，Q^ij，S^ijk可以离线计算，因此在气动机器人降阶后的仿真过程中节省了大量的计算时间。

第三步：计算降阶矩阵U

利用一组线性模态组成降阶矩阵实现模态截断，对于线性系统而言是一种常用的模型降阶方法。然而，对于非线性系统降阶，这种方法通常不是很精确。为了获得适合非线性系统的降阶基，本发明采用多体动力学领域中应用的模态导数方法形成非线性模态，利用非线性模态组成的降阶基从而描述气动软体机器人系统变形。

(a)线性模态求解

对于公式(1),对应的对称线性特征值问题可以表示为:

其中，ω_i是固有频率，Ф_i是第i阶模态，Ф＝[Ф_i,…Ф_3n]是由模态组成的矩阵，

是特征值，Λ＝diag([λ₁,…λ_3n])是特征值矩阵，一般地，特征值按照从小到到大排序，即：λ₁<λ₂<…<λ_3n。模态归一化：

因此固有模态为

空间的一组正交基，所以对于位移u可以表示为u＝Фq，q被称为广义坐标。把它带入静力学平衡方程Ku＝f_ext中，并且左右同乘Ф^T可得：

代入模态正交归一化公式可得外力与广义坐标q之间的关系为：

f_ext＝MΦΛq (12)

变换广义坐标，即当q＝p时，坐标转换公式为u＝Фp右端项也可以表示为f_ext＝MФΛp。

对于瞬态问题，它的基本方程为：

在简谐波外力f_ext(t)＝f_extcos(ωt)作用下,位移表示为u(t)＝ucos(ωt)。当u＝Фq时，把公式(10)带入动力学方程公式(13)，即可得：

-ω²Iq+Λq＝Φ^Tf_ext＝Λp (14)

所以：

这一结果表明当外部荷载频率固定在某一频率时，其可简化为接近外部荷载频率的多阶模态振型的和，而不是使用计算复杂的全阶模态。在实际工程中，外部荷载往往处于低频，因此可将模态截断为低阶模态，即取前r阶模态组成降阶矩阵，这种降阶方法又称一阶线性模态降阶。

(b)非线性模态求解

对于大变形问题，系统的变形往往是非线性的，因此对系统的位移u(q)近似取二阶麦克劳林展开：

显然，u(0)＝0，

所以上式可以简化为：

为了求解

利用公式公式(12)对广义坐标q求偏导可得：

公式(18)再次对广义坐标q求导可得：

已知刚度阵

刚度阵的海森阵

可得：

KΨ＝-(H:Φ)Φ (20)

被称为模态导数，因此公式(16)可以写为：

由第二步可知，一般取前r阶线性模态构成线性子空间，因此由上式可以看出Ф_i和Ψ_ij构成了一个r+r(r+1)维非线性运动空间，考虑到它的维数会影响计算效率，因此根据对应的线性模态的特征值来缩放模态导数：

模态导数缩放让低频模态和它们的导数具有更大的权重，防止它们可能被高频模态和模态导数掩盖。之后通过对A进行奇异值分解得到非线性子空间U，其中U的第i列就是这个系统的第i阶非线性模态。把U带入第二步，对动力学系统进行Newmark积分就可得到气动机器人的变形。

本发明的有益效果为：

(1)本发明在传统基于位置动力学方法的基础上，采用连续介质力学中StVD材料的应变能函数建立应变约束，构建气动机器人的动力学方程，简化了模型，与只采用几何约束的传统基于位置动力学方法相比，引入实际物理参数，可以直观反映应力、应变等物理特性，对气动机器人的变形做定性和定量的分析。

(2)本发明在构建基于位置动力学方法模拟气动机器人系统的基础上，对其进行了模型降阶处理。该方法在保证一定精度的情况下有更高的计算效率，可以实现实时仿真计算，为气动机器人的设计仿真方面提供了参考。

附图说明

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

图2是气动机器人动力学模型。

图3是气动机器人前四阶模态；其中，图(a)为前四阶的线性模态；图(b)前四阶线性模态对应的模态导数；图(c)为前四阶的非线性模态；

图4是气动机器人模型降阶仿真与全局仿真、实验照片的对比图；其中，图(a)为气动机器人实验图；图(b)气动机器人在有限元方法下全局仿真图；图(c)为气动机器人模型降阶仿真图

具体实施方式

以下结合具体实施例对本发明做进一步说明。

结合图1，以哈佛大学Whitesides实验室设计的PneuNet软体驱动器为例，对本发明所涉及的方法进行验证。PneuNet软体驱动器是由小通道的气动网络(“气动网”)充气驱动，利用简单的控制制造就可以实现复杂的运动。为了节约设计成本，PneuNet软体驱动器在进行实验前都会先进行大量的物理仿真从而验证可行性。传统的物理仿真基于全局模型，耗时长，因此本专利采用模型降阶方法进行降维计算，提高计算效率。具体实施方案如下：

(1)利用三维建模软件Solidworks对PneuNet软体驱动器进行建模如图2所示，材料的弹性模量和泊松比分别为1e7Pa和0.45，在构型上对它采取斜痕设计。该气动机器人长度w＝259.5mm，高度h＝40mm，宽度d＝40mm。对于内部每一个充气腔体，腔体右壁l₁和左壁l₃厚度相同，l₁＝l₃＝5mm，腔体宽度l₂＝20mm，相邻两个腔体之间的距离l₄＝5mm。腔体上层厚度r₁＝3mm，下层厚度r₅＝7mm，腔体与腔体连接处的厚度r₃＝3mm，腔体高度r₂＝15mm。气动机器人一端固定，腔体内部通过施加压力模拟充气过程。

(2)利用划分网格软件Hypermesh选取四面体单元对上述模型进行剖分，提取网格信息诸如：节点坐标、单元编号，利用公式(2)构建应变能约束，同时采取基于位置动力学方法构建气动软体机器人的动力学方程。

(3)对PneuNet软体驱动器模型进行降阶。根据上述建立的动力学方程计算该模型初始状态下的质量阵M和刚度阵K，利用广义特征值问题求解方法计算线性模态Ф，并根据公式(20)计算该模型的模态导数Ψ。利用线性模态和模态导数构成运动子空间A，对A进行奇异值分解获得非线性模态，具体如图3所示。以取前4阶线性模态为例，模态导数由于存在对称性形所以共有10个不同的模态导数。这4个线性模态和10个模态导数共同构成了气动机器人运动子空间，对其进行SVD分解便可得到非线性模态。利用非线性模态构成降阶矩阵U，带入公式(3)、(4)和(5)中，完成动力学方程降阶。根据公式(7)和(8)计算降阶后内力与刚度阵中每一项多项式的常系数，并且存储数据。

(4)将上一步计算的多项式的常系数直接带入动力学方程(3)中，给定气压值大小，通过数值积分方法求解PneuNet软体驱动器模型在充气力作用下的变形。图4是PneuNet软体驱动器模型在实验、有限元全局模型仿真和本专利选取前20阶非线性模态进行降阶仿真结果的对比图。由图可知本专利计算PneuNet软体驱动器降阶模型各个时刻的变形与实验结果相近，可以为PneuNet软体驱动器的设计提供一种数值仿真参考。对于计算时间，本专利完成数值仿真时间仅需34.45s，但是利用有限元全局仿真则需要1200s左右，因此本发明在保证一定精度的情况下极大的提高了计算效率。

以上所述实施例仅表达本发明的实施方式，但并不能因此而理解为对本发明专利的范围的限制。应当指出，对于本领域的技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些均属于本发明的保护范围。

Claims (1)

1.一种用于气动软体机器人快速模拟的模型降阶方法，其特征在于，首先根据已有的气动机器人模型，通过剖分网格、建立应变能约束构建其动力学方程；其次，基于动力学方程，根据初始状态计算线性模态以及模态导数，并由此构建非线性模态；再次，根据非线性模态建立转换矩阵，对动力学方程进行降阶；最后，通过数值积分方法，求解气动机器人变形；包括以下步骤：

第一步：建立气动机器人基于位置动力学方法的动力学方程

基于位置动力学方法通过建立各种约束建立节点之间的关系，动力学方程由n个节点和Q个约束组成：

其中，f_int＝▽_uE(u)代表系统内力，是内能E对位移u的梯度；能量函数E(u)＝1/2C(u)^Tα^-1C(u)则由约束函数C＝[C₁(u)C₂(u)...C_Q(u)]^T构成；α是对角柔度阵，即刚度阵逆矩阵；u＝[u₁u₂...u_n]^T是节点位移列向量,

是节点加速度，M是对角矩阵,f_ext是外力的列向量；

在位置动力学方法基础上采用连续介质力学中的StVD超弹性材料构建内能从而建立应变约束，引入实际工程物理量，实现对模型定性和定量的分析；引入的应变约束表示如下：

其中，V是初始构型下四面体单元的体积；λ和μ是拉梅常数，它们均由材料的弹性模量和泊松比决定；E是系统内能；刚度阵K_c与基于位置动力学方程中的α^-1等效，属于单元刚度阵；约束C＝[ε_xxε_yyε_zzε_xyε_xzε_xy]^T由应变ε组成；

对基于位置动力学方法构建的动力学方程进行降阶，在保证精度的同时，提升计算效率，可以实现实时计算；

第二步：建立气动机器人模型降阶动力学方程

在降阶模型中，位移向量u可以表示为u＝Uq，其中

是全局位移空间下r维子空间的基矩阵，不随时间发生变化，

为广义坐标向量；对于U，选择一组正交基，使得U^TMU＝I_r,其中

为单位矩阵；将u＝Uq带入公式(1)，并左乘U^T，得到降阶动力学方程，它的变量是q，但也同时描述全局变形u(t)＝Uq(t)；

其中，

分别是降阶后的内力和外力：

f_ext(q)＝U^Tf_ext(Uq) (5)

同样的，整个降阶系统的刚度阵表示为：

对于StVK材料，它的应变能函数由格林应变构成，因此应变能可以表达为位移变量u的四阶多元多项式函数；内力f_int(u)作为应变能对于位移变量u的导数，每一个节点的内力都可以表示为一个三阶多元多项式函数；因此，降阶后的内力

中的每一项都是关于q的三次多项式：

其中，Pⁱ，Q^ij，S^ijk属于常向量系数；降阶后的刚度阵

作为

的雅可比矩阵，它的每一项都是关于q的二次多项式，

的第l列可以表示为：

在材料参数和降阶矩阵U确定后多项式系数Pⁱ，Q^ij，S^ijk可以离线计算；

第三步：计算降阶矩阵U

利用一组线性模态组成降阶矩阵实现模态截断，为了获得适合非线性系统的降阶基，采用多体动力学领域中应用的模态导数方法形成非线性模态，利用非线性模态组成的降阶基从而描述气动软体机器人系统变形；

(a)线性模态求解

对于公式(1)，对应的对称线性特征值问题表示为：

其中，ω_i是固有频率，Ф_i是第i阶模态，Ф＝[Ф_i,…Ф_3n]是由模态组成的矩阵，

是特征值，Λ＝diag([λ₁,…λ_3n])是特征值矩阵，特征值按照从小到大排序，即：λ₁<λ₂<…<λ_3n；模态归一化：

因此固有模态为

空间的一组正交基，所以对于位移u可以表示为u＝Фq，q被称为广义坐标；

外力与广义坐标q之间的关系为：

f_ext＝MΦΛq (11)

对于瞬态问题，基本方程为：

在简谐波外力f_ext(t)＝f_extcos(ωt)作用下，位移表示为u(t)＝ucos(ωt)；当u＝Фq时，把公式(10)带入动力学方程公式(12)，即可得：

-ω²Iq+Λq＝Φ^Tf_ext＝Λp (13)

所以：

在实际工程中，外部荷载处于低频，因此可将模态截断为低阶模态，即取前r阶模态组成降阶矩阵；

(b)非线性模态求解

对于大变形问题，对系统的位移u(q)近似取二阶麦克劳林展开：

公式(15)可进一步表示为：

由第二步可知，根据对应的线性模态的特征值来缩放模态导数：

模态导数缩放让低频模态和它们的导数具有更大的权重，防止它们可能被高频模态和模态导数掩盖；之后通过对A进行奇异值分解得到非线性子空间U，其中U的第i列就是这个系统的第i阶非线性模态；把U带入第二步，对动力学系统进行Newmark积分就可得到气动机器人的变形。

Images