介电弹性体（dielectric elastomers, DEs）是典型的电致活性聚合物，因其在电场激励下能够产生显著变形而备受关注. 典型的DE结构由1层柔软的弹性聚合物薄膜（如硅橡胶、丙烯酸酯）夹在2层柔性电极（如碳纳米管、银纳米线或导电聚合物）之间构成. 当在电极间施加外部电压时，正负电荷在电极表面积聚，产生强大的静电吸引力，即麦克斯韦应力（Maxwell stress）. 该应力在厚度方向表现为压应力，导致DEs薄膜在厚度方向收缩，同时在平面方向产生双轴变形，产生显著的面积膨胀应变，从而实现大变形驱动^[1-2]. 这种独特的力电转换特性使DEs具有响应速度快（可达毫秒级）、能量密度高、轻量化、低噪音运行、优异的生物相容性以及可承受大应变和适应复杂曲面等显著优势. 基于这些特性，DEs在软体机器人、仿生设备、可变性光学器件、触觉反馈装置、能量收集以及生物医学工程领域展现出极其广阔的应用场景，尤其适用于需要柔性、适应性强且高效驱动的场景^[3-4]. 然而，目前该技术大多仍局限于实验室应用阶段，在实际工程场景中，DEs 技术的落地还面临诸多挑战，例如驱动电压高（高达数千伏）、材料性能与稳定性欠佳（存在起皱、蠕变、分层等问题），以及非线性控制难度较大等.

为了深入理解并预测DEs在力电耦合场作用下的复杂力学行为，建立本构模型至关重要. 在DE理论建模领域，本构模型的研究始于Pelrine等^[5]的开创性工作. 随后，Suo等^[6]基于能量原理，系统性地建立了介电弹性体的非线性理论框架. Henann等^[7]在此基础上进一步验证了该框架下介电弹性体力电耦合理论的合理性，并推导了相应的应力、电场及电位移表达式，为后续的数值模拟与工程应用奠定了理论基础. 此后，众多学者致力于发展各类介电弹性体本构模型，用以刻画其电致失稳、大变形响应、非线性硬化等关键力学现象^[8-9]. 这类模型^[10-26]通常将描述材料弹性行为的超弹性模型（如neo-Hookean、Mooney-Rivlin、Ogden模型）与理想介电响应理论相耦合. 近年来，相关学者^[27-38]提出更为完善的模型，例如考虑材料黏弹性效应、预拉伸影响以及失效机理的模型，以期更全面地表征DEs的实际性能. 如Xiao等^[23]提出的本构模型，能够成功描述双峰网络介电弹性体的特性^[24]. 该领域研究得到了广泛关注，新理论与新模型不断涌现.

在DEs的本构模型中，精准表征弹性体力学行为是构建可靠力电耦合模型的前提. 目前广泛采用的超弹性力学模型，如基于网络统计理论的八链模型和Gent模型，在描述单一加载模式（如单轴拉伸）下的超弹性行为时表现良好. 然而，这些经典模型在同时预测材料不同加载条件（如单轴拉伸、等双轴拉伸）下的实验数据时，往往存在明显的不足，特别是难以精确描述复杂多轴变形下的应力-应变响应^[8-9]，而这对于介电弹性体性能准确预测至关重要. 上述局限性的根源在于，这些模型主要关注由化学交联点决定的网络弹性，而未能充分计入聚合物链之间缠结所产生的拓扑约束效应^[19]. 这种缠结效应在高分子材料，尤其是未充分交联或处于大变形状态的弹性体中，对力学行为（如模量硬化、变形模式依赖性）起着至关重要的作用. 为了解决这一问题，相关学者在Doi和Edwards^[17]发展的经典管模型（Tube Model）基础上，发展了一系列考虑缠结效应的本构模型. 这些模型通过将系统的总自由能分解为交联和缠结2部分，显著提升了对聚合物材料复杂多轴力学行为的预测能力.

基于上述分析，现有介电弹性体本构模型因忽略缠结效应，难以准确预测材料在复杂多轴加载下的力学响应及电致变形行为. 为此，本研究基于Doi-Edwards理论，在力电耦合框架中引入缠结效应，通过在超弹性自由能中同时考虑交联与缠结贡献，构建了一个介电弹性体力-电耦合本构模型，旨在更精准地预测材料在力电载荷下的响应行为，为高性能驱动器和软体机器人设计提供理论支撑.

