基于超弹性的纤维增强材料的准静态变形 亚历山德·弗拉努斯（Aleksander Franus）和斯坦尼斯瓦夫·耶米奥洛（Stanis?aw Jemio?o）

《Materials》：Quasi-Static Deformations of Fiber-Reinforced Materials Based on Hyperelasticity Aleksander Franus and Stanis?aw Jemio?o

　　本文的重点是基于多凸性的公式化，确保了超弹性问题解的存在性。对于中等大的变形，广义Ciarlet模型将参数与初始体积模量和Lamé常数联系起来。标量内部变量使得损伤类过程的可逆描述成为可能。对于纤维增强固体，提出了具有耗散效应的各向同性模型。基于UMAT的实现支持在Abaqus中进行有限元分析。

摘要：本研究探讨了基于超弹性公式的纤维增强材料的准静态行为，该公式通过增加粘性和损伤类效应来描述材料的响应。在内部标量变量框架内开发了一个各向同性本构模型，该模型能够实现材料损伤的可逆描述，同时确保了存储能量函数的对象性、热力学可接受性和多凸性。各向同性贡献来自于广义Ciarlet模型，而各向异性部分则考虑了嵌入在粘弹性矩阵中的一族弹性纤维，通过简单的混合理论进行解释。所得到的本构方程通过UMAT子程序在Abaqus/Standard中实现，并且其速率形式与Zaremba–Jaumann目标应力速率保持一致。通过有限元模拟检验了模型的性能，包括单轴应变和简单剪切的均匀测试、松弛和蠕变问题以及类似膨胀的问题。结果证明了该模型能够捕捉应变率敏感性、蠕变、应力松弛和能量耗散，以及非均匀变形模式，同时也指出了其在表示永久变形方面的当前局限性。

1. 引言
在土木工程领域的工程解决方案中，广泛使用了表现出非弹性行为的材料。其中一些材料类型对变形速率敏感，在标准服役条件下会经历超出小应变理论准确捕捉范围的有限弹性应变。这种情况特别适用于主要用作减振器组件的弹性体材料。例如，一项工作[1]展示了在硬质和软质弹性体产品样品上进行的典型实验测试的结果。名义应力-应变图上的加载路径与卸载路径不一致，表明了能量耗散。此外，在不同位移控制变形速率下进行的加载-卸载循环产生了定量不同的结果。这种材料特性被称为粘性[2]，这是流体固有的属性，但同样也适用于固体。一个表现出较大变形的结构的例子是由织物膜构成的屋顶[3,4]。在[3,4]中分析了这种类型的结构，但未考虑耗散效应。关于材料的流变特性，即它们对应变率的敏感性（粘性）[5]和Mullins效应[6]，在大量变形固体力学文献中已经有多种表述[7,8]。Mullins效应通常通过引入一个参数来表示，该参数定义了加载过程中达到的最大弹性能量[6,9]，或者通过一个参数来指定最大变形强度[1,10]。这构成了在内部变量理论框架内制定的模型的一个特例[11,12]，该理论为符合热力学第二定律的非弹性效应的纳入提供了正式结构[13,14]。另一种方法借鉴了具有记忆的材料理论[15,16]，在这种理论中引入了对所谓响应函数的近似。其中一种近似采用积分形式[17,18]，从而导致了准线性粘弹性（QLV）模型的出现[19,20]。这类模型也已为各向同性和纤维增强材料开发[21,22]。由于传统的描述Mullins效应的模型没有考虑应变率敏感性，因此可以在QLV框架内对其进行推广[23,24]。

在软生物组织和技术复合材料的背景下，纤维增强超弹性的一个主导框架是Holzapfel–Gasser–Ogden（HGO）模型[25]，它通过与一个或多个纤维族相关的伪不变量来表示各向异性贡献。然而，应该注意的是，标准的有限元代码（包括Abaqus/Standard[26]）中实现的HGO模型是用修改后的（偏应变）不变量来表述的，而不是本文中使用的全不变量。这种区别对于体积-偏应变耦合以及在有限应变下几乎不可压缩材料的行为有非平凡的影响。一些作者提出了将HGO框架扩展到有限应变粘弹性区域的扩展[27]。本文开发的模型可以被视为一种补充公式，在这种公式中，耗散效应完全由各向同性矩阵承担，而完整未修改的不变量在整个过程中得到保持。

