COMSOL-based simulation of pneumatic soft actuators

Part of this work has been submitted.

Copyright @典同学

气动软体机器人是近年来机械方向的研究热点，本文包含了对一类Origami inspired气动软体人工肌肉OPAM的仿真。

以下内容禁止转载

We created a new kind of pneumatic soft actuators which is called OPAM（Origami Pneumatic Artificial Muscles）.

FEM SIMULATION

1.Hyperelastic Theory

Most silicone-based elastomers are considered to be hyperelastic. That is, materials that are capable of experiencing large deformations under small loads and then return to their original shape without any significant plastic deformation once that load is removed. Hyperelastic materials, such as solid rubber, are close to ideally elastic. When deformed at a constant temperature or adiabatically, stress is solely a function of the current strain. It is independent of the history or rate of loading The stress-strain behavior of a hyperelastic material is very nonlinear, meaning that the material’s elastic modulus is not enough to characterize its elastic behavior (Figure 1). A constitutive mathematical model is, therefore, necessary to represent the real behavior of a hyperelastic material. ^[1]

Polymers are materials of very high molecular weight and they have considerable commercial importance. The polymer PDMS can be classified as a hyperelastic material, in which the strains may be large. So does the Ecoflex.

非线性弹性力学（机械本科鶸，这块没学过，若有错误，欢迎指正）

PDMS/Ecoflex作为超弹性材料，它们的应力应变曲线应该由应变能密度函数导出。

当一种材料由于外部负载而变形时，它倾向于在其整个体积内部内部存储能量。 该内部能量与材料应变有关，被称为应变能。 考虑处于拉伸状态的物体，其中体积元素呈单向应力状态。

一旦发生位移，该应力在元件的顶面和底面上产生力。 这些力的差值被定义为

其中 是单向应力，而Δ 是长度为Δ 的单元经过垂直位移后，单元的横截面积与应力垂直的变化。 功由力与相关位移的乘积定义。 由于当达到此位移时，力从零均匀增加到其最终大小，因此，作用在元件上的力功等于平均力大小与位移的乘积。 在没有能量损失的假设下，外部功变得等于内部能量，也称为内部能量。 因此，应变能可以定义为

其中 是应变，而 ∆ 是材料微元的垂直位移。 可以通过公式[2]来规范化单位体积材料的应变能。 这称为应变能密度。

现在考虑相同的体积元素，现在要求微元处于三向应力状态。 可以得到上述公式，以使每个主应力占总应变能密度的一部分。

对于线弹性材料的行为，适用胡克定律，并且可以根据材料的杨氏模量来重写这些公式。 但是，像PDMS这样的橡胶状材料的应力-应变行为是高度非线性的，因此，该材料的弹性模量不足以表征其在应对大变形时的行为。 需要一个本构模型来模拟超弹性材料和相关结构在高应变下的真实行为]。 可以从应变能密度函数（W）得出合适的本构模型，该函数表示原始几何结构的单位体积中存储在材料中的能量，是该材料在该点处的应变的函数。 在数学上，此关系定义为

其中

是格林变形张量的三个应变不变量。 各个应变不变量是主要拉伸比

(principal stretch ratios)（

， = 1、2和3）的函数。

拉伸比是垂直于材料线元素的拉伸应变的量度。 它可以在自然状态或变形状态下定义为材料线最终长度和初始长度之间的比率。 也就是，

其中 是最终长度，

是原始长度， 是材料应变。 拉伸不变量是超弹性模型可以具有非零解的属性。

等式

给出的第三个常数项与弹性体积比（ ）的关系为

＝

。 引入体积比项可以为几乎不可压缩的材料开发更方便的本构模型：

假设

和

在纯体积变化下保持恒定。

等式

给出的材料的应变能密度（ ）可以写成两部分之和。 第一部分是等容项

，它表示在保持一定体积的同时使元件变形或扭曲所需的能量。 第二部分表示导致元件的体积发生变化而不发生形状变化所需的能量，称为体积项 （ ）。 因此，应变能密度可以表示为

这里

，假设不可压缩的材料， =1。

在这次仿真中，我们使用Mooney-Rivlin模型以模拟气动软执行器中PDMS材料的超弹性特性。 由Mooney-Rivlin定义的应变能密度的一般形式为

其中 和 1分别是固体的剪切模量和体积模量。 在这种情况下，

，

。 相应地，剪切模量为 = 2(

+

)。设定材料常数

和

以便满足类似橡胶的材料的近似值(

)。

还可以使用其他超弹性模型。 根据Holzapfel ^[2]在教科书中描述的实验，将Mooney-Rivlin模型与NeoHookean，Varga和Ogden模型进行了比较，用于对橡胶球囊进行充气的评估。 简化的Neo-Hookean模型和Varga模型能够再现小应变时变形的超弹性形状的真实行为。 但是，这两个模型都无法跟踪压力变化下的局部最大和最小位移值（称为极限点）。 Mooney-Rivlin和Ogden模型能够显示这些点。 奥格登（Ogden）模型可以更逼真地估算出气球的变形； 但是，它需要超出现有范围之外的其他参数。 因此，这里选择了Mooney-Rivlin模型。

