INVESTIGATION OF VISCO-ELASTIC MECHANICAL BEHAVIORS OF UV-CURING RESIN BASED ON FINITE VOLUME METHOD
-
摘要: 光固化聚合物体系成分复杂, 随时间变化的黏弹性力学行为具有不确定性. 因此, 以代表性光固化树脂中的聚氨酯丙烯酸树脂(polyurethane acrylate, PUA)为例, 基于二阶位移展开有限体积法及分数阶黏弹性本构关系发展光固化树脂的力学模型. 首先, 根据边界条件确定胞元面的连续性条件和数值模型的总刚矩阵, 显式计算所有胞元的未知面位移. 为探究光固化树脂在不同应变率影响下的黏弹性力学行为, 引入分数阶有限变形Kelvin-Voigt黏弹性力学模型, 建立应力-应变之间的本构关系. 然后, 通过应变率分别为10−4, 10−3, 10−2以及10 s−1的单向拉伸试验完成本构模型中弹性项以及黏性项的参数识别. 最后, 开展不同应变率下的单向拉伸试验及数字图像散斑(digital image correlation, DIC)分析, 评估数值模型的计算误差. 结果表明, 文章所建立的数值模型能有效预测不同应变率下光固化树脂的力学行为, 并且平均预测误差仅为1.98%.
-
关键词:
- 聚氨酯丙烯酸树脂 /
- 有限体积法 /
- 分数阶Kelvin-Voigt黏弹性力学模型 /
- 光固化树脂 /
- 参数识别
Abstract: The complexity of the ultraviolet-curing (UV-curing) polymer system, comprised of various resins, contributes to the uncertainty in its viscoelastic mechanical behaviors as it evolves over time. By using polyurethane acrylate (PUA) as a representative material, the fractional order viscoelastic constitutive is incorporated into the second-order displacement expanded finite volume model to explore the mechanical behaviors of UV-curing resin in this paper. Firstly, the continuity conditions of cell surface and the global stiffness matrix of the numerical model are determined in accordance with loading conditions, enabling the explicit calculation of unknown surface displacements for all cells. To investigate the viscoelastic behavior of UV-curing resin with respect to varying strain rates, a fractional-order finite deformation Kelvin-Voigt viscoelastic model is introduced to establish the stress-strain constitutive relationship. Secondly, parameters of the elastic and viscous components in the constitutive model are determined through uniaxial tensile tests at strain rates of 10−4, 10−3, 10−2, and 10 s−1. Finally, digital image correlation (DIC) analysis is performed in conjunction with uniaxial tensile tests at various strain rates to compare the experimental results and evaluate the computational accuracy of the numerical model. It is indicated that the established numerical model can effectively predict the visco-elastic mechanical behaviors of light-cured resin under different strain rates, with an average prediction error of 1.98%. -
引 言
树脂的光固化技术是一种重要的成型工艺, 其基本原理是以液态光敏树脂为原料, 利用材料逐层叠加的原理, 在紫外光束的精确控制下, 触发液体光敏树脂中的预聚物发生聚合和交联等反应, 从而在极短时间内转化为稳定的网状高分子聚合物[1]. 光固化成型工艺具有许多优点, 如成型速度快、加工精度高、材料利用率高和能制作具有复杂结构的几何模型等. 然而, 在长期使用时会引起由树脂的黏弹性特性而导致的蠕变变形或应力松弛, 其力学响应将影响打印件的刚度、强度以及尺寸稳定性[2-3]. 具体而言, 光固化树脂的弹性响应体现在它能够在外界载荷作用下发生变形, 并在一定范围内恢复原状, 且此过程不伴随能量损耗[4]. 同时, 光固化聚合物的非线性黏性行为则通过应变能函数得以描述, 其在加载过程中展现出了明显的率相关性[5]. 因此, 准确评估树脂的黏弹性响应对于预测光固化结构件的非线性力学行为以及破坏强度至关重要.
Du等[6]借助Prony 级数理论研究了丙烯酸材料随时间变化的机械响应以及黏弹性损伤演变机制. 然而, 建立的本构模型需要多达几十项不等的Prony级数刻画黏弹性行为[7], 为参数的识别带来较大的困难. Alizadeh等[8]利用广义Maxwell模型拟合光固化树脂的宏观力学行为, 通过试验分析光固化聚合物中丙烯酸共聚物和聚氨酯的比例关系对宏观率相关力学响应的影响. 然而由于试验较为复杂, 影响因素较多, 参数识别误差较大. 与整数阶理论相比, 分数阶理论在精确描述黏弹性力学行为时需要的本构参数较少, 且参数的物理意义明确[9-10]. 分数阶模型具有记忆效应[11], 综合考量了历史加载的影响[12], 适用于聚合物的黏弹性力学建模[13-14]. Drozdov[15]针对有限应变黏弹性介质构建了Kelvin-Voigt, Maxwell以及Maxwell-Weichert 等分数阶微分本构模型, 通过杆件的单轴拉伸试验验证了本构模型的预测精度; 与其他模型相比, 该模型可以显著减少可调参数的数量. Sumelka等[16]采用分数阶形式表述黏弹性材料的损伤演化规律, 基于历史变量的损伤信息可以快速评估加载过程中的材料强度. 尹耀得等[17]基于分数阶有限变形Kelvin-Voigt 流变学模型建立弹性体的三维张量本构, 推导了单向拉伸下的本构关系. 因此, 分数阶黏弹性力学模型在聚合物的非线性力学行为预测上具有较强的优势.
材料的本构模型建立了应力-应变之间的力学响应关系, 高效的数值模型可以利用分数阶黏弹性力学模型的特点, 精确预测整体结构变形和局部力学响应. 有限体积法(finite volume method, FVM)可以快速预测连续介质体内的力学场物理量分布, 虽然在固体力学中使用FVM并不常见, 但已有研究证明, 该方法在精度上可与流行的有限元方法相媲美[18]. Demirdzic等[19]首次将FVM用以解决固体力学问题, 提出将胞元中心作为积分点的隐式求解方法. FVM强调在胞元内部变量的强平衡形式, 确保了局部守恒的严格性, 进而在全局范围内实现了守恒[20-21]. 相比之下, 有限元方法使用弱形式的积分平衡以获得单元的偏微分方程近似解, 实现了数值模型平均意义上的局部守恒, 而非每个单元的直接守恒[22]. 因此, FVM可以采用较少的离散胞元, 实现高精度局部场的物理量预测[23]. 此外, Cardiff等[24]通过二维线弹性数值模型验证了FVM的计算效率优势, 证明了FVM在解决固体力学的大模型计算时极具潜力[25]. Scolaro等[26]提出了一种直接包含摩擦力和非正交边界修正向量的法向接触应力隐式离散化方案, 通过FVM解决接触仿真的问题. Soar等[27]将元胞自动机方法与二维FVM相结合, 模拟了结晶过程的微观力学行为. Ivanova等[22]在有限容积模型中增加了描述黏弹性材料动态和准静态行为的源项, 进一步拓展了FVM的应用范围.
为评估光固化树脂的非线性力学行为, 本文以光固化树脂之一的聚氨酯丙烯酸树脂(polyurethane acrylate, PUA)为例, 采用分数阶黏弹性本构模型建立FVM数值仿真框架. 首先, 建立三维数值模型并引入分数阶黏弹性本构模型; 然后, 通过单向拉伸试验获得应变率分别为10−4, 10−3, 10−2和10 s−1情况下的真实应力-应变曲线, 并基于试验结果识别PUA树脂弹性及黏性阶段的本构模型参数; 最后, 采用DIC技术捕捉PUA标准样件在单向拉伸载荷下的位移云图, 评估理论模型的准确性.
1. 理论模型
1.1 基于二阶位移展开的FVM数值计算框架
本文采用胞心法, 即未知的物理量位于胞元的表面, 并将胞元的中心作为局部变量的平衡中心. 在Castrillo等[28]关于一阶位移展开工作的基础上, 采用高阶的位移展开提高大变形问题中数值计算的精度. 因此, 本研究对三维胞元的面位移采用二阶Legendre多项式展开
$$\begin{split} & u_i^{(n)}(\zeta ,\eta ,\xi ) = u_{i(000)}^{(n)} + \zeta u_{i(100)}^{(n)} + \eta u_{i(010)}^{(n)} + \\ &\qquad \xi u_{i(001)}^{(n)} + \frac{1}{2}(3{\zeta ^2} - 1)u_{i(200)}^{(n)} + \frac{1}{2}(3{\eta ^2} - 1)u_{i(020)}^{(n)} + \\ &\qquad \frac{1}{2}(3{\xi ^2} - 1)u_{i(002)}^{(n)}{\text{ }}\quad (i = 1,2,3) \end{split} $$ (1) 式中, ζ, η和ξ分别为参数化坐标系的坐标, 取值范围[−1,1]. 上标n为单元的编号. 下标i为对应的坐标系方向. 对六面体胞元采用一阶线性插值, 建立胞元顶点在笛卡尔坐标系与参数化坐标系之间的映射关系
$$\begin{split} & {\varepsilon _{ij}}\left(y\right) = \frac{1}{2}\left(\frac{{\partial {u_i}}}{{\partial {y_j}}} + \frac{{\partial {u_j}}}{{\partial {y_i}}}\right) = \frac{1}{2}\left({J_{i1}}\frac{{\partial {u_i}}}{{\partial \zeta }} + {J_{i2}}\frac{{\partial {u_i}}}{{\partial \eta }} + \right. \\ &\qquad \left.{J_{i3}}\frac{{\partial {u_i}}}{{\partial \xi }}\right) + \frac{1}{2}\left({J_{j1}}\frac{{\partial {u_j}}}{{\partial \zeta }} + {J_{j2}}\frac{{\partial {u_j}}}{{\partial \eta }} + {J_{j3}}\frac{{\partial {u_j}}}{{\partial \xi }}\right) \end{split} $$ (2) 式中, 上标p表示对应的顶点标号. 型函数N的展开式如下
$$ \left. \begin{aligned} & {{N_{1,2}}(\zeta ,\eta ,\xi ) = \frac{1}{8}(1 - \zeta )(1 \mp \eta )(1 - \xi )} \\ & {{N_{3,4}}(\zeta ,\eta ,\xi ) = \frac{1}{8}(1 - \zeta )(1 \pm \eta )(1 + \xi )} \\ & {{N_{5,6}}(\zeta ,\eta ,\xi ) = \frac{1}{8}(1 + \zeta )(1 \mp \eta )(1 - \xi )} \\ & {{N_{7,8}}(\zeta ,\eta ,\xi ) = \frac{1}{8}(1 + \zeta )(1 \pm \eta )(1 + \xi ) } \end{aligned} \right\} $$ (3) 根据式(1), 在笛卡尔坐标系下对n-胞元的左、右、后、前、上和下面位移分别做积分处理
$$\left.\begin{split} & \bar u_i^{({\mathrm{left}},{\mathrm{right}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {{u_i}( \mp 1,\eta ,\xi )} } {\mathrm{d}}\eta {\mathrm{d}}\xi = \\ &\qquad {u_{i(000)}} \mp {u_{i(100)}} + {u_{i(200)}} \\ & \bar u_i^{({\mathrm{rear}},{\mathrm{front}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {{u_i}(\zeta , \mp 1,\xi )} } {\mathrm{d}}\zeta {\mathrm{d}}\xi = \\ &\qquad {u_{i(000)}} \mp {u_{i(010)}} + {u_{i(020)}} \\ & \bar u_i^{({\mathrm{lower}},{\mathrm{upper}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {{u_i}(\zeta ,\eta , \mp 1)} } {\mathrm{d}}\zeta {\mathrm{d}}\eta = \\ &\qquad {u_{i(000)}} \mp {u_{i(001)}} + {u_{i(002)}} \end{split}\right\} $$ (4) 由此, n-胞元面位移的一阶项和二阶项可以表示为
$$ \begin{split} &{\left[ {\begin{array}{*{20}{c}} {{u_{i(100)}}} \\ {{u_{i(010)}}} \\ {{u_{i(001)}}} \\ {{u_{i(200)}}} \\ {{u_{i(020)}}} \\ {{u_{i(002)}}} \end{array}} \right]^{(n)}} = \\ &\frac{1}{2}\left[ {\begin{array}{*{20}{c}} 1&{ - 1}&0&0&0&0 \\ 0&0&1&{ - 1}&0&0 \\ 0&0&0&0&1&{ - 1} \\ 1&1&0&0&0&0 \\ 0&0&1&1&0&0 \\ 0&0&0&0&1&1 \end{array}} \right]{\left[ {\begin{array}{*{20}{c}} {\bar u_i^{{\mathrm{left}}} - {u_{i(000)}}} \\ {\bar u_i^{{\mathrm{right}}} - {u_{i(000)}}} \\ {\bar u_i^{{\mathrm{rear}}} - {u_{i(000)}}} \\ {\bar u_i^{{\mathrm{front}}} - {u_{i(000)}}} \\ {\bar u_i^{{\mathrm{lower}}} - {u_{i(000)}}} \\ {\bar u_i^{{\mathrm{upper}}} - {u_{i(000)}}} \end{array}} \right]^{(n)}}\end{split} $$ (5) 为了获取n-胞元的面应变, 需建立面位移在笛卡尔坐标系与参数化坐标系下的映射关系, 构建Jacobi矩阵
$$ \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial \zeta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \eta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \xi }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {y_1}}}{{\partial \zeta }}}&{\dfrac{{\partial {y_2}}}{{\partial \zeta }}}&{\dfrac{{\partial {y_3}}}{{\partial \zeta }}} \\ {\dfrac{{\partial {y_1}}}{{\partial \eta }}}&{\dfrac{{\partial {y_2}}}{{\partial \eta }}}&{\dfrac{{\partial {y_3}}}{{\partial \eta }}} \\ {\dfrac{{\partial {y_3}}}{{\partial \xi }}}&{\dfrac{{\partial {y_3}}}{{\partial \xi }}}&{\dfrac{{\partial {y_3}}}{{\partial \xi }}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial {y_1}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_2}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_3}}}} \end{array}} \right] = {{{\boldsymbol{J}}}}\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial {y_1}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_2}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_3}}}} \end{array}} \right] $$ (6) 由式(6)得
$$ \left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial {y_1}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_2}}}} \\ {\dfrac{{\partial {u_i}}}{{\partial {y_3}}}} \end{array}} \right] = {{\boldsymbol{J}}^{ - 1}}\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial \zeta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \eta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \xi }}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{J_{11}}}&{{J_{12}}}&{{J_{13}}} \\ {{J_{21}}}&{{J_{22}}}&{{J_{23}}} \\ {{J_{31}}}&{{J_{32}}}&{{J_{33}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {\dfrac{{\partial {u_i}}}{{\partial \zeta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \eta }}} \\ {\dfrac{{\partial {u_i}}}{{\partial \xi }}} \end{array}} \right] $$ (7) 根据式(7)可得n-胞元在笛卡尔坐标系下的面应变为
$$ \begin{split} & {\varepsilon _{ij}}\left(y\right) = \frac{1}{2}\left(\frac{{\partial {u_i}}}{{\partial {y_j}}} + \frac{{\partial {u_j}}}{{\partial {y_i}}}\right) = \frac{1}{2}\left({J_{i1}}\frac{{\partial {u_i}}}{{\partial \zeta }} + {J_{i2}}\frac{{\partial {u_i}}}{{\partial \eta }} + \right. \\ &\qquad \left.{J_{i3}}\frac{{\partial {u_i}}}{{\partial \xi }}\right) + \frac{1}{2}\left({J_{j1}}\frac{{\partial {u_j}}}{{\partial \zeta }} + {J_{j2}}\frac{{\partial {u_j}}}{{\partial \eta }} + {J_{j3}}\frac{{\partial {u_j}}}{{\partial \xi }}\right) \end{split} $$ (8) 根据黏弹性力学本构模型, 建立n-胞元的面应力与面应变以及面应变率之间函数关系
$$ {{\boldsymbol{\sigma}} } = {{{\boldsymbol{\sigma}} }_e}\left( {{{\boldsymbol{\varepsilon}} }} \right) + {{{\boldsymbol{\sigma}} }_v}\left( {{{\dot {\boldsymbol{\varepsilon}} }}} \right) $$ (9) 式中, $ {{\boldsymbol{\sigma}} _e} $和$ {{\boldsymbol{\sigma}} _v} $分别应力的弹性项和黏性项. 根据胞元的面应力$ {\boldsymbol{\sigma}} $以及对应面的方向向量m, 可计算对应胞元面的牵引力
$$ {\boldsymbol{t}} = {\boldsymbol{m}}{\boldsymbol{\sigma}} $$ (10) 式中, 方向向量m为对应面的外法线向量与笛卡尔坐标系夹角的余弦值. 同样, 对n-胞元的面牵引力分别进行积分, 可得平均面牵引力为
$$ \left.\begin{split} & \bar t_i^{({\mathrm{left}},{\mathrm{right}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {t_i^{({\mathrm{left}},{\mathrm{right}})}( \mp 1,\eta ,\xi )} } {\mathrm{d}}\eta {\mathrm{d}}\xi \\ & \bar t_i^{({\mathrm{rear}},{\mathrm{front}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {t_i^{({\mathrm{rear}},{\mathrm{front}})}(\zeta , \mp 1,\xi )} } {\mathrm{d}}\zeta {\mathrm{d}}\xi \\ & \bar t_i^{({\mathrm{lower}},{\mathrm{upper}})} = \frac{1}{4}\int\nolimits_{ - 1}^1 {\int\nolimits_{ - 1}^1 {t_i^{({\mathrm{lower}},{\mathrm{upper}})}(\zeta ,\eta , \mp 1)} } {\mathrm{d}}\zeta {\mathrm{d}}\eta \end{split}\right\} $$ (11) 根据高斯散度定理[29], n-胞元的面应力之和为零, 即
$$\begin{split} & \left[ {\bar t_i^{{\mathrm{left}}} + \bar t_i^{{\mathrm{right}}} + \bar t_i^{{\mathrm{upper}}} + \bar t_i^{{\mathrm{lower}}} + \bar t_i^{{\mathrm{roar}}}} \right. + \\ &\qquad {\left. { \bar t_i^{{\mathrm{front}}}} \right]^{{\text{(}}n{\text{)}}}} = \int_\varOmega {\nabla \cdot {\boldsymbol{\sigma}} } {\mathrm{d}}\varOmega = \oint\nolimits_\varGamma {{\boldsymbol{m}} \cdot {\boldsymbol{\sigma}} } {\mathrm{d}}\varGamma = 0\end{split} $$ (12) 将式(5)、式(8)和式(10)代入上式, 可以求得n-胞元的面位移零阶分量的表达式
$$ {\left[ {\begin{array}{*{20}{c}} {{u_{1(000)}}} \\ {{u_{2(000)}}} \\ {{u_{3(000)}}} \end{array}} \right]^{(n)}} = {{\boldsymbol{\varPhi}} ^{ - 1}}{\boldsymbol{\varTheta}} {\left[ {\begin{array}{*{20}{c}} {\bar u_1^{{\mathrm{left}}} + \bar u_1^{{\mathrm{right}}}} \\ {\bar u_2^{{\mathrm{left}}} + \bar u_2^{{\mathrm{right}}}} \\ {\bar u_3^{{\mathrm{left}}} + \bar u_3^{{\mathrm{right}}}} \\ {\bar u_1^{{\mathrm{rear}}} + \bar u_1^{{\mathrm{front}}}} \\ {\bar u_2^{{\mathrm{rear}}} + \bar u_2^{{\mathrm{front}}}} \\ {\bar u_3^{{\mathrm{rear}}} + \bar u_3^{{\mathrm{front}}}} \\ {\bar u_1^{{\mathrm{lower}}} + \bar u_1^{{\mathrm{upper}}}} \\ {\bar u_2^{{\mathrm{lower}}} + \bar u_2^{{\mathrm{upper}}}} \\ {\bar u_3^{{\mathrm{lower}}} + \bar u_3^{{\mathrm{upper}}}} \end{array}} \right]^{(n)}} $$ (13) 式中, 矩阵$ {\boldsymbol{\varPhi}} $和$ {\boldsymbol{\varTheta}} $的展开式参见参考文献[29]. 将式(13)代入式(12)中可以得到n-胞元面位移的零阶、一阶以及二阶分量的表达式. 由此, 建立n-胞元面平均位移与面平均牵引力之间的局部刚度矩阵关系
$$ {\bar {\boldsymbol{t}}^{(n)}} = {\boldsymbol{K}}_{{\mathrm{cell}}}^{(n)}{\bar {\boldsymbol{u}}^{(n)}} $$ (14) 根据试验样件的几何结构采用映射网格的方式建立数值模型, 数值模型分别沿笛卡尔坐标系的y1, y2与y3方向分别被等分为Nα, Nβ以及Nγ个胞元. 如图1所示, 根据真实的加载情况对蓝色区域内胞元的面位移进行完全约束, 对黄色区域内的胞元沿y1方向施加拉伸载荷. 根据胞元之间的承载情况确定对应位移与牵引力的连续性条件, 建立数值模型的总刚矩阵. 已知数值模型的胞元总数为Nα × Nβ × Nγ个, 则未知胞元面位移为3(Nα + 1) × Nβ × Nγ + 3Nα × (Nβ + 1) × Nγ + 3Nα × Nβ × (Nγ + 1) 个. 如图2所示, 则可针对数值模型内部的胞元面建立面位移连续性条件, 即
![]()
图 2 相邻胞元间面位移的连续性条件Figure 2. Continuous conditions of surface displacement between adjacent cells$$ \left. \begin{split} & {\bar u_i^{{\mathrm{left}}(\alpha \beta \gamma )} = \bar u_i^{{\mathrm{right}}(\alpha - 1,\beta ,\gamma )}{\text{ = }}u_i^{1(\alpha \beta \gamma )}} \\ & {\bar u_i^{{\mathrm{rear}}(\alpha \beta \gamma )} = \bar u_i^{{\mathrm{front}}(\alpha ,\beta - 1,\gamma )}{\text{ = }}u_i^{2(\alpha \beta \gamma )}} \\ & {\bar u_i^{{\mathrm{lower}}(\alpha \beta \gamma )} = \bar u_i^{{\mathrm{upper}}(\alpha ,\beta ,\gamma - 1)}{\text{ = }}u_i^{3(\alpha \beta \gamma )}} \\ &\qquad {(i = 1,2,3)} \end{split} \right\} $$ (15) 对于数值模型表面的面位移, 由于样件受到沿y1方向的单向拉伸载荷, 因此数值模型表面胞元沿y2与y3方向的面位移可表示为
$$ \left. \begin{split} & u_1^{2(\alpha ,1,\gamma )} = u_1^{2(\alpha ,{N_\beta } + 1,\gamma )}\\ & u_1^{3(\alpha ,\beta ,1)} = u_1^{3(\alpha ,\beta ,{N_\gamma } + 1)} \\ & u_2^{2(\alpha ,1,\gamma )} + u_2^{2(\alpha ,{N_\beta } + 1,\gamma )} = 0\\ &u_2^{3(\alpha ,\beta ,1)} = u_2^{3(\alpha ,\beta ,{N_\gamma } + 1)} \\ & u_3^{2(\alpha ,1,\gamma )} = u_3^{2(\alpha ,{N_\beta } + 1,\gamma )}\\ &u_3^{3(\alpha ,\beta ,1)} + u_3^{3(\alpha ,\beta ,{N_\gamma } + 1)} = 0 \end{split} \right\} $$ (16) 由于数值模型端面垂直于y1方向的面位移被约束, 可得
$$ \left.\begin{split} & {u_1^{1(1,\beta ,\gamma )} = 0,{\text{ }}u_1^{1({N_\alpha } + 1,\beta ,\gamma )} = U} \\ & {u_2^{1(1,\beta ,\gamma )} = u_2^{1({N_\alpha } + 1,\beta ,\gamma )}} \\ & {u_3^{1(1,\beta ,\gamma )} = u_3^{1({N_\alpha } + 1,\beta ,\gamma )} } \end{split} \right\} $$ (17) 式中, U为样件端部的位移. 根据式(16)和式(17), 未知的胞元面位移的数量减少为9Nα × Nβ × Nγ. 同样, 需要确定相邻胞元面牵引力的连续性条件. 由于样件在y2与y3方向的面牵引力不受外载作用, 可建立数值模型内部相邻胞元的面牵引力以及外表面胞元的牵引力关系
$$ \left.\begin{split} & {\bar t_i^{{\mathrm{rear}}(\alpha \beta \gamma )} + \bar t_i^{{\mathrm{front}}(\alpha ,\beta - 1,\gamma )}{\text{ = }}0} \\ & {\bar t_i^{{\mathrm{rear}}(\alpha ,{N_\beta },\gamma )} + \bar t_i^{{\mathrm{front}}(\alpha ,1,\gamma )}{\text{ = }}0} \\ & {\bar t_i^{{\mathrm{lower}}(\alpha \beta \gamma )} + \bar t_i^{{\mathrm{upper}}(\alpha ,\beta ,\gamma - 1)}{\text{ = }}0} \\ & {\bar t_i^{{\mathrm{lower}}(\alpha ,\beta ,{N_\gamma })} + \bar t_i^{{\mathrm{upper}}(\alpha ,\beta ,1)}{\text{ = }}0} \\ &\qquad {(i = 1,2,3)} \end{split}\right\} $$ (18) 通过图3所示的胞元面牵引力的受力分析可知, 对于加载方向的相邻胞元间面牵引力的合力应与外载的作用力方向相同, 可得
$$ \left.\begin{split} & {\bar t_i^{{\mathrm{left}}(\alpha \beta \gamma )} + \bar t_i^{{\mathrm{right}}(\alpha - 1,\beta ,\gamma )}{\text{ = }}T_i^{1(\alpha - 1,\beta ,\gamma )}{\text{ = }}2T} \\ & {\bar t_i^{{\mathrm{left}}({N_\alpha },\beta ,\gamma )} + \bar t_i^{{\mathrm{right}}(1,\beta ,\gamma )}{\text{ = }}T_i^{1(1,\beta ,\gamma )}{\text{ = }}2T} \\ &\qquad {(i = 1,2,3)} \end{split}\right\} $$ (19) 式中, $ {\text{2}}T_i^{1(\alpha \beta \gamma )} $为相邻胞元之间面牵引力的合力, T为对样件施加的外载荷.
联立式(15) ~ 式(19), 基于总刚矩阵$ {{\boldsymbol{K}}_{{\text{global}}}} $建立数值模型的未知胞元面位移$ {\boldsymbol{U}} $与外载荷$ {\boldsymbol{T}} $之间的关系为
$$ {{\boldsymbol{K}}_{{\text{global}}}}{\boldsymbol{U}} = {\boldsymbol{T}} $$ (20) ![]()
图 3 沿加载方向相邻胞元之间面牵引力的连续性条件Figure 3. Continuous conditions of surface traction force between adjacent cells along the loading direction1.2 针对光固化树脂的黏弹性本构模型
本构模型中采用并联的黏壶和弹簧描述介质变形过程中的黏弹性力学行为[30], 如图3所示. 变形历程中$\tau $时刻的变形梯度张量Fτ定义为
$$ {{\boldsymbol{F}}_\tau } = {\boldsymbol{I}} + \frac{{\partial {\boldsymbol{u}}}}{{\partial {{\boldsymbol{x}}_\tau }}} = {\boldsymbol{I}} + {\boldsymbol{\varepsilon}} (\tau ) = {\boldsymbol{F}}({{t}}){{\boldsymbol{F}}^{ - 1}}(\tau ) $$ (21) 式中, ${{{\boldsymbol{F}}}}(t)$为现时构型的变形梯度张量, ${\boldsymbol{I}}$为单位矩阵, ${\boldsymbol{\varepsilon}} $为$\tau $时刻的应变. 变形速度张量规定为
$$ {\boldsymbol{D}} = {\mathrm{sym}}({\boldsymbol{L}}) = \frac{1}{2}({\boldsymbol{L}} + {{\boldsymbol{L}}^{\text{T}}}),\quad {\boldsymbol{L}} = \dot {\boldsymbol{F}}{{\boldsymbol{F}}^{ - 1}} $$ (22) 因此, 变形速度张量${\boldsymbol{D}}$和Cauchy-Green张量B分别为
$$ {\boldsymbol{D}} = \frac{1}{2}\left[ {\dot {\boldsymbol{F}}{{\boldsymbol{F}}^{ - 1}} + {{\left( {\dot {\boldsymbol{F}}{{\boldsymbol{F}}^{ - 1}}} \right)}^{\text{T}}}} \right],\quad {\boldsymbol{B}}{\text{ = }}{\boldsymbol{F}}{{\boldsymbol{F}}^{\text{T}}} $$ (23) 式(9)中n-胞元Cauchy应力中的弹性应力${\boldsymbol{\sigma}} _e^n$表示为
$$ {{\boldsymbol{\sigma}} }_e^n = {\left[ { - {\boldsymbol{pI}} + 2\left( {{\chi _1}{\boldsymbol{B}} + {\chi _2}{{\boldsymbol{B}}^2}} \right)} \right]^n} $$ (24) 式中, ${\boldsymbol{I}}$为单位张量; 系数${\chi _1}$和${\chi _2}$与应变能密度函数W有关, 本文选用双参数Mooney-Rivlin模型[31]的应变能密度函数$ W = {C_{10}}\left( {{I_1} - 3} \right) + {C_{01}}\left( {{I_2} - 3} \right) $,将其代入下式
$$ \left.\begin{split} & {{\chi _1} = \frac{{\partial W}}{{\partial {I_1}}} + {I_1}\frac{{\partial W}}{{\partial {I_2}}} = {C_{10}} + {C_{01}}{I_1}} \\ & {{\chi _2} = - \frac{{\partial W}}{{\partial {I_2}}} = - {C_{01}}} \end{split}\right\} $$ (25) 式中, I1和I2分别为Cauchy-Green张量的第一不变量和第二不变量. 分数阶阻尼器的黏性应力${{\boldsymbol{\sigma}} }_{{{{v}}}}^n$与速度梯度张量${\boldsymbol{D}}$成正比, 式(9)中n-胞元Cauchy应力中的弹性应力${\boldsymbol{\sigma}} _v^n$表示为
$$ {\boldsymbol{\sigma}} _v^n = 2\eta {\left( {{{\boldsymbol{D}}^{\left\{ \alpha \right\}}}} \right)^n} $$ (26) 式中η表示为黏性系数, ${{\boldsymbol{D}}^{\left\{ \kappa \right\}}}$为变形梯度张量的分数阶导数, 表示为
$$ {{\boldsymbol{D}}^{\left\{ \kappa \right\}}} = \int_{{l_t}}^t {\frac{{\left[ {{{\left( {t - \tau } \right)}^{ - \kappa }}{\boldsymbol{F}}_\tau ^{\mathrm{T}}{\boldsymbol{D}}(\tau ){{\boldsymbol{F}}_\tau }} \right]}}{{\Gamma \left( {1 - \alpha } \right)}}} {\mathrm{d}}\tau $$ (27) 式中$\Gamma $为Gamma函数. 记忆时间${l_t} \in \left[ {0,t} \right)$, 认为其黏性行为仅与近期的材料历史行为有关. 基于以上分析, 分数阶黏弹性本构模型需要通过不同的率相关试验识别关键弹性参数$ {C_{10}} $, $ {C_{01}} $以及黏性参数η和$\kappa $.
2. 参数识别
2.1 单向拉伸试验
本研究采用苏州铼赛智能科技有限公司的shape 1+®光固化打印机, 依据ASTM D638-14 [32]中的І型试验标准, 制备PUA拉伸标准样件. 使用计算机辅助设计软件Solidworks及Shapeware创建几何模型, 设置打印参数. 打印的聚合物试样的长度和厚度分别为165 ± 0.5 mm和3.2 ± 0.4 mm, 样件的几何尺寸如图4所示.
制备流程如图5所示, 首先通过表面投影数字光处理(digital light processing, DLP)技术自下而上的逐层连续制模. 为获取较高的打印精度, 打印层厚设置为50 μm, 像素尺寸为100 μm, 曝光时间为3 s, LCD紫外线光源波段为405 nm, 光强为3.7 mW/cm2. 打印完成后, 使用超声波清洗机ShapeWash 020S对预打印样件清洗2 min去除其表面的残余材料, 再将打印样件置于ShapeCure UV固化干燥箱中40 min, 箱内温度保持在30 °C. 固化后, 在室温下保存48 h. 按照上述流程共制备24个样件, 将所打印的PUA样件分为4组, 分别执行10−4, 10−3, 10−2以及10 s−1情况下的单向拉伸试验. 使用Zwick/Roell万能试验机对不同应变率工况下的3个样件执行拉伸试验, 记录真实应力-应变曲线, 识别材料参数.
为验证理论模型的预测精度, 对其余12个试验样件执行相同工况下的拉伸测试, 采用DIC (digital image correlation)技术对PUA试件的局部变形进行观测. 执行验证试验的样件需要在待观测区域内喷涂散斑, 并使用楔形夹具将其固定在试验机上. 同时采用强光对试件表面进行照射, 光源上下两侧的相机实时捕捉图像并将其传递到数据采集软件GOM Correlate中. 通过TESTXPER用户界面软件与控制计算机交互, 实现对Zwick/Roell试验机的状态控制及加载速率的设置. 万能试验机中内附的力传感器将载荷信号通过传输信道传递到电脑接口, 集成到DIC数据采集软件GOM Correlate中. 通过图像捕捉及数据传输, 可以得到拉伸试验中样件的力学响应及局部形变情况.
2.2 参数识别
根据图6所示的PUA树脂在不同拉伸速率下测得的弹性模量, 发现其弹性模量随着拉伸速率的提升而缓慢地增加, 证明了光固化PUA树脂的弹性模量存在率相关性. 为获取更为精准的模型参数, 需要将试验测得不同拉伸应变速率下的弹性模量, 作为本构模型中弹性项的输入. 其中, 弹性项部分基于双参数Mooney-Rivlin模型的应变能密度方程, 其杨氏弹性模量E为[33]
![]()
图 6 试验测得不同拉伸速率下的弹性模量Figure 6. The elastic modulus determined by tests at various tensile rates$$ E = 3G = 6({C_{10}} + {C_{01}}) $$ (28) 式中G为材料的剪切模量, 具体结果如表1所示. 在弹性项识别的基础上, 根据率相关的单向拉伸试验结果, 分别对应变率为10−4, 10−3, 10−2以及10 s−1的试验结果进行拟合, 得到真实应力与真实应变的曲线. 通过调整黏性系数η和分数阶参数$\kappa $, 使拟合曲线最大程度接近于试验结果的力学响应曲线. 通过拟合, 得到应变率分别为10−4, 10−3, 10−2以及10 s−1的分数阶$\kappa $分别为0.60, 0.62, 0.63以及0.58. Jin等[34]的工作证明, 选用固定分数阶导数的阶数对模型参数识别影响并不大, 同时相较于可变的$\kappa $值, 固定的$\kappa $值可以提高模型的拟合效率. 因此, 本文通过对不同拉伸速率的拟合分数阶进行处理, 模型选用$\kappa = 0.61$. 此时黏性系数η作为黏性项中的唯一变量, 通过式(29)计算得到不同加载速率下的最优解, 具体的拟合结果如表2所示.
表 1 本构模型的弹性参数Table 1. Elastic parameters of constitutive model$\dot \varepsilon $/s−1 E/GPa C01, C10/GPa 10−4 2.16 C10 = −4.73, C01 = 5.09 10−3 2.17 C10 = −4.73, C01 = 5.09 10−2 2.22 C10 = −4.73, C01 = 5.10 10−1 2.27 C10 = −4.73, C01 = 5.11 表 2 不同速率下固定阶数的分数阶模型材料参数Table 2. Material parameters of the fractional order model with constant order for different strain rates$\dot \varepsilon $/s−1 η $\kappa $ $R_F^2$ 10−4 68 0.61 0.9877 10−3 65 0.61 0.9912 10−2 63 0.61 0.9848 10−1 66 0.61 0.9822 为对比黏性行为的拟合结果, 基于Mooney-Rivlin模型描述弹性响应, 采用广义Maxwell模型模拟黏性行为. 线性黏弹性模型的本构方程可写成
$$ \sigma \left( {\varepsilon ,t} \right) = {\sigma _0}\left( {\varepsilon ,t} \right) - \int_0^t {g\left( {t - s} \right)} \frac{{{\mathrm{d}}{\sigma _0}\left( {\varepsilon ,\tau } \right)}}{{{\mathrm{d}}\tau }}{\mathrm{d}}s $$ (29) 式中, $ {\sigma _0}\left( {\varepsilon ,t} \right) $为弹性响应部分, $ g\left( t \right) $为衰减函数, 可写成Prony级数的形式
$$ g\left( t \right) = 1 - \sum\limits_{i = 1}^N {\left\{ {{g_i}\left[ {1 - \exp \left( {\frac{{ - t}}{{{\tau _i}}}} \right)} \right]} \right\}} $$ (30) 其中, N为Prony级数项数, gi为无量纲系数, τi为松弛时间, 具体的拟合结果如表3所示.
表 3 不同速率下的Prony系数Table 3. Prony coefficients at different rates$\dot \varepsilon $/s−1 g1 τ1/s g2 τ2/s g3 τ3/s 10−4 0.14 1.09 0.08 6.83 0.78 364.64 10−3 0.11 0.39 0.08 2.53 0.77 129.59 10−2 0.13 0.84 0.47 2.53 — — 10−1 0.97 8.71 — — — — 图7为4种拉伸速率下两种黏弹性本构模型的拟合结果. 其中虚线代表同种拉伸速率下的3个试验结果, 实线为通过参数识别得到的拟合结果, 拟合结果的极限应变值为3次试验测得破坏强度对应应变的平均值. 根据不同拉伸应变率下的真实应力-真实应变试验结果, 评估确定系数$ {R^2} $
![]()
图 7 不同拉伸速率下的分数阶模型拟合曲线Figure 7. Fitting curves of the fractional order model at different tensile strain rates$$ {R^2} = 1 - \frac{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{d_i} - {s_i}} \right)}^2}} }}{{\displaystyle\sum\limits_{i = 1}^n {{{\left( {{d_i} - {{\bar d}_i}} \right)}^2}} }} $$ (31) 式中, n表示数据点个数, di表示各个试验点的真实应力值, $ {\bar d_i} $表示其平均值, si表示模型预测点的真实应力值. 此外, 图7中$ R_F^2 $和$ R_L^2 $分别表示分数阶黏弹性模型的确定系数和线性黏弹性模型的确定系数. 通过对比发现, 两黏弹性模型在保证拟合误差的前提下, 线性黏弹性模型需要更多的拟合参数; 当试验的应变率较小时, 线弹性模型需要更多的拟合参数, 且精度低于分数阶黏弹性模型的拟合结果. 因此, 选用分数阶黏弹性模型将提高数值模型的预测效率.
3. 结果与讨论
结合DIC技术, 分别对PUA样件执行应变率为10−4, 10−3, 10−2以及10 s−1的单向拉伸试验, 不同应变率下的待测样件数量均为3个. 对样件的DIC观测区进行标定, 左右两端分别被标记为line 1和line 2. 根据实验结果, 计算得到试验样件在拉伸速率为10−4, 10−3, 10−2以及 10 s−1的平均极限拉伸长度分别为4.6, 4.3, 4.1和3.8 mm. 图8为当端部位移达到平均极限拉伸长度时, 不同拉伸应变率工况下的数值模型仿真结果与部分样件沿y1方向的位移测量结果. 数值模型中line 1和line 2的位移结果被标识在图中. 通过实验结果发现, 光固化样件的极限伸长率随加载应变率的提高而减少. 图8所示的4个试验样件在4种拉伸速率下line 1处的极限位移分别为1.10, 1.09, 1.07以及1.04 mm; 而line 2处的极限位移分别为3.44, 3.16, 3.09以及2.83 mm. 整体数值计算的结果要略大于试验结果. 通过与试验结果对比, line 1和line 2位置的平均误差分别为2.60%和1.37%, 整体的仿真误差为1.98%.
![]()
图 8 y1方向位移云图试验结果与仿真结果的对比Figure 8. Comparison of experimental results and numerical results with respect to the y1-directional displacement如图9所示, 通过DIC分析得到样件在单向拉伸过程中的横向变形情况. 由于光固化样件存在蠕变变形, 因此样件沿y2方向的位移测量结果远小于沿加载方向的变形量. 试验样件在拉伸速率为10−4, 10−3, 10−2以及 10 s−1的最大横向位移分别为0.054, 0.052, 0.048以及0.047 mm. 如图10所示, 相较于拉伸方向的计算误差, 横向数值结果的预测误差更大, 平均预测误差为6.49%. 由于样件在不同应变率下的横向变形量较小, 预测误差在可接受的数值范围内. 因此, 数值模型可以有效预测不同应变率下光固化树脂的黏弹性力学响应.
![]()
图 9 y2方向位移云图试验结果与仿真结果的对比Figure 9. Comparison of experimental results and numerical results with respect to the y2-directional displacement![]()
图 10 试验样件在测试位置的位移测试结果、仿真结果以及误差Figure 10. Experimental results and numerical results of the y1-directional displacement, as well as the errors of samples at target positions4. 结论
基于有限体积法的固体力学框架, 本文建立三维PUA光固化聚合物的力学模型, 根据试验样件的几何尺寸和边界条件确定数值模型的加载与约束. 通过对位移变量执行二阶Legendre多项式展开, 保留了位移分量的高阶项, 以提高模型的仿真精度. 为提高模型的预测效率, 数值模型引入分数阶有限变形Kelvin-Voigt构建胞元的面应变与面应力之间的本构关系. 通过不同应变率下的单向拉伸试验, 对本构模型中的弹性项和黏性项进行参数识别. 弹性项部分使用不同拉伸速率下的弹性模量作为输入参数; 黏性项选用固定阶的分数阶黏弹性模型以提高拟合的效率. 通过对比PUA聚合物样件在极限载荷下的数值仿真云图与DIC分析结果, 发现整体的仿真结果略大于试验的测量结果, 沿拉伸方向的整体仿真误差为1.98%, 证明三维数值模型具有良好的预测精度. 本文所提出力学模型为进一步探究光固化结构件的非线性力学响应创造了必要的前提条件.
-
![]()
图 2 相邻胞元间面位移的连续性条件
Figure 2. Continuous conditions of surface displacement between adjacent cells
![]()
图 3 沿加载方向相邻胞元之间面牵引力的连续性条件
Figure 3. Continuous conditions of surface traction force between adjacent cells along the loading direction
![]()
图 6 试验测得不同拉伸速率下的弹性模量
Figure 6. The elastic modulus determined by tests at various tensile rates
![]()
图 7 不同拉伸速率下的分数阶模型拟合曲线
Figure 7. Fitting curves of the fractional order model at different tensile strain rates
![]()
图 8 y1方向位移云图试验结果与仿真结果的对比
Figure 8. Comparison of experimental results and numerical results with respect to the y1-directional displacement
![]()
图 9 y2方向位移云图试验结果与仿真结果的对比
Figure 9. Comparison of experimental results and numerical results with respect to the y2-directional displacement
![]()
图 10 试验样件在测试位置的位移测试结果、仿真结果以及误差
Figure 10. Experimental results and numerical results of the y1-directional displacement, as well as the errors of samples at target positions
表 1 本构模型的弹性参数
Table 1 Elastic parameters of constitutive model
$\dot \varepsilon $/s−1 E/GPa C01, C10/GPa 10−4 2.16 C10 = −4.73, C01 = 5.09 10−3 2.17 C10 = −4.73, C01 = 5.09 10−2 2.22 C10 = −4.73, C01 = 5.10 10−1 2.27 C10 = −4.73, C01 = 5.11 
下载: 导出CSV
表 2 不同速率下固定阶数的分数阶模型材料参数
Table 2 Material parameters of the fractional order model with constant order for different strain rates
$\dot \varepsilon $/s−1 η $\kappa $ $R_F^2$ 10−4 68 0.61 0.9877 10−3 65 0.61 0.9912 10−2 63 0.61 0.9848 10−1 66 0.61 0.9822
下载: 导出CSV
表 3 不同速率下的Prony系数
Table 3 Prony coefficients at different rates
$\dot \varepsilon $/s−1 g1 τ1/s g2 τ2/s g3 τ3/s 10−4 0.14 1.09 0.08 6.83 0.78 364.64 10−3 0.11 0.39 0.08 2.53 0.77 129.59 10−2 0.13 0.84 0.47 2.53 — — 10−1 0.97 8.71 — — — —
下载: 导出CSV
-
[1] Bakhshandeh E, Sobhani S, Croutxé-Barghorn C, et al. Siloxane- modified waterborne UV-curable polyurethane acrylate coatings: Chemorheology and viscoelastic analyses. Progress in Organic Coatings, 2021, 158: 106323 doi: 10.1016/j.porgcoat.2021.106323
[2] Hu R, Zhang X, Chen Y, et al. Characterization and prediction of the nonlinear creep behavior of 3D-printed polyurethane acrylate. Additive Manufacturing, 2022, 50: 102583 doi: 10.1016/j.addma.2021.102583
[3] 张毅, 薛世峰, 韩丽美等. 半结晶聚合物损伤演化的试验表征与数值模拟. 力学学报, 2021, 53(6): 1671-1683 (Zhang Yi, Xue Shifeng, Han Limei, et al. Experimental characterization and numerical simulation of damage evolution in semi-crystalline. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(6): 1671-1683 (in Chinese) Zhang Yi, Xue Shifeng, Han Limei, et al. Experimental characterization and numerical simulation of damage evolution in semi-crystalline. Chinese Journal of Theoretical and Applied Mechanics, 2021, 53(6): 1671-1683 (in Chinese)
[4] 田传帅, 詹林, 肖锐. 基于超弹性模型的玻璃态聚合物应变强化行为研究. 力学学报, 2023, 55(7): 1473-1483 (Tian Chuanshuai, Zhan Lin, Xiao Rui. Modelling strain hardening of glassy polymers based on hyperelastic models. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(7): 1473-1483 (in Chinese) Tian Chuanshuai, Zhan Lin, Xiao Rui. Modelling strain hardening of glassy polymers based on hyperelastic models. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(7): 1473-1483 (in Chinese)
[5] Hossain M, Navaratne R, Perić D. 3D printed elastomeric polyurethane: Viscoelastic experimental characterizations and constitutive modelling with nonlinear viscosity functions. International Journal of Non-Linear Mechanics, 2020, 126: 103546 doi: 10.1016/j.ijnonlinmec.2020.103546
[6] Du Z, Yang Y, Wang Z, et al. A finite strain visco-hyperelastic damage model for rubber-like materials: Theory and numerical implementation. Acta Mechanica Sinica, 2023, 39(3): 222473 doi: 10.1007/s10409-023-22473-x
[7] Ho CH, Romero P. Characterizing the low-temperature viscoelastic behavior of asphalt mixtures: A comparative study. International Journal of Pavement Research and Technology, 2013, 6(5): 479
[8] Alizadeh N, Celestine AN, Auad ML, et al. Mechanical characterization and modeling stress relaxation behavior of acrylic-polyurethane-based graft-interpenetrating polymer networks. Polymer Engineering & Science, 2021, 61(5): 1299-1309
[9] Meng R, Yin D, Drapaca CS. A variable order fractional constitutive model of the viscoelastic behavior of polymers. International Journal of Non-Linear Mechanics, 2019, 113: 171-177 doi: 10.1016/j.ijnonlinmec.2019.04.002
[10] Shariyat M, Mohammadjani R. 3D nonlinear variable strain-rate-dependent-order fractional thermoviscoelastic dynamic stress investigation and vibration of thick transversely graded rotating annular plates/discs. Applied Mathematical Modelling, 2020, 84: 287-323 doi: 10.1016/j.apm.2020.03.023
[11] Burlon A, Alotta G, Di Paola M, et al. An original perspective on variable-order fractional operators for viscoelastic materials. Meccanica, 2021, 56: 769-784 doi: 10.1007/s11012-021-01316-4
[12] Meng R, Cao L, Zhang Q. Study on the performance of variable-order fractional viscoelastic models to the order function parameters. Applied Mathematical Modelling, 2023, 121: 430-444 doi: 10.1016/j.apm.2023.05.017
[13] Adolfsson K, Enelund M, Olsson P. On the fractional order model of viscoelasticity. Mechanics of Time-dependent Materials, 2005, 9: 15-34 doi: 10.1007/s11043-005-3442-1
[14] 顾建平, 孙慧玉, 方建士等. 热致非晶态形状记忆聚合物的热黏弹性参数及回复行为. 高分子材料科学与工程, 2016, 32(6): 107-112 (Gu Jianping, Sun Huiyu, Fang Jianshi, et al. Thermovisoelastic parameters and recovery behaviors of thermally activated amorphous shape memory polymers. Polymer Materials Science & Engineering, 2016, 32(6): 107-112 (in Chinese) Gu Jianping, Sun Huiyu, Fang Jianshi, et al. Thermovisoelastic parameters and recovery behaviors of thermally activated amorphous shape memory polymers. Polymer Materials Science & Engineering, 2016, 32(6): 107-112 (in Chinese)
[15] Drozdov AD. Fractional differential models in finite viscoelasticity. Acta Mechanica, 1997, 124(1): 155-180
[16] Sumelka W, Voyiadjis GZ. A hyperelastic fractional damage material model with memory. International Journal of Solids and Structures, 2017, 124: 151-160 doi: 10.1016/j.ijsolstr.2017.06.024
[17] 尹耀得, 赵德敏, 刘建林等. 丙烯酸弹性体的率相关分数阶黏弹性模型研究. 力学学报, 2022, 54(1): 154-162 (Yin Yaode, Zhao Demin, Liu Jianlin, et al. Study on The Rate Dependency of Acrylic Elastomer Based Fractional viscoelastic model. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 154-162 (in Chinese) Yin Yaode, Zhao Demin, Liu Jianlin, et al. Study on The Rate Dependency of Acrylic Elastomer Based Fractional viscoelastic model. Chinese Journal of Theoretical and Applied Mechanics, 2022, 54(1): 154-162 (in Chinese)
[18] Taylor GA, Bailey C, Cross M. Solution of the elastic/visco-plastic constitutive equations: A finite volume approach. Applied Mathematical Modelling, 1995, 19(12): 746-760 doi: 10.1016/0307-904X(95)00093-Y
[19] Demirdzic I, Martinovic P, Ivankovic A. Numerical simulation of thermal deformation in welded workpiece. Zavarivanje, 1988, 31(5): 209-219
[20] Cavalcante M, Pindera MJ, Khatam H. Finite-volume micromechanics of periodic materials: Past, present and future. Composites Part B: Engineering, 2012, 43(6): 2521-2543 doi: 10.1016/j.compositesb.2012.02.006
[21] Sevilla R, Giacomini M, Huerta A. A face-centred finite volume method for second-order elliptic problems. International Journal for Numerical Methods in Engineering, 2018, 115(8): 986-1014 doi: 10.1002/nme.5833
[22] Ivanova EA, Jatar Montaño LE. A new approach to solving the solid mechanics problems with matter supply. Continuum Mechanics and Thermodynamics, 2021, 33(4): 1829-1855 doi: 10.1007/s00161-021-01014-2
[23] Hassan OI, Ghavamian A, Lee CH, et al. An upwind vertex centred finite volume algorithm for nearly and truly incompressible explicit fast solid dynamic applications: Total and updated lagrangian formulations. Journal of Computational Physics: X, 2019, 3: 100025
[24] Cardiff P, Tuković Ž, Jasak H, et al. A block-coupled finite volume methodology for linear elasticity and unstructured meshes. Computers & Structures, 2016, 175: 100-122
[25] Cardiff P, Demirdžić I. Thirty years of the finite volume method for solid mechanics. Archives of Computational Methods in Engineering, 2021, 28(5): 3721-3780 doi: 10.1007/s11831-020-09523-0
[26] Scolaro A, Fiorina C, Clifford I, et al. Development of a semi-implicit contact methodology for finite volume stress solvers. International Journal for Numerical Methods in Engineering, 2022, 123(2): 309-338 doi: 10.1002/nme.6857
[27] Soar P, Kao A, Pericleous K. The impact of two and three dimensional assumptions on coupled structural mechanics and microstructure solidification modelling. Crystals, MDPI, 2023, 13(1): 114 doi: 10.3390/cryst13010114
[28] Castrillo P, Schillaci E, Rigola J. High-order cell-centered finite volume method for solid dynamics on unstructured meshes. Computers & Structures, 2024, 295: 107288
[29] Cai H, Ye J, Shi J, et al. A new two-step modeling strategy for random micro-fiber reinforced composites with consideration of primary pores. Composites Science and Technology, 2022, 218: 109122 doi: 10.1016/j.compscitech.2021.109122
[30] Askarian AR, Permoon MR, Shakouri M. Vibration analysis of pipes conveying fluid resting on a fractional Kelvin-Voigt viscoelastic foundation with general boundary conditions. International Journal of Mechanical Sciences, 2020, 179: 105702 doi: 10.1016/j.ijmecsci.2020.105702
[31] Yao Q, Dong P, Zhao Z, et al. Temperature dependent tensile fracture strength model of rubber materials based on Mooney-Rivlin model. Engineering Fracture Mechanics, 2023, 292: 109646 doi: 10.1016/j.engfracmech.2023.109646
[32] Braconnier DJ, Jensen RE, Peterson AM. Processing parameter correlations in material extrusion additive manufacturing. Additive Manufacturing, 2020, 31: 100924 doi: 10.1016/j.addma.2019.100924
[33] 马建章. 基于 Workbench 的导电橡胶 Mooney-Rivlin 参数拟合与应用. 无线电工程, 2017, 47(10): 79-82 (Ma Jianzhang. Mooney-rivlin parameter fitting and application of conductive rubber based on workbench. Radio Engineering, 2017, 47(10): 79-82 (in Chinese) Ma Jianzhang. Mooney-rivlin parameter fitting and application of conductive rubber based on workbench. Radio Engineering, 2017, 47(10): 79-82 (in Chinese)
[34] Jin L, Zhao D, Liu J. A visco-hyperelastic constitutive model for rubber considering the strain level and one case study in the sealing packer. Acta Mechanica Solida Sinica, 2023, 36(5): 710-723 doi: 10.1007/s10338-023-00397-w
计量
- 文章访问数: 90
- HTML全文浏览量: 28
- PDF下载量: 19