当前广泛使用的基于分子链统计理论的超弹性模型，包括高斯链模型（即neo-Hookean模型），James-Guth 三链模型、Arruda-Boyce 八链模型以及全链网络模型，建模均基于统计力学，从单根分子链的自由能密度函数出发，并结合不同的宏微观映射，建立本构理论.

上述模型仅包括交联部分对于自由能的贡献，其中单链模型采用Kuhn等^[16]提出的Langevin统计模型，其形式可以表示为

(1)$ {W}_{\text{c}}=kTN\left(\dfrac{\lambda }{\sqrt{N}}\beta +\ln\;\left(\dfrac{\beta }{\text{sinh}\;\;\beta }\right)\right)+\text{const} , $

(2)$ \beta ={{L}}^{-1}\left({\lambda }/{\sqrt{N}}\right) . $

式中：${W}_{\text{c}} $为交联部分自由能密度，$ k $为玻尔兹曼常数，$ T $为绝对温度，$ N $为链段数，$ \lambda $为拉伸比，$ \text{const} $为确保初始无应力状态下自由能为零的常数， $ {{L}}^{-1} $为Langevin逆函数. 其中，拉伸比$ \lambda $与链的端到端距离$ r $有关，可表示为$ \lambda =r/(\sqrt{N}b) $，$ b $为库恩长度.

交联部分的自由能密度由拉伸过程中交联聚合物链的构象变化所决定. 在传统模型中，通常假定模型中的分子链的变形不受周围分子链的影响. 但实际上，链的变形会受到周围分子链缠结的约束效应. Doi等^[17]采用管模型来描述缠结效应，链的运动被限制在管内部，因此缠结效应随着管的收缩而增加，其假设单根分子链的运动都被限制在管内，当管道收缩时，缠结部分增加，缠结自由能密度W_e与管截面积收缩比呈线性关系：

(3)$ {W}_{\text{e}}=\alpha kTN{\left(\dfrac{b}{{d}_{0}}\right)}^{2}v +\text{const} . $

式中：$ \alpha $为恒定的预因子，$ {d}_{0} $为管的初始直径，$ v $为变形后管道与变形前管道的面积之比.

建立微宏观关联的本构模型，关键在于确定聚合物链/管道模型的微观形变与宏观变形之间的关联关系. Khiem等^[18]提出结合微球模型和仿射假说的网络平均模型，以获得宏观和微观变形之间的关系. 假设初始未变形的分子链均匀分布在微球表面^[19]. 如图1所示，变形后，微球沿主应变方向转变为椭球体. 分子链在初始构型下的单位矢量$ {\boldsymbol{R}}_{{0}} $和当前构型下的$ \boldsymbol{R} $可建立如下关系：

(4)$ \boldsymbol{R}=\boldsymbol{F}{\boldsymbol{R}}_{{0}} . $

式中：$ \boldsymbol{F} $为变形梯度张量.

由于微球初始分子链在表面上近似均匀分布并且微观尺度远小于宏观变形尺度，假设微球表面所有分子链长度的均方根值可作为宏观变形条件下分子链的等效长度$\bar{\lambda} $^[19]：

(5)$ \left\{\begin{split}\bar{\lambda} =&\sqrt{\oiint \lambda^2 \mathrm{d} S}=\sqrt{\frac{1}{4 {\text{π}}} \int_0^{2 {\text{π}}} \int_0^{\text{π}} l \sin \;\varphi \;\mathrm{d} \varphi \mathrm{d} \theta}= \\& \sqrt{\frac{\lambda_1^2+\lambda_2^2+\lambda_3^2}{3}} = \sqrt{\frac{I_1}{3}} ,\\ l= &\lambda_1^2(\sin \;\varphi \cos \;\theta)^2 + \lambda_2^2(\sin \;\varphi \sin \;\theta)^2 + \lambda_3^2(\cos \;\varphi)^2.\end{split}\right. $

式中：θ为初始构型中分子链在XOY平面的投影与X轴的夹角，$\varphi $为分子链与Z轴的夹角，λ₁、λ₂、λ₃为3个主方向上的拉伸比.

同理，基于Nanson公式，参考构型下的面积矢量和当前构型下的面积矢量可建立如下关系：

(6)$ \text{d}\boldsymbol{S}=J{\boldsymbol{F}}^{{-{\mathrm{T}}}}\text{d}{\boldsymbol{S}}_{{0}} . $

式中：$ J=\mathrm{det}\;\boldsymbol{F} $，通常假设体积不可压缩，故$ J=1 $；上标$ {}^{{-{\mathrm{T}}}} $表示逆的转置运算. 定义面积之比$v = \text{d}S/\text{d}{S}_{0} $. 假设变形后微球所有单元的面积与初始面积比的均方根值$\overline{v } $为宏观变形管的面积收缩之比：

(7)$ \left.\begin{split}\overline{v }=& \sqrt{{\int}_{\mathrm{d}S}{v}^{2}\mathrm{d}S}=\sqrt{\frac{1}{4{\text{π}} }\int\limits_{0}^{2{\text{π}} }\int\limits_{0}^{{\text{π}} }\eta \sin\; \varphi \mathrm{d}\varphi \mathrm{d}\theta }=\\& \sqrt{\frac{\lambda _{2}^{2}\lambda _{3}^{2}+\lambda _{3}^{2}\lambda _{1}^{2}+\lambda _{1}^{2}\lambda _{2}^{2}}{3}}=\sqrt{\frac{{I}_{2}}{3}},\\\eta =& \lambda _{2}^{2}\lambda _{3}^{2}{\left(\sin\; \varphi \cos\; \theta \right)}^{2} + \lambda _{3}^{2}\lambda _{1}^{2}{\left(\sin\; \varphi \sin\; \theta \right)}^{2} +\\&\lambda _{1}^{2}\lambda _{2}^{2}{\left(\cos\; \varphi \right)}^{2}\end{split} \right\} $