各向同性QLV公式已在[21]或其他文献中发展。然而，我们注意到[21]中的公式处理了各向同性材料的速率依赖性粘弹性，但没有纳入Mullins类型的损伤或应力软化效应；因此，它是本文框架的补充，而不是直接竞争。模型[23,24]采用的本构公式——分别基于主应变和现象学软化函数——在结构上与这里采用的内部变量框架不同，后者基于不变量且符合热力学第二定律[11,12]。只有在从同一实验数据集仔细识别参数后，对这些方法与本文提出的公式进行直接定量比较才有意义。因此，这样的比较将推迟到未来的工作中进行。

在目前的工作中，重点是为使用作者先前工作中描述的纤维增强材料模型的选定的边界值问题实现有限元解[28,29]。强调该模型满足以下要求：与伽利略变换相关的对象性要求以及材料在初始配置中的给定对称性[30]；通过Clausius–Duhem不等式得出的热力学可接受性限制[31]；存储能量函数的多凸性和增长条件[32]。在本构关系的背景下，对象性要求主要涉及确定应力张量和Helmholtz自由能函数的函数。采用材料描述时，必须满足以下条件：(1) 式中的符号表示相应的量。除了与运动和当前配置的观察者相关的（欧拉的）对象性要求外，还施加了材料在初始配置中对称性的（拉格朗日的）对象性要求。因此，这里的旋转张量指的是材料坐标系统的旋转。如果满足条件(2)，则称该群为材料对称群；否则，我们称材料为各向异性的。我们可以通过参数张量来描述材料的对称性，因此，作为例如对称二阶张量的各向异性不变量的函数[33,34]。如果参数张量是对称的二阶张量[35]，则对称群定义了物体在其自然状态下的各向同性材料[34,36]。

在目前的工作中，关注点限于超弹性框架内的各向同性材料的特殊情况。各向同性材料模型包含一个参数张量，其单位向量与材料在初始配置中区分出的纤维方向一致[37,38]。因此，存储能量函数可以用五个不变量来表示：(4) 从现在开始，我们限制自己使用纯机械理论，即忽略所有热效应。在平衡方程中省略了温度和熵等场。尽管如此，这种方法通过引入所谓的机械耗散[30]，仍与热力学第二定律保持一致。因此，由热力学原理得出的要求简化为满足Clausius–Duhem不等式(5)。随后，通过引入一个内部变量来确定自由能的材料时间导数：(6) 将(6)代入(5)，我们得到(7) 不等式应该对所有可接受的值成立。因此，(8) 将张量内部变量替换为标量函数g后，耗散不等式简化为(9) 除了由热力学第二定律、对象性和材料对称性得出的要求外，还应该引入额外的超弹性本构关系要求，以便运动方程的初始-边界值问题有解[39]。这里，我们主要要求存储能量函数是多凸的，并满足适当的增长条件[32,40]。由于不变量对于是对的凸函数，而不变量对于是对的凸函数，因此函数具有以下形式(10) 如果函数和对于张量是凸函数，而函数和对于是凸函数，且函数对于是凸函数，则该函数是多凸的。然而，众所周知，该条件过于严格，出于几个原因应该用比多凸性更弱的条件，即准凸性[41]来替代。然而，在局部增长条件下，即作为的情况下，还没有已知的准凸性有用描述。

本文的结构如下：第2节发展了具有标量内部变量的纤维增强超弹性材料的本构框架，包括材料模型及其准静态特性。第3节描述了在Abaqus/Standard中通过UMAT公式化和相关的速率形式本构方程的有限元实现。第4节展示了均匀测试和非均匀膨胀类问题的数值结果，说明了耗散、速率效应和不稳定性类响应。第5节总结了主要结论，并概述了未来工作的方向，特别是关于永久变形和动态应用。

