对超弹性材料应用损伤演化定律

今天的特邀博主是来自Lightness by Design 的Björn Fallqvis，他分两部分讨论了超弹性材料的分析方法，首先通过搭接剪切试验估计超弹性材料参数，随后建立了一个损伤模型来验证材料的软化行为。本文为第二部分内容，您可以点击此处阅读第一部分内容。

损伤模型的数学背景

为了建立超弹性材料的损伤模型，我们用一个损伤变量来修改应变能密度函数 \omega;也就是说，应变能密度现在被定义为 \omega \psi。根据定义，\omega = 1 表示未损伤的状态；\omega = 0 时，材料完全断裂。可以注意到，这与习惯上统一减去损伤变量的定义不同，但却是等效的。演化定律也可以写成这样的形式，尽管符号发生了变化。

根据之前定义的演化定律(参考文献1)，我们使用下列方程对损伤演化的速率形式进行模拟:

其中，\beta和 b 是材料参数。

该方程的物理基础是阿累尼乌斯定律(Arrhenius’ law)，该定律指出反应速率(在本文中是指化学交联的断裂)是一个指数，表述为：

其中，k 是反应速率常数，A 是每个化学反应的指数前因子，R 是通用气体常数，T 为绝对温度，E_b 为反应活化能(交联势阱的深度)。

如果要考虑由于施加变形而产生的额外能量，上式则可以改写为：

反应速率 k 是交联键断裂的速率，假设损伤率与它成正比。此外，我们推测，损伤率与材料的当前状态 \omega 成正比，材料状态指材料在某时刻能够断裂的交联键数目。以此我们推导出在本小节中定义的第一个演化定律，如图 1 所示(参考文献 2)。参数 b 可以看作是用于破坏交联键的应变能（单位浓度）的一部分。

图1. 交联分子的结合状态。为了在加载状态下逃离势阱，断裂键所需的能量为 E_b-b\psi。图片来自参考文献 2。

内部耗散 D_{int} 必须大于或等于零才能在热力学上有效，根据下列方程：

可知，此演化定律是有效的，但有以下限制：\beta \geq 0。

在 COMSOL Multiphysics^® 中模拟损伤模型

由于损伤演化定律是一阶微分方程，因此可以在 COMSOL Multiphysics^® 软件中通过插入黏性域的常微分方程对象（如图 2 所示）轻松定义。将初始值 \omega 设置为 1 (无损伤)。参考温度设置为 20°C (293 K)。

图2.损伤演化的域 ODE 定义。

如上一篇博文中所述，超弹性材料定义中使用了损伤变量“omeg”，并且 \psi_{adh} 被定义为解变量，如图 3 所示。

图3. ODE 定义中使用的应变能密度函数的定义。

根据之前的描述进行优化设置，不同的是现在最小二乘目标选择了更多的点进行优化。同时，添加了两个新的损伤模型参数作为控制变量，如图4所示。

图4.设置优化，包括损伤演变。

在初始分析中，我们选择使用 2、6、10 和 13.5 mm 处的点。在第二次分析中，我们使用了所有点，将来自初始分析中的优化值作为起点。呈现的参数值(和结果)取自第二次分析的所有数据点。

结果：优化有损伤的超弹性模型的材料参数

优化参数和力-位移曲线

损伤模型的初始和优化材料参数，见表 1。

该分析的力-位移曲线如图 5 所示。

图5.有损伤的超弹性模型的力-位移曲线。

尽管实验数据和模型预测值之间存在一些差异，但在较大的应变下存在明显的软化，并且损伤演变描述了材料的软化行为。参数 b 是每个交联浓度的分数，并且对于每摩尔浓度 M，大约有 41% 的总应变能用于断裂交联键。

应变率依赖性

根据演变速率形式，演化定律捕捉到了应变率效应。为了比较，前面的分析以 130 mm/min 的末端位移速率运行，如图 6 所示。

图6. 变形率的影响。

较高的变形率意味着损伤发生的时间和机率较少。这也适用于阿累尼乌斯定律，它给出了每单位时间内发生化学反应的概率，这里是指交联分离。

滞后和材料循环

下面，我们使用表 1 中定义的材料参数进行分析。

在软材料中观察到的一个常见现象是滞后。即由于材料中的能量耗散，加载和卸载路径不重合，表现为应力-应变曲线在每个周期中的变化，并且通常在一定程度的软化(或硬化)后趋于稳定。本文提出的损伤模型捕获了滞后现象，如果我们在分析中还包括卸载，则滞后很明显。见图 7。

图7. 滞后曲线。卸载与加载以相同的速率发生。