结合单链模型和微观-宏观映射关系，可以得到宏观力学本构模型. 假设单位体积内分子链数目为$ n $，并使用平均化的$ \overline{\lambda } $、$ \overline{\nu } $代替$ \lambda $、 $ \nu $，可以得到

(8)$ W={W}_{\text{c}}+{W}_{\text{e}}+\text{const} , $

(9)$ {W}_{\text{c}}={G}_{\text{c}}N\left(\frac{\sqrt{{I}_{1}/3}}{\sqrt{N}}\beta +\ln\; \frac{\beta }{\text{sinh}\;\beta }\right) , $

(10)$ {W}_{\text{e}}=3{G}_{\text{e}}\sqrt{\dfrac{{I}_{2}}{3}}+\text{const}. $

式中:$ W $为总自由能密度； $ {W}_{\text{c}} $为交联部分的自由能密度，描述交联点间分子链段构象熵引起的弹性形变；$ {W}_{\text{e}} $为分子链缠结约束引起的自由能；$ \text{const} $为积分常数，其值的确定须满足初始未变形状态下材料的应变能密度为零的条件；$ {{G}}_{\text{c}}=nkT $为交联部分的剪切模量；$ {{G}}_{\text{e}}=\alpha nkT{\left(\dfrac{b}{{d}_{0}}\right)}^{2}\bigg/3 $为缠结部分的剪切模量.

须指出的是，部分学者也提出了其他类型的缠结效应模型，较有代表性的包括Davidson等^[25]所建立的基于非仿射网络假设的本构模型，以及 Xiang等^[26]基于三链模型所建立的缠结模型. 不同的缠结形式在模型拟合单双轴以及预测力电耦合效应时，其本构模型的预测效果相近，因此本研究采用上述平均网络模型.

在电场作用下，介电弹性体会发生显著的电致变形行为（见图2）. 自由能包括变形能和电场能. 本研究采用理想介电模型来表示电场自由能，并考虑不可压缩条件，则系统自由能可以表示为

(11)$ W={W}_{\text{s}}+{W}_{\text{elec}}-p\left(J-1\right) , $

(12)$ {W}_{\text{s}} = {G}_{\text{c}}N\left( \frac{\sqrt{{I}_{1}/3}}{\sqrt{N}}\beta + \ln \;\frac{\beta }{\text{sinh}\;\beta } \right) + 3{G}_{\text{e}}\sqrt{\frac{{I}_{2}}{3}} + \text{const}, $

(13)$ {W}_{\text{elec}}=\frac{1}{2\varepsilon }{J}^{-1}\tilde{\boldsymbol{D}}\boldsymbol{C}\tilde{\boldsymbol{D}} , $

(14)$ {I}_{1}={\mathrm{tr}}\;(\boldsymbol{C}),{I}_{2}=\frac{1}{2}[{({\mathrm{tr}}\;\boldsymbol{C})}^{2}-{\mathrm{tr}}\;({\boldsymbol{C}}^{2})] . $

式中：$ {W}_{\text{s}} $为变形贡献的自由能密度（其具体形式与上章节的超弹性模型一致），$ {W}_{\text{elec}} $为电场贡献的自由能密度，$ p $为不可压缩条件参数，$ \varepsilon $为绝对介电常数，$ \boldsymbol{C}\mathbf{=}{\boldsymbol{F}}^{\mathrm{T}}\boldsymbol{F} $为右柯西-格林应变张量，$ \tilde{\boldsymbol{D}} $为名义电位移矢量.