2. 纤维增强材料模型
2.1. 一般框架
纤维增强材料模型通过对Helmholtz能量的加法分解来制定[42,43]，形式为(11) 非弹性行为通过一个标量量来表示，类似于在连续损伤力学（CDM）框架内制定的材料模型的本构关系[44]。标准CDM模型与这里提出的模型之间的主要区别在于能够描述材料损伤过程的可逆性[45]。我们提出了一个Helmholtz自由能函数，类似于在[46]中采用的函数，其中使用各向同性Murnaghan本构模型研究了波传播问题。如果我们有(12)，那么根据(9)，我们得到(13) 和 (14) 假设变量g的演化方程满足不等式(13) (15) 实际上，如果，则 (16) 这清楚地表明两种情况都可能发生，这意味着标量内部变量使得损伤类过程的可逆描述成为可能。在或的情况下，模型不预测耗散，即它描述了一个热力学上可逆的过程。从热力学要求来看，还得出(17) 我们假设，并且进一步假设在极限情况下恢复超弹性关系。因此。鉴于初始条件为零耗散的自然状态假设和，(15) 导致(18) 对于这类势能，可以简洁地陈述存在弱解和物理上真实的波传播的结构要求。假设纯弹性贡献W是多凸的，因此在感兴趣的域中是一阶凸的（强椭圆的）。降解函数必须是非负的、凸的、非增强的()，并且满足。势能被假设为严格凸的，随着g接近允许区间的下限。在这些条件下，对于每个固定的g，Helmholtz函数保持多凸，而对于每个固定的，g是凸的，同时通过的屏障属性保持了刚性。此外，由于瞬时弹性刚度与正因子成比例，冻结损伤响应继承了W的一阶凸性。因此，在每个允许的损伤水平下都保持了强椭圆性。对于的速率型演化律提供了粘性正则化，防止了在速率无关软化中经常观察到的病态局部化。