为了进一步说明，我们将指定的 7.5 mm 末端位移指定为正弦波，频率为 0.25/s (即 4 秒为一个完整周期)。5 个初始周期的力-位移曲线如图 8 所示。

图8. 5 个初始周期的滞后曲线。

可以观察到明显的软化，并且加载-卸载路径不重合。根据损伤率曲线的定义方程(这里即指原始形式)，最后滞后曲线没有稳定。如果我们施加一个相同频率的 6000 N 正弦力，也可以看到这一点。在图 9 中，末端位移不断增加。

图9.等幅正弦力的位移。

在图 10 中，材料的位移出现在时间 1、13、25 和 37s。

图10. 循环过程中 1、13、25 和 37s 的位移。

如果希望在循环后实现材料曲线的稳定，我们可以通过使用指数衰减函数修改损伤率方程来实现：

其中，\omega_t 是材料损伤状态的阈值，低于该阈值不应发生损伤。

系数 \omega_{t,rate} 控制损伤率衰减到阈值以下的速度。这种形式是需要的，因为简单地在阈值以下将速率设置为零会导致不连续收敛问题。使用与图 8 相同的分析，但是 \omega_t=0.6，\omega_{t,rate}=1000，循环滞后曲线如图11所示。

图11. 5 个初始循环的滞后曲线，其中 \omega_t=0.6，\omega_{t,rate}=1000。

现在可以在第二个循环中看到材料响应明显稳定。

所述损伤演化定律具有三个参数(\omega_{t,rate}
只能设置为 1000)，并包含以下效应：

• 滞后性
• 单调加载期间的材料软化/破裂
• 循环材料软化
• 应变率相关性
• 蠕变性

结论和最终想法

在今天的博客中，我们展示了如何结合使用 COMSOL Multiphysics 中的固体力学 和优化 接口，通过搭接剪切试验确定材料参数。通过直接对用户定义的应变能函数进行修改，任何基于速率的损伤定律都可以用类似的方式轻松实现。

我们还看到，用于优化的点数以及我们使用的点会影响结果。在初始值相同的情况下，对所有数据点进行优化会导致曲线拟合不佳。如前所述，每隔一个点优化的最终值被用作第二次分析中的初始值。记录力最高的数据点也包含在第二次分析中，用于略微改善拟合。对于均匀分布的点更容易优化，因此如果曲线的最小二乘拟合存在问题，可以尝试选择较少的大致均匀分布的点来获得良好的初始参数值。

由于操作简单，本文演示的损伤演化法非常具有吸引力。它建立了施加的变形(应变)率和应力，循环和单调软化以及蠕变效应之间的关系。

材料响应曲线的稳定性已被证明可以通过包括损伤演变的阈值来复制。这可以被认为是一个渗透阈值，并且交联网络已经被破坏到这样一种程度，以至于聚合物链之间没有足够的连接性来完全转移交联之间的机械载荷并破坏它们。其他变形模式控制宏观特性，例如长丝缠结和熵刚度。包含阈值主要被用在材料将被循环的应用中。对于单调加载，材料将在足够大的变形下断裂(例如本博客系列中的测试模型)，两种机制之间的区别是模糊的，损伤定律可用于现象学意义。

另一个阈值定义可以是应变能密度，尽管这并不容易被用于微观结构的物理表示，并且选择阈值并不简单。

对于某些类别的材料，例如生物组织，我们还希望实现可逆的黏弹性行为。这可以通过用将超弹性参数中的等容第一不变量替换为弹性对应式，并为黏性变形指定演化定律(参考文献2)。最终，我们得到一种具有不同时间尺度的可逆和不可逆变形的材料参数。

由于演化速率表达式中存在温度变量，热效应可以包含在损伤模型中。但这主要是恒温下的蠕变应用中需要考虑的，因为材料参数也可以预期取决于此。

关于作者

Björn Fallqvist 是 Lightness by Design 公司的一名顾问，从事基于数值分析的产品开发工作。他于 2016 年获皇家理工学院博士学位，致力于开发用于捕捉生物细胞机械行为的本构模型。主要研究兴趣是材料表征领域，专长是和使用各种材料模型捕捉物理现象。

参考文献

1. B. Fallqvist, M. Kroon, “Constitutive modelling of composite biopolymer networks”, Stockholm: Journal of Theoretical Biology, vol. 395, 2016.
2. B. Fallqvist, M. Kroon, “A chemo-mechanical constitutive model for cross-linked actin networks and a theoretical assessment of their viscoelastic behaviour”, 2, Stockholm: Biomechanics and Modeling in Mechanobiology, vol. 12, 2013.