在连续介质力学框架下，介电弹性体的热力学第二定律可表示为

(15)$ \dot{W}-\boldsymbol{S}\colon \frac{1}{2}\dot{\boldsymbol{C}}\mathbf-\tilde{\boldsymbol{E}}\colon \dot{\tilde{\boldsymbol{D}}}\leqslant 0 . $

式中：$ \boldsymbol{S} $为第二Piola-Kirchhoff应力，$ \tilde{\boldsymbol{E}} $为名义电场. 自由能函数对时间的导数$ \dot{W} $可以在链式法则下表达为

(16)$ W=W(\boldsymbol{C},\tilde{\boldsymbol{D}}) , $

(17)$ \dot{W}=\frac{\partial W}{\partial \boldsymbol{C}}\colon \dot{\boldsymbol{C}}+\frac{\partial W}{\partial \tilde{\boldsymbol{D}}}\colon \dot{\tilde{\boldsymbol{D}}} . $

将式(17)代入到式(15)中，得到

(18)$ \left[2 \frac{\partial W}{\partial \boldsymbol{C}}-\boldsymbol{S}\right]: \frac{1}{2} \dot{\boldsymbol{C}}+\left[\frac{\partial W}{\partial \tilde{\boldsymbol{D}}}-\tilde{\boldsymbol{E}}\right]: \dot{\tilde{\boldsymbol{D}}} \leqslant 0 . $

热力学框架要求式(18)恒成立，则须满足以下关系：

(19)$ \boldsymbol{S}=2\frac{\partial W}{\partial \boldsymbol{C}},\;\tilde{\boldsymbol{E}}=\frac{\partial W}{\partial \tilde{\boldsymbol{D}}} . $

相应的第一Piola-Kirchhoff应力和Cauchy应力可表示为

(20)$ \boldsymbol{P}=\boldsymbol{F}\boldsymbol{S},\;{{\bf\textit{σ}} }=\frac{1}{J}\boldsymbol{F}\boldsymbol{S}{\boldsymbol{F}}^{\rm{T}} . $

将式(11)~ (13)代入式(19)~(20)，得到

(21)$ \begin{split} \boldsymbol{S} =&2 G_{\mathrm{c}} \sqrt{\frac{N}{3 I_1}} \beta \boldsymbol{I}+2 G_{\mathrm{e}} \frac{I_1 \boldsymbol{I}-\boldsymbol{C}}{\sqrt{I_2 / 3}} +\\& \frac{1}{\varepsilon} J^{-1} \tilde{\boldsymbol{D}} \otimes \tilde{\boldsymbol{D}}-J p \boldsymbol{C}^{-1},\end{split} $

(22)$ \begin{split} \boldsymbol{P} =&2 G_{\mathrm{c}} \sqrt{\frac{N}{3 I_1}} \beta \boldsymbol{F}+2 G_{\mathrm{e}} \frac{I_1 \boldsymbol{F}-\boldsymbol{F} \boldsymbol{F}^{\rm{T}} \boldsymbol{F}}{\sqrt{I_2 / 3}} +\\& \frac{1}{\varepsilon} J^{-1} \boldsymbol{F} \tilde{\boldsymbol{D}} \otimes \tilde{\boldsymbol{D}}-J p \boldsymbol{F}^{\rm{T}},\end{split} $

(23)$ \begin{split} \boldsymbol{\sigma}= & 2 J^{-1} G_{\mathrm{c}} \sqrt{\frac{N}{3 I_1}} \beta \boldsymbol{b}+2 J^{-1} G_{\mathrm{e}} \frac{I_1 \boldsymbol{b}-\boldsymbol{b}^2}{\sqrt{I_2 / 3}}+ \\& \frac{1}{\varepsilon} J^{-1} \boldsymbol{F} \tilde{\boldsymbol{D}} \otimes J^{-1} \boldsymbol{F} \tilde{\boldsymbol{D}}-p \boldsymbol{I},\end{split} $

(24)$ \tilde{\boldsymbol{E}}=\frac{1}{\varepsilon }{J}^{-1}\boldsymbol{C}\tilde{\boldsymbol{D}} . $

在不可压材料的双轴拉伸（沿$ {\lambda }_{1}、{\lambda }_{2} $方向）实验中，有