与能量各向异性部分相关的函数假设为形式(19)，其中。如文章前一节所解释的，参数张量由与纤维族的不同方向相关的单位向量定义。对于函数（11），它导出了一个横向各向同性材料的简化模型[47,48,49]。所提出模型类的本构关系可以解释为由粘弹性基质和弹性纤维组成的复合材料。应当注意的是，以类似的方式，我们也可以建立描述具有弹性基质和粘弹性纤维模型的方程。2.2. 材料模型前面方程中出现的函数与超弹性有关。在这里，我们采用了一种各向同性的函数形式，即广义Ciarlet模型[50,51]：(20)其中是一个凸函数。根据自然状态的假设，可以得出(21)关于基本Ciarlet模型，即当(52)时，函数（20）预测了对给定剪切变形的不同定性响应。我们假设函数的形式为(22)对于参数值和，储存能量函数定义了一种基于多凸性的公式，确保在生长条件下超弹性的解的存在。如果我们考虑中等程度的变形，函数（20）可以相对于应变张量进行二次近似[52]。该模型将参数与初始体积模量和Lamé常数联系起来，使得(23)如果我们还假设和(24)，那么我们得到(25)我们还可以假设参数c的形式。因此，模型参数可以用弹性常数来表示。如果我们还假设，那么储存能量函数可以通过进行调整，我们可以使用一组无量纲参数，因为(26)参数值必须满足约束(27)从多凸性的角度来看[35]，我们考虑函数如下：(28)其中可以解释为纤维家族的初始杨氏模量。根据上述要求，与材料模型的耗散部分相关的函数可以采用[45]的形式(29)其中参数是正常数标量，具有初始配置中单位体积的能量维度。通过与Ciarlet模型参数类似的归一化方法，对于此处考虑的中等程度变形计算，引入无量纲参数和是方便的，通过关系和来定义。在本文的其余部分，“FRD模型”一词将表示由Helmholtz函数（11）以及纤维描述（28）定义的本构关系。在该模型中，Helmholtz函数的耗散部分被理解为（12）与广义Ciarlet模型（20）和（29）结合。2.3. 具有恒定变形率的准静态问题为了说明所考虑模型在各向同性部分的一些基本特性，本小节检查了一个具有恒定变形率的准静态变形问题。所考虑的变形类型是三轴压缩/拉伸，即。在这种情况下，Cauchy应力张量也是各向同性的，表示为。拉伸比根据图1所示的图表定义。图1. 随时间尺度变化的预设拉伸：(a)两个间隔和(b)四个内部间隔。问题在假设参数值、、和的情况下得到解决。与模型耗散部分相关的参数取为、、。这些值的选择是为了获得不同于文献中报告的损伤力学[53,54]中的解决方案。由于变量g的演化方程是非线性的，解决方案是使用Wolfram Mathematica环境中的NDSolve程序确定的[55]。图2中的图表显示了在两个间隔上定义的拉伸情况下，分别对于的重新缩放应力值，即单个加载/卸载循环。结果的明显差异不仅来自于随着的增加，模拟了更硬的材料响应，还因为较大的弹性能量显著影响内部变量g的解决方案。图2. 在具有恒定变形率的两个间隔的问题中，的(a)和(b)的缩放应力。图3和图4展示了图1b中所示的预设拉伸的内部变量g和重新缩放应力作为的函数所获得的值。从变量g的演化方程的形式来看，解决方案直接取决于储存能量函数的值，这可以在g的图表行为中观察到。对于，弹性能量的值显著大于对于，这导致在时间间隔和之间获得的应力大大减小。这可以解释为模拟材料的更明显的损伤。图3. 在具有分段恒定速率的问题中，内部变量(a)g和(b)。图4. 在具有分段恒定速率的问题中，内部变量(a)g和(b)。3. 材料模型实现3.1. Abaqus/Standard用户子程序Abaqus提供了多种实现用户定义材料模型本构关系的方法[56]。主要工具是所谓的用户子程序，它们是用Fortran或C/C++编写的，通过明确定义的接口与Abaqus核心模块交互。这些子程序可以修改标准程序或引入自定义材料模型等。从超弹性的角度来看，UHYPER和UANISOHYPER_INV子程序特别重要，因为它们允许通过指定储存能量密度函数及其相对于适当变形不变量的导数来定义各向同性和各向异性材料的本构关系[57]。基于这些量，程序确定Cauchy应力张量的分量和切线刚度张量的分量，以满足客观性要求。(30)在Abaqus/Standard中，使用了Zaremba–Jaumann速率[26]。对于Kirchhoff应力张量，有(31)对于包含耗散效应的本构关系，使用UMAT子程序接口。它需要定义Cauchy应力张量的分量和由本构关系产生的切线刚度张量的分量。如果将切线刚度张量理解为(30)中的，则实现与超弹性方程一致，这种方法被称为总公式[26]。在这种情况下，可以使用具有混合公式的有限元，尽管需要以类似于UHYPER程序的方式定义适当的三分量导数。3.2. 本构方程的速率形式FRD模型的本构关系的增量形式，包括Zaremba–Jaumann目标速率，可以表示为(32)其中(33)与广义Ciarlet模型和不变量的多项式相关联的切线刚度张量已在[43]中确定。因此，评估(32)中的张量还需要确定表达式。从内部变量的演化方程采用的形式来看，可以直接将视为变量的函数。因此，随之而来的是(34)连同(35)方程（35）的实现是通过在积分步骤上应用前向欧拉方法来执行的，使用(36)如果，则导数是常数，即。在这种情况下，得到一个直接的关系(37)完整的本构模型，包括与标量内部变量相关的耗散贡献，通过UMAT子程序在Abaqus/Standard中实现，并在下一节中介绍的有限元研究中使用。4. 结果4.1. 数值验证4.1.1. 超弹性首先通过单轴应变和简单剪切试验在单个C3D8R有限元上验证了通过UMAT子程序实现FRD模型本构关系的正确性，没有考虑非弹性效应。所有边界条件都指定为位移控制。超弹性模型的参数采用为。计算出的应力值与从解析表达式获得的值在规定的数值精度范围内一致。图5和图6展示了Cauchy应力张量重新缩放分量的图表。图5. 单轴应变问题：(a)纤维方向，(b)纤维方向。图6. 简单剪切问题：(a)纤维方向，(b)纤维方向。4.1.2. 具有耗散效应的纤维增强模型验证FRD模型实现的第二步是包含粘性效应。为此，进行了松弛和蠕变试验，前者采用参数值，后者采用。松弛问题涉及具有位移边界条件的一轴应变状态。变形梯度张量采用的形式，纤维方向为。在蠕变问题中，考虑平面应力，即应力状态由张量定义，变形状态由定义。假设通过在中未变形的有限元素表面上规定牵引向量来控制“11”应力分量。图7展示了松弛问题中内部变量和Cauchy应力张量重新缩放的非零分量的图表。值得注意的是，在考虑的时间尺度范围内，分量的值没有变化，因为各向异性部分中的本构关系不提供松弛。图7. (a) 内部变量和(b)松弛问题中Cauchy应力张量的非零分量作为时间尺度的函数—FRD模型。在蠕变问题中，纤维方向被假设与向量对齐。图8和图9显示了作为时间函数的内部变量和拉伸的结果值。使用Abaqus获得的结果与通过Mathematica中的NDSolve程序计算的结果一致。图8. (a) 内部变量和(b)在蠕变问题中作为时间尺度的函数的缩放应力—FRD模型。图9. 在蠕变问题中，作为时间函数的拉伸(a)、(b)。4.2. 非均匀变形考虑了一个半径为R、厚度为的圆形几何体的膨胀问题，该问题受到边界条件的约束。分析是针对参数值的FRD模型进行的。假设纤维家族方向与x方向对齐。由于实现的UMAT程序专用于三维应力状态，因此使用了C3D8R类型的矩形、八节点、简化积分的砖元素。这种元素类型不与壳理论公式相关[58]。通过在圆形的上表面规定压力来施加加载和卸载，如图10所示。离散域由44个有限元素组成。图10. (a) 有限元素网格；(b) 膨胀问题中的重新缩放压力。图11显示了作为区域中点位移函数的归一化压力的图表。值得注意的是，在最终卸载阶段，位移值出现了“跳变”到正值，其中也展示了下一个加载周期。这可以解释为一种不稳定性；见[59]。图11显示了作为拉伸函数的重新缩放应力值。这些值是从其中一个中心元素提取的，清楚地表明了能量耗散。在几个选定的元素中评估的内部变量在0.15到0.30之间达到最大值；见图12。图11. (a) 作为中点位移函数的归一化压力；(b) 重新缩放的应力。图12. 在选定的有限元素中作为时间函数的内部变量。5. 结论提出了一种超弹性模型本构关系的推广，其中包括纤维增强和粘性效应，涉及引入一个以标量量形式表示的内部变量，类似于在连续损伤力学（CDM）框架内制定的材料模型的本构关系。这些模型属于纤维增强材料模型[38]的类别，并在简单混合理论[42]的框架内额外制定。亥姆霍兹自由能函数各向同性分量的形式借鉴自[45]，随后被扩展以涵盖广义Ciarlet模型和纤维增强效应，从而形成一个具有耗散效应的各向同性材料模型[28]。CDM模型与本文所提出的模型类别之间的主要区别在于后者能够描述材料损伤过程的可逆性。正如本文讨论的边界值问题所示，这类模型不仅可以模拟材料损伤现象，还能捕捉应变率敏感性以及蠕变和应力松弛现象。这些模型的主要局限性在于它们无法表征永久变形。这一问题，以及基于所提出的材料模型对初始边界值问题中动态效应的研究，将在未来的工作中得到解决。