改变橡胶状材料的形状比改变其体积容易得多。因此，在计算模型中通常的做法是将它们视为不可压缩的，引入弹性体积比 = 1的约束条件，如2.1.1节中所述。在此假设下，仿真过程中设定PDMS的泊松比 为0.5。然而，由于在诸如等体积模量的性质方程的分母中引入泊松比，这样会导致计算误差。比如，体积模量

为了避免被零除，泊松比设为0.49。 但是，这引入了基于位移的有限元分析由于体积固定而产生偏斜结果的可能性。 固定的体积会通过产生过于僵化的响应来影响模型。 通过确保在COMSOL超弹性工具栏中选择“几乎不可压缩的材料”选项，可以减小体积固定带来的负面效果。 在这种情况下，材料的可压缩性足够小以至于可以忽略不计。

类似过程可以参考官方的tutorial，其中研究包含了Mooney-Rivlin，Neo-Hookean，Ogden等模型

球形橡胶气球的膨胀​cn.comsol.com

可以看到官方pdf里面表现最好的是Ogden，这个修正后的Ogden模型需要的参数太多实际操作起来有点麻烦，希望有空能补充。

2. PDMS Material Properties

COMSOL’s Material Library Module includes data for PDMS. Initial bulk and shear moduli for PDMS have listed as 3.333 × 107 Pa and 6.67 × 105 Pa, respectively.

为了满足上节中所述的剪切模量关系和近似值，

。 在标准拉伸试样上进行了单独的测试，以评估原型开发中使用的PDMS材料的机械性能。一般来说，测量结果与Johnston发表的趋势一致。

^[3]

我们需要将上述两个参数输入材料特性。

3. COMSOL Simulations and Results

• 新建项目->结构力学->固体力学（solid）
• 材料->PDMS和聚甲基丙烯酸甲酯（poly(methylmethacrylate)，简称PMMA）
• 边界载荷施加到腔室的主变形壁

导入完成后如图所示。

域选择

正如上一节所述。此时选择超弹性材料为不可压缩。

约束

在 “物理场”工具栏中，单击 “边界” 并选择 “固定约束”。

其他具体操作请参考-结构力学模块COMSOL

https://cdn.comsol.com/translated-documentation/cn/5.3a/IntroductionToStructuralMechanicsModule.pdf​cdn.comsol.com

到这一步先测试一下双参数Mooney-Rivlin

没有报错，存在收敛解，证明基本配置ok，因为还没有设置！鸭力等参数，此时的解并没有什么用。

在OPAM内壁增加一个均匀的边界载荷模拟大气压差，然后开始计算。

结果：

计算自由度过大，还没跑完。

第一次求解失败，历时5h，内存比较爆炸，到这一步基本就是参数或者模型可能有问题。

第二次求解，12h，没有收敛。

修正参数后结果。

其他步骤待补充。

网格生成

如上所述，COMSOL Multiphysics主要使用有限元方法（FEM）解决任何给定的研究。因此，创建网格是必要的，COMSOL提供了详细的系统来自动和手动创建网格。网格的配置通过将模型分成几何简单形状的小元素来确定创建离散化模型的网格的分辨率。在这些元素的每一个中使用一组偏微分方程（PDE）来近似结构位移场。也就是说，模型在（ ， ， ）方向上变形了多少。分析三个坐标方向意味着每个元素都有自己的体积。理想情况下，COMSOL可以快速解决难以置信的精细网格问题，只需几秒钟的计算即可提供准确的结果。不幸的是，情况并非如此，因为制作得太细的网格可能永远不会收敛到解。再次重要的是要注意，解的准确性与创建的网格的大小直接相关。所有计算都受到有限的计算资源和时间的限制。因此，尝试求解网格大小接近或等于零的精确解是不合理的。相反，计算必须依赖于实际解的近似值。因此，至关重要的是，最小化精确解及其计算出的近似值之间的差异或误差。这称为网格细化或独立性。在本研究论文中，这是通过在从非常粗糙到非常精细的有限网格范围内执行相同的模拟来确保的。一旦计算的数据停止受网格粗度的影响，就被认为与网格无关。另外，重要的是要考虑网格本身的质量。网格质量表示网格元素的长宽比。对于薄膜几何形状，一般的经验法则是保持最小网格质量至少为0.1。对于本研究，求解的最小网格质量为0.18，平均网格质量为0.612。（或者采用默认网格其实也行）

参考

1. ^Galley, Alexandre, "Pneumatic Hyperelastic Robotic End-Effector for Grasping Soft Curved Organic Objects" (2019). Electronic Thesis and Dissertation Repository. 6392. https://ir.lib.uwo.ca/etd/6392
2. ^G. A. Holzapfel, “Some Forms of Strain-energy Functions,” in Nonlinear Solid Mechanics: A Continuum Approach For Engineering, 1st ed., John WIley & Sons, LTD, 2010, p. 238.
3. ^I. D. Johnston, D. K. McCluskey, C. K. L. Tan, and M. C. Tracey, “Mechanical characterization of bulk Sylgard 184 for microfluidics and microengineering,” J. Micromechanics Microengineering, vol. 24, no. 3, p. 035017, Mar. 2014.