(25)$ \boldsymbol{F}={\lambda }_{1}{\boldsymbol{e}}_{1}\otimes {\boldsymbol{e}}_{1}+{\lambda }_{2}{\boldsymbol{e}}_{2}\otimes {\boldsymbol{e}}_{2}+\frac{1}{{\lambda }_{1}{\lambda }_{2}}{\boldsymbol{e}}_{3}\otimes {\boldsymbol{e}}_{3} , $

(26)$ \tilde{\boldsymbol{D}}=\varepsilon \left(\mathit{\Phi }/{L}_{3}\right)\lambda _{1}^{2}\lambda _{2}^{2}{\boldsymbol{e}}_{3} . $

将式(25)、(26)代入式(22)，退化为双轴拉伸模式下的名义应力表达式，其中沿$ {\lambda }_{1} $方向的名义应力可以表达为

(27)$ {P}_{1}=\frac{\partial {W}_{\text{s}}}{\partial {\lambda }_{1}}-\varepsilon {\left(\frac{\mathit{\Phi }}{{L}_{3}}\right)}^{2}{\lambda }_{1}\lambda _{2}^{2} . $

电压（$ \mathit{\Phi } $）和拉伸比（$ \lambda $）之间的关系如下：

(28)$ \mathit{\Phi }={L}_{3}\times \sqrt{\left({\dfrac{\partial {W}_{\text{s}}}{\partial {\lambda }_{1}}-{P}_{1}}\right)\bigg/\left({\varepsilon {\lambda }_{1}\lambda _{2}^{2}}\right)} . $

综上所述，基于考虑缠结效应的本构模型推导得到对应的第一、第二Piola-Kirchhoff应力及Cauchy应力的张量表达式. 将其代入双轴拉伸工况后，进一步获得电压与拉伸比的关系，可进一步验证本构模型. 其中令模型中缠结部分的剪切模量$ {G}_{\text{e}}=0 $，力学模型退化为八链模型，与考虑缠结效应的模型进行分析对比.

为了验证所建模型的可靠性，采用上述本构理论对文献[12]中的聚丙烯酸酯弹性体（3M^TM VHB4905）力电耦合效应进行预测. 采用文献[20]中关于该材料的单轴、纯剪切、等双轴力学加载的实验数据来验证力学模型的准确性. 须指出的是，文献[20]中的力电耦合测试采用了分级加载-弛豫的方式，以确保材料达到热力学平衡态；而基础力学性能测试则在低应变率（$ 0.003\;{\mathrm{s}}^{-1} $）下进行. 尽管上述2种加载模式在黏性耗散的衰减程度上存在差异，但均严格满足准静态假设.

为了量化模型预测结果与实验数据的吻合程度，采用非线性最小二乘法对模型参数进行迭代优化，目标函数为实验应力与预测应力之差的平方和最小化.

对于沿$ {\lambda }_{1} $方向的单轴拉伸实验，假设材料不可压缩，横向变形由不可压缩条件决定. 此时变形梯度张量$ \boldsymbol{F} $及其对应的柯西应力张量${\boldsymbol{ \sigma }}$可以表示为

(29)$ \boldsymbol{F}={\lambda }_{1}{\boldsymbol{e}}_{1}\otimes {\boldsymbol{e}}_{1}+\frac{1}{\sqrt{{\lambda }_{1}}}{\boldsymbol{e}}_{2}\otimes {\boldsymbol{e}}_{2}+\frac{1}{\sqrt{{\lambda }_{1}}}{\boldsymbol{e}}_{3}\otimes {\boldsymbol{e}}_{3} , $

(30)$ {\sigma }_{1} = \frac{{G}_{\text{c}}}{3}\frac{\sqrt{N}}{\sqrt{\left({I}_{1}/3\right)}}\beta \left(\lambda _{1}^{2} - \frac{2}{{\lambda }_{1}}\right) + \frac{{G}_{\text{e}}}{\sqrt{{I}_{2}/3}}\left({\lambda }_{1} - \frac{1}{\lambda _{1}^{2}}\right). $

针对沿$ {\lambda }_{1} $和$ {\lambda }_{2} $方向的双轴拉伸实验，变形梯度张量$ \boldsymbol{F} $对应式(25)，与其对应的柯西应力张量$ \boldsymbol{\sigma } $在$ {\lambda }_{1} $方向上的分量可以表示为

(31)$ {\sigma }_{1} = \dfrac{{G}_{\text{c}}}{3}\dfrac{\sqrt{N}}{\sqrt{{I}_{1}/3}}\beta \left( \lambda _{1}^{2} - \dfrac{1}{\lambda _{1}^{2}\lambda _{2}^{2}} \right) + \dfrac{{G}_{\text{e}}}{\sqrt{{I}_{2}/3}}\lambda _{2}^{2}\left( \lambda _{1}^{2} - \frac{1}{\lambda _{1}^{2}\lambda _{2}^{2}} \right) . $

分别采用八链模型和网络平均模型对聚丙烯酸酯弹性体（3M^TM VHB 4905）的单轴拉伸试验数据进行拟合，结果如图3所示. 其中，$ p_{1}^{\exp } $为实验测得的第1主方向（加载方向）名义应力；$ p_{1}^{\text{NAT}} $和$ p_{1}^{\text{8-chain}} $分别为网络平均管模型（NAT）与八链模型预测的第 1 主方向名义应力，其中网络平均管模型可进一步分解为交联网络贡献（$ p_{1}^{\text{c}} $）与缠结网络贡献（$ p_{1}^{\text{e}} $）；$ \Delta {p} $为模型预测的误差带. 可以看出，在大变形区域（$ \lambda \gt 4 $），2种模型均能较好地拟合单轴拉伸的曲线. 然而在小变形时，八链模型的预测值显著低于实验数据，未能准确反映材料的初始模量特征. 考虑缠结效应的本构模型相比于传统的八链模型具有更高的拟合精度.

为了确保本构模型在多轴复杂应力状态下的预测精度，采用联合参数优化策略. 首先基于单轴拉伸数据确定模型参数的初始范围，随后引入等双轴拉伸数据，构建包含2种加载模式的联合误差目标函数.

参数优化的联合误差目标函数定义如下：

(32)$ {R}=\dfrac{1}{n}\sum\limits_{i=1}^{n}\dfrac{\left| {\sigma }_{\mathrm{Exp}}-{\sigma }_{\mathrm{Sim}}\right| }{{\sigma }_{\mathrm{Exp}}} . $

式中：$ n $为数据点数量，$ {\sigma }_{\mathrm{Exp}} $为实验测得的柯西应力，$ {\sigma }_{\mathrm{Sim}} $为预测的柯西应力. ${R} $的全局最小值对应最优材料参数.

基于上述参数拟合准则，选取一组最优参数，使模型能够同时较好地拟合单轴拉伸与等双轴拉伸实验数据. 图4展示了本研究所建立的模型对聚丙烯酸酯材料的单轴拉伸（见图4（a））和等双轴拉伸（见图4（b））试验的拟合曲线. 结果表明，当使用一组参数同时对单轴与等双轴实验数据进行拟合时，八链模型的预测曲线显著低于实验数据点. 这一结果证明了仅依赖第一应变不变量$ {I}_{1} $的经典八链模型在描述多轴变形时具有局限性，即无法通过单一参数组同时预测单轴与双轴的应力响应. 本研究建立的考虑缠结效应的模型（红色实线）在引入第二不变量 $ {I}_{2} $ 相关的缠结约束效应后，无论是在单轴大变形的硬化阶段还是在等双轴的高应力区域，均展现出较好的预测精度.

为了探究缠结效应在模型拟合中的作用机制，进一步评估交联部分和缠结部分对材料应力响应的贡献，结果如图5所示. 可以看出，在单轴拉伸模式中，小变形和中等应变阶段（$ \lambda \lt 1.5 $），缠结网络贡献占比超过40%，对力学响应有显著影响. 这一现象与图3的结论一致，即在小应变和中等应变时，考虑缠结效应的本构模型性能优于八链模型. 在多轴拉伸状态下，缠结部分与交联部分对于整体应力的贡献大致相当，表明交联、缠结的耦合机制是能准确预测软材料力学响应的关键因素. 忽略缠结效应的贡献将导致本构模型的拟合参数出现系统性偏差.

如表1所示为分子统计理论模型的材料参数及单轴拉伸拟合误差. 可以看出，八链模型的链段参数$ N $在多轴拟合和单轴拟合时，出现较大差异. 这主要是在本构模型中，缠结效应的贡献被强制归入交联项. 相较之下，考虑了缠结效应的模型无论是在单轴拟合还是在同时拟合多轴拉伸试验中，均表现出显著的参数稳定性.

Mathew等^[12]对厚度为0.5 mm的聚丙烯酸酯弹性体VHB 4905开展了纯剪切模式下力电驱动的试验. 其实验设计通过约束实现面内预拉伸，固定横向位移边界条件，在纵向方向上施加恒定载荷实现预拉伸，如图6所示. 通过离散方式增加电压来获得实验数据点，保持电压5 min后记录相应实验数据，从而确保在准静态条件下测量不同预拉伸比下（$ {\lambda }_{\text{1p}}= 0.5、0.9、1.8、3.0、3.8 $，$ {\lambda }_{\text{2p}}=6.0 $）的驱动电压-拉伸比的关系. 本研究采用的力电耦合数据集源于此试验体系.

考虑纯剪切实验（沿$ {\lambda }_{1} $方向拉伸，$ {\lambda }_{2} $方向固定）. 同样假定材料不可压缩，则在纯剪切模式下，变形梯度张量及力电响应关系可表示为

(33)$ \boldsymbol{F}={\lambda }_{1}{\boldsymbol{e}}_{1}\otimes {\boldsymbol{e}}_{1}+{\lambda }_{2\text{p}}{\boldsymbol{e}}_{2}\otimes {\boldsymbol{e}}_{2}+\dfrac{1}{{\lambda }_{1}{\lambda }_{2\text{p}}}{\boldsymbol{e}}_{3}\otimes {\boldsymbol{e}}_{3} , $

(34)$ {P}_{\text{1p}}={\dfrac{\partial {W}_{\text{s}}}{\partial {\lambda }_{\text{1p}}}}_{{{\lambda }_{2}}={{\lambda }_{\text{2p}}}}-\varepsilon {\left(\dfrac{\mathit{\Phi }}{{L}_{3}}\right)}^{2}{\lambda }_{\text{1p}}\lambda _{2\text{p}}^{2} , $

(35)$ \varPhi=L_3 \times \sqrt{\dfrac{\dfrac{\partial W_{\mathrm{s}}}{\partial \lambda_1}_{ \lambda_2=\lambda_{2 \mathrm{p}}}-P_{1 \mathrm{p}}}{\varepsilon \lambda_1 \lambda_{2 \mathrm{p}}^2}} . $

基于理想介电理论框架，文献[12]在EV模型中参考Jiang等 ^[21]的工作直接采用经验性相对介电常数$ {\varepsilon }_{\text{r}}=3 $. 在本工作中，选取主方向上不同预拉伸比实验曲线中的一条为基准，通过该曲线拟合标定介电常数$ {\varepsilon }_{\text{r}} $的值. 具体而言，本研究拟合的是工况$ {\lambda }_{1\text{p}}=0.9 $，$ {\lambda }_{\text{2p}}=6 $下的力电响应曲线，得到八链模型和网络平均管模型的介电常数$ {\varepsilon }_{\text{r}} $的值分别为4.4、3.7. 采用不同的模型对不同预拉伸工况下的力电耦合效应进行预测，结果如图7所示. 其中，$ {\mathit{\Phi }}^{\exp } $为实验记录的驱动电压，$ {\mathit{\Phi }}^{8\text{-chain}} $和$ {\mathit{\Phi }}^{\text{NAT}} $分别为八链模型和网络平均管模型预测的理论驱动电压. 随着驱动电压升高，麦克斯韦应力诱导薄膜厚度减小，导致侧向拉伸比$ {\lambda }_{1} $呈现非线性增长趋势.

采用决定系数$ {R}^{2} $来衡量拟合和预测的精度：

(36)$ {R}^{2}=1-{\dfrac{\displaystyle\sum\limits_{i=1}^{n}({\mathit{\Phi }}_{\mathrm{Exp}}-{\mathit{\Phi }}_{\mathrm{Sim}})^{2}}{\displaystyle\sum\limits_{i=1}^{n}{({{\mathit{\Phi }}_{\mathrm{Exp}}}-{{\overline{\mathit{\Phi }}}_{\mathrm{Exp}}})}^{2}}} . $

式中： $ {\mathit{\Phi }}_{\mathrm{Exp}} $为实验测量驱动电压，${{\overline{\mathit{\Phi }}}_{\mathrm{Exp}}} $为实验测量平均驱动电压，$ {\mathit{\Phi }}_{\mathrm{Sim}} $为预测驱动电压. 对应的模型的决定系数如表2所示. 考虑到标准$ {R}^{2} $指标在曲线的高斜率区段对纵向绝对误差过度敏感，因此，计算拟合优度时，$ {\lambda }_{1\text{p}}=0.5 $的$ {R}^{2} $不计算在内. 在熵弹性模型的基础上，引入缠结约束自由能项可以更准确地表征文献[12]的实验数据. 相比于文献报道的EV模型和Gent模型以及八链模型，考虑缠结效应的本构模型在预测介电弹性体的力电耦合效应时表现出更好的趋势一致性与更小的误差.

研究结果表明，与Gent模型、EV模型相比，本研究建立的考虑缠结效应的弹性体本构模型，利用单轴拉伸实验拟合得到的参数，能够较好地预测横向约束模式下的力电耦合响应. 本研究中，虽然缺乏聚丙烯酸酯弹性体在等双轴拉伸模式下的力电耦合试验数据，但可基于此理论框架，将纯剪切试验预测良好的本构模型合理外推至等双轴拉伸. 具体而言，在等双轴拉伸试验中，变形梯度张量$ \boldsymbol{F} $表达式为

(37)$ \boldsymbol{F}={\lambda }_{\text{p}}{\boldsymbol{e}}_{1}\otimes {\boldsymbol{e}}_{1}+{\lambda }_{\text{p}}{\boldsymbol{e}}_{2}\otimes {\boldsymbol{e}}_{2}+\frac{1}{\lambda _{\text{p}}^{2}}{\boldsymbol{e}}_{3}\otimes {\boldsymbol{e}}_{3} . $

将式(37)代入到式(34)，并选取预拉伸比分别为$ {\lambda }_{\text{p}}=1,2,3,4,5 $，得到：

(38)$ {P}_{1\text{p}}={\frac{\partial {W}_{{\mathrm{s}}}}{\partial {\lambda }_{\text{1p}}}}_{{{\lambda }_{2}}={{\lambda }_{\text{p}}}}-\varepsilon {\left(\frac{\mathit{\Phi }}{{L}_{3}}\right)}^{2}\lambda _{\text{p}}^{3} . $

在施加恒定载荷$ {P}_{1{\mathrm{p}}} $并移除位移约束后，材料在电场作用下产生进一步的等双轴变形. 此时，定义当前的等双轴拉伸比为$ \lambda $，则此时的变形梯度张量应改写为一般形式；

(39)$ \boldsymbol{F}=\lambda {\boldsymbol{e}}_{1}\otimes {\boldsymbol{e}}_{1}+\lambda {\boldsymbol{e}}_{2}\otimes {\boldsymbol{e}}_{2}+\frac{1}{{\lambda }^{2}}{\boldsymbol{e}}_{3}\otimes {\boldsymbol{e}}_{3} , $

(40)$ \mathit{\Phi }={L}_{3}\times \sqrt{\left({\dfrac{\partial {W}_{\text{s}}}{\partial \lambda }-{P}_{1\text{p}}}\right)\bigg/\left({\varepsilon {\lambda }^{3}}\right)} . $

其力电响应的预测结果如图8所示. 可以看出，预拉伸可以有效地提高介电弹性体的驱动性能，且预拉伸比越大，实现相同面内拉伸比所需施加的驱动电压（或电场强度）显著降低. 这一预测结果为通过机械预加载优化介电弹性体驱动器的能效提供了理论依据.

本研究建立的模型基于超弹性与理想介电的平衡态假设，未考虑VHB 4905材料固有的黏弹性行为. 尽管实验主要通过低应变率加载及应力弛豫过程来获取准静态条件的数据^[12]，但材料在实际加载中仍可能存在黏弹性效应. 该简化有助于理解缠结效应对于力电耦合响应的贡献. 因此，本模型主要适用于预测稳态或平衡态下的力电响应，未针对不同应变速率下的材料行为开展研究. 未来工作将在该模型基础上引入黏弹性效应，特别是考虑缠结效应的黏弹性机制.

本研究建立了考虑缠结约束效应的介电弹性体力-电耦合本构模型，该理论的核心在于通过非高斯分布的Langevin统计模型来描述分子链构象熵对弹性自由能的贡献，并基于Doi-Edwards提出的管模型来描述由缠结约束效应引起的自由能贡献. 利用网络平均方法建立微观链段变形与宏观应变量的定量关系，建立三参数的超弹性本构模型，在此基础上引入理想介电模型，最终建立四参数的力-电耦合模型. 与采用经典超弹性模型的力-电耦合模型相比，本研究所建立的考虑缠结效应的介电弹性体本构模型能够更好地表征文献中的介电弹性体中的力学响应和力电耦合效应，即考虑了缠结效应的介电弹性体本构模型在同时预测单轴拉伸、双轴拉伸、力电耦合等情况下均表现出较为优异的预测结果.

本研究所建立的力电耦合模型，仅包含4个参数，但能同时较好地描述介电弹性体的力学变形行为和电致变形行为.后续可在该模型的基础上进一步引入黏弹性效应，用以预测介电弹性体材料在动态加载下的响应行为.