弹性有限变形

发表于 2024-09-02 更新于 2025-03-25 分类于本构理论阅读次数：

弹性有限变形的一般本构方程

在本节推导中假设如下:

材料在参考构型中的内能体积密度 \(W = \rho_0 e\) 可以表示为变形梯度的函数, \(W= W(F)\);
不考虑材料与外部的热交换, 因此能量守恒定律表述为

\[ \begin{equation} \frac{\mathrm{D} W}{\mathrm{D} t} = J \sigma : v \nabla \label{eq:energy_conservation} \end{equation} \]

式中, \(J \triangleq \det F = \rho_0/\rho\). 上式左端应用链式法则得到 \[ \begin{equation} \frac{\mathrm{D} W}{\mathrm{D} t} = \frac{\partial W}{\partial F} : \frac{\mathrm{D} F}{\mathrm{D} t} = \left( \frac{\partial W}{\partial F} \cdot F^{\top} \right) : v \nabla \label{eq:matl_deriv_of_e_density} \end{equation} \] 对比式 \(\ref{eq:energy_conservation}\) 与式 \(\ref{eq:matl_deriv_of_e_density}\), 再由速度梯度场的任意性可得等式 \[ \sigma = \frac{1}{J} \frac{\partial W}{\partial F} \cdot F^{\top} \] 上式从形式看很简单, 但需要确定 \(W\) 关于变形梯度 \(F\) 9 个分量的函数形式. 因此在实际操作中不可用. 接下来通过本构的客观性原理和材料的对称性减少 \(W\) 的函数独立变量数量.

客观性原理对内能体积密度函数的简化

考虑变形 1 和变形 2, 变形 1 给出变形梯度场 \(F\), 变形 2 在变形 1 的基础上作施加旋转场 \(Q\), 因此变形 2 的变形梯度场为 \[ F' = Q \cdot F \] 经旋转之后, 内能密度在每个物质点 \(X\) 处应保持不变, 因此有 \[ W(F) \equiv W(Q \cdot F) \] 又因为变形梯度可以通过极分解原理分解为正交张量 \(R\) 和对称正定张量 \(U\), 作变量替换, 得到关于 \(U\) 的内能密度函数 \(\widetilde W\) \[ \widetilde W(U, R) \triangleq W(R \cdot U) \] 在任意旋转下 \(\widetilde W\) 保持不变, 所以 \[ \widetilde W(U, R) \equiv W(Q \cdot R \cdot U) \] 因此总可以选取一个旋转场: \(Q \equiv R^{-1}\), 上式转化为 \[ \widetilde W(U, R) \equiv W( U) \] 因此, 内能密度函数可以表示为关于右伸长张量 \(U\) 的函数, 并且函数形式与 \(W(U)\) 相同. 这样就将分量的个数简化为 6 个. 又因为 \(U = \sqrt{C}\), 所以 \[ \hat{W}(C) = W(\sqrt{C}) \] 在上述推导中, 为了区分不同的函数表达形式, 使用记号 \(\widetilde{W}\) 或 \(\hat{W}\), 之后则不加区别地使用函数记号 \(W\).

如果将内能密度考虑为右 Cauchy-Green 张量 \(C\) 的函数, 那么 \(W\) 的物质导数为 \[ \frac{\mathrm{D} W}{\mathrm{D} t} = \left[ F \cdot \left( \frac{\partial W}{\partial C} + \frac{\partial W}{\partial C^{\top}} \right) \cdot F^{\top} \right] : v \nabla \] 如果 \(W\) 可以表示为关于分量 \(C_{IJ}\) 的对称函数, 也即 \(W(C_{IJ}) \equiv W(C_{JI})\), 那么上式又可以化简为 \[ \frac{\mathrm{D} W}{\mathrm{D} t} = \left( F \cdot 2 \frac{\partial W}{\partial C} \cdot F^{\top} \right) : v \nabla \] 再比较上式和式 \(\ref{eq:energy_conservation}\), 就得到 \[ \begin{equation} \boxed{ \sigma = \frac{1}{J} F \cdot 2 \frac{\partial W}{\partial C} \cdot F^{\top} } \label{eq:sigma_rep_as_psi_to_C} \end{equation} \]

[!NOTE]

内能密度函数与 PK2 应力的关系

Cauchy 应力与 PK2 应力之间的关系为 \[ F \cdot S \cdot F^{\top} = J \sigma \quad \] 再根据式 \(\ref{eq:sigma_rep_as_psi_to_C}\), 就得到 \[ \boxed{ S = 2 \frac{\partial W}{\partial C} } \]

[!NOTE]

小变形下的近似

在小变形下, \[ C \sim I + 2 \varepsilon \] 因此 \[ \frac{\partial W}{\partial \varepsilon} \sim 2 \frac{\partial W}{\partial C} \] 将式 (4) 展开为位移梯度的函数 \[ \sigma \sim \frac{1}{1 + \mathrm{tr}(u\nabla)} (I + u \nabla) \cdot \frac{\partial W}{\partial \varepsilon} \cdot (I + \nabla u) \\ \sim \frac{\partial W}{\partial \varepsilon} + 2\frac{\partial W}{\partial \varepsilon} \cdot \varepsilon + \mathcal{O}(\varepsilon^2) \] 忽略掉关于 \(\varepsilon\) 的一阶和二阶项, 得到 \[ \boxed{ \sigma = \frac{\partial W}{\partial \varepsilon} } \]

各向同性对内能密度函数的简化

考虑旋转矩阵 \(Q\), 若材料对 \(Q\) 具有旋转不变性, 那么 \[ W(C) \equiv W(Q^{\top} \cdot C \cdot Q) \] 特别的, 对于各向同性材料, \(Q\) 是任意的, 上式等价于 \(W\) 在任意坐标变化下形式不变, 因此内能密度函数 \(W\) 可以写成应变不变量 \(\ref{eq:invariants_1}\) 的函数: \[ W = W(\mathscr{I}_1, \mathscr{I}_2, \mathscr{I}_3) \] 根据链式法则, \(W\) 对右 Cauchy-Green 张量 \(C\) 的导数为 \[ \frac{\partial W}{\partial \boldsymbol{C}} = \frac{\partial W}{\partial \mathscr{I}_1} \frac{\partial \mathscr{I}_1}{\partial \boldsymbol{C}} + \frac{\partial W}{\partial \mathscr{I}_2} \frac{\partial \mathscr{I}_2}{\partial \boldsymbol{C}} + \frac{\partial W}{\partial \mathscr{I}_3} \frac{\partial \mathscr{I}_3}{\partial \boldsymbol{C}} \] 式中, 应变不变量对 \(\boldsymbol{C}\) 的导数为 \[ \boxed{ \frac{\partial \mathscr{I}_1}{\partial \boldsymbol{C}} = \boldsymbol{I}, \quad \frac{\partial \mathscr{I}_2}{\partial \boldsymbol{C}} = \mathscr{I}_1\boldsymbol{I} - \boldsymbol{C}, \quad \frac{\partial \mathscr{I}_3}{\partial \boldsymbol{C}} = \mathscr{I}_3\boldsymbol{C}^{-1} } \] 代入 \(\partial W / \partial \boldsymbol{C}\), \[ \frac{\partial W}{\partial \boldsymbol{C}} = W_{,1} \boldsymbol{I} + W_{,2} (\mathscr{I}_1\boldsymbol{I} - \boldsymbol{C}) + W_{,3}\mathscr{I}_3 \boldsymbol{C}^{-1} \] 式中, \(W_{,i} = {\partial W}/{\partial \mathscr{I}_i}\). 再根据式 \(\ref{eq:sigma_rep_as_psi_to_C}\) , 并应用 \(\mathscr{I}_3 = J^2\) , 就得到应力的表达式 \[ \boldsymbol{\sigma} = \frac{1}{\sqrt{\mathscr{I}_3}} \boldsymbol{F} \cdot 2 \left( W_{,1} \boldsymbol{I} + W_{,2} (\mathscr{I}_1\boldsymbol{I} - \boldsymbol{C}) + W_{,3}\mathscr{I}_3 \boldsymbol{C}^{-1} \right) \cdot \boldsymbol{F}^{\top} \] 注意到 \(\boldsymbol{B}\) 与 \(\boldsymbol{C}\) 之间的关系满足 \[ \boldsymbol{F} \cdot \boldsymbol{C}^{n} \cdot \boldsymbol{F}^{\top} = \boldsymbol{F} \cdot ( \boldsymbol{F}^{\top} \cdot \boldsymbol{F} \cdot \boldsymbol{F}^{\top} \cdots \boldsymbol{F}^{\top} \cdot \boldsymbol{F} ) \cdot \boldsymbol{F}^{\top} = \boldsymbol{B}^{n+1}, \quad n \in \mathbb{Z} \] 所以应力 \(\sigma\) 又可以写作关于 \(B\) 的形式 \[ \boxed{ \boldsymbol{\sigma} = \frac{1}{\sqrt{\mathscr{I}_3}} 2 \left( W_{,3}\mathscr{I}_3 I +(W_{,1} + W_{,2} \mathscr{I}_1) B - W_{,2} B^2 \right) } \] 根据 Cayley-Hamilton 定理, 得到关于 \(B\) 的代数恒等式 \[ \boldsymbol{B}^2 \equiv \mathscr{I}_1 \boldsymbol{B} - \mathscr{I}_2 \boldsymbol{I} + \mathscr{I}_3 \boldsymbol{B}^{-1} \] 代回之前用 \(\boldsymbol{B}\) 表示的应力表达式中, 得到使用 \(\{ \boldsymbol{B}^{-1}, \boldsymbol{B} \}\) 表示的应力 \[ \boxed{ \boldsymbol{\sigma} = \frac{1}{\sqrt{\mathscr{I}_3}} 2 \left( W_{,1} \boldsymbol{B} + (W_{,2} \mathscr{I}_2 + W_{,3}\mathscr{I}_3) \boldsymbol{I} - W_{,2} \mathscr{I}_3 \boldsymbol{B}^{-1} \right) } \]

[!TIP]

使用 \(\{\bar{\mathscr{I}}_1, \bar{\mathscr{I}}_2, J \}\) 表示各向同性弹性材料

通过链式法则, 可以得到 \[ \begin{equation} \boxed{ \boldsymbol{\sigma} = \frac{2}{J} \left[ \frac{1}{J^{2/3}} \left( W_{,1} + \bar{\mathscr{I}}_1 W_{,2} \right) \boldsymbol{B} - \frac{1}{3} (\bar{\mathscr{I}}_1 W_{,1} + 2\bar{\mathscr{I}}_2 W_{,2}) \boldsymbol{I} - \frac{1}{J^{4/3}} W_{,2} \boldsymbol{B}^2 \right] + W_{,3} \boldsymbol{I} } \label{eq:cons_rep_as_B_inva_3} \end{equation} \]

考虑材料不可压缩的情况

当材料不可压缩时 \[ \mathscr{I}_3 \equiv 1 \] 对于不可压缩情况, 不可以直接令 \(\mathscr{I}_3 = 1\), 因为这会导致式中函数 \(W\) 导数项爆破. 使用 Lagrange 乘子将不可压缩的特性考虑进去: \[ \widetilde W = W(\mathscr{I}_1, \mathscr{I}_2) - \frac{1}{2}p (\mathscr{I}_3 -1) \]

\[ \sigma = - pI + (2W_{,1} + 2W_{,2}\mathscr{I}_1) B - 2W_{,2} B^2 \]

动力学约束中, \(p\) 只能通过平衡方程求解, 无法通过本构关系得到

一些超弹本构

Hypoelasticity 材料

在绝大多数情况下, 只有当材料受剪时才会表现出非线性的响应 (与 \(\mathscr{I}'_2\) 相关), 同时应力与体积变化呈线性关系 (与 \(\mathscr{I}_1\) 相关). 一种定义内能密度函数的方式为 \[ W(\mathscr{I}_1, \mathscr{I}'_2) = \frac{1}{6} K \mathscr{I}_1^2 + \frac{2n \sigma_0 \varepsilon_0}{n+1} \left( \frac{\mathscr{I}'_2}{\varepsilon_0^2} \right)^{(n+1)/2n} \] 对应变 \(\varepsilon\) 求导数得到 \[ \boldsymbol{\sigma} = \frac{ \partial W }{ \partial \varepsilon} = \frac{1}{3}K \mathscr{I}_1 \mathbf{I} + \sigma_0 \left( \frac{\mathscr{I}'_2}{\varepsilon_0^2} \right)^{(1-n)/2n} \frac{\boldsymbol{\varepsilon}'}{\varepsilon_0} \] 或者表示成关于应力的表达式: \[ \boldsymbol{\varepsilon} = \frac{1}{3K} \mathscr{J}_1 \mathbf{I} + \varepsilon_0 \left( \frac{\mathscr{J}_2}{\sigma_0^2} \right)^{(n-1)/2} \frac{\boldsymbol{\sigma}'}{\sigma_0} \] 式中, \(\mathscr{J}_1\), \(\mathscr{J}_2\) 为应力的第一和第二不变量: \[ \mathscr{J}_1 \triangleq \mathrm{tr}~\sigma, \quad \mathscr{J}_2 \triangleq \frac{1}{2} \left( \mathrm{tr}~ \sigma^2 - \frac{1}{3}(\mathrm{tr}~\sigma)^2 \right) \] 由此可以得到以下三种工况的应力应变关系:

单轴拉伸. 应力分量为 \(\sigma_1\). 此时没有剪切变形, 应变分量为 \[ \varepsilon_{11} = \frac{1}{3K} \sigma_1 + \frac{2 \varepsilon_0}{\sqrt{3}} \left( \frac{\sigma_1}{\sqrt{3} \sigma} \right)^{n} \]
承受静水压力 \(-p\). 此时 \(\mathscr{J}_2 = 0\), 因此变形为 \[ \boldsymbol{\varepsilon} = -\frac{p}{K} \mathbf{I} \]
承受剪切力 \(\tau\). 此时应变只有对应剪应力分量处非零, 并满足幂形式: \[ \varepsilon = \varepsilon_0 \left( \frac{\tau}{\sigma_0} \right)^n \]

Neo-Hookean 材料

内能密度函数为 \[ \boxed{ W(\bar{\mathscr{I}}_1, J) = \frac{1}{2} \mu (\bar{\mathscr{I}}_1 - 3) + \frac{1}{2} K (J - 1)^2 } \]

式中, \(\mu\) 和 \(K\) 为材料参数, 在小变形的情况下分别对应剪切和 bulk 模量. 对于近不可压缩材料, \(K \gg \mu\). 将 \(W\) 对不变量的导数带入到式 \(\ref{eq:cons_rep_as_B_inva_3}\) 中得到 \[ \boxed{ \boldsymbol{\sigma} = \frac{\mu}{J^{5/3}} ~\mathrm{dev}\boldsymbol{B} + K(J-1)\boldsymbol{I} } \]

Mooney-Rivlin 材料

内能密度函数为 \[ \boxed{ W(\bar{\mathscr{I}}_1, \bar{\mathscr{I}}_2, J) = \frac{1}{2} \mu_1 (\bar{\mathscr{I}}_1 - 3) + \frac{1}{2} \mu_2 (\bar{\mathscr{I}}_2 - 3) + \frac{1}{2} K (J - 1)^2 } \] 式中, \(\mu_1\), \(\mu_2\) 和 \(K\) 为材料参数, 对于近不可压缩材料, \(K \gg \mu_1 , \mu_2\). 应力应变关系为 \[ \boxed{ \boldsymbol{\sigma} = \frac{\mu_1}{J^{5/3}} ~\mathrm{dev}\boldsymbol{B} + \frac{\mu_2}{J^{7/3}} \left( \mathscr{I}_1 \boldsymbol{B} - \frac{2}{3} \mathscr{I}_2 \boldsymbol{I} - \boldsymbol{B}^2 \right) + K(J-1)\boldsymbol{I} } \]

附录: 应变的不变量

我们使用右 Cauchy-Green 张量 \(C\) 的不变量定义应变的不变量 \[ \begin{equation} \begin{gathered} \mathscr{I}_1 = \mathrm{tr}~C \\ \mathscr{I}_2 = \frac{1}{2} \left( (\mathrm{tr}~C)^2 - \mathrm{tr}~C^2 \right) \\ \mathscr{I}_3 = \det C \label{eq:invariants_1} \end{gathered} \end{equation} \] 如果定义 \(C\) 的特征值为 \(\lambda_1^2\), \(\lambda_2^2\), \(\lambda_3^2\), 那么上述应变不变量又可以表述为 \[ \mathscr{I}_1 = \lambda_1^2 + \lambda_2^2 + \lambda_3^2 \\ \mathscr{I}_2 = \lambda_1^2 \lambda_2^2 + \lambda_1^2 \lambda_3^2 + \lambda_1^2 \lambda_3^2 \\ \mathscr{I}_3 = \lambda_1^2 \lambda_2^2\lambda_3^2 \]

注意到 \(C\) 的不变量与左 Cauchy-Green 张量 \(B\) 相同, 因此式 (5) 中将 \(C\) 替换为 \(B\) 同样成立.

Cayley-Hamilton 定理给出关于 \(C\) \(B\) 以及不变量作为系数的恒等式: \[ C^3 - \mathscr{I}_1 C^2 + \mathscr{I}_2 C - \mathscr{I}_3 I = 0 \\ B^3 - \mathscr{I}_1 B^2 + \mathscr{I}_2 B - \mathscr{I}_3 I = 0 \]

[!NOTE]

如果考虑仿射量 \(A\) 的 deviatoric 部分: \[ A' = A - \frac{1}{3} (\mathrm{tr}A)~ I \] \(A'\) 的不变量定义与式 (0) 一致, 那么 \[ \begin{equation} \begin{gathered} \mathscr{I}'_1 = \mathrm{tr}A' = 0 \\ \mathscr{I}'_2 = \frac{1}{2} \left( (\mathrm{tr}A')^2 - \mathrm{tr}A'^2 \right) \\ \mathscr{I}'_3 = \det A' \end{gathered} \label{eq:invariants_2} \end{equation} \] 不变量 \((\mathscr{I}'_1, \mathscr{I}'_2, \mathscr{I}'_3)\) 也可以通过 \((\mathscr{I}_1, \mathscr{I}_2, \mathscr{I}_3)\) 表示, 或者 \((\mathrm{tr}A, \mathrm{tr}A^2, \mathrm{tr}A^3)\) 表示. 例如, \(A'\) 的第二不变量可以表示为 \[ \mathscr{I}'_2 = \frac{1}{2} \left( (\mathrm{tr}A')^2 - \mathrm{tr}A'^2 \right) = -\frac{1}{2} \mathrm{tr}\left( A^2 - \frac{2}{3} (\mathrm{tr}A)~ A + \frac{1}{9} (\mathrm{tr}A)^2~ I \right) = -\frac{1}{2} \mathrm{tr}A^2 + \frac{1}{6} (\mathrm{tr}A)^2 = \mathscr{I}_2 -\frac{1}{3} \mathscr{I}_1^2 \] 可以证明, 上述定义的 \(\mathscr{I}'_2\) 总是负值, 如果 \(A\) 的特征值总是大于 0. 因此我们重定义 \(\mathscr{I}'_2\) 为原来的负数, 使得它总是大于 0: \[ \mathscr{I}'_2 := -\mathscr{I}'_2 > 0 \]

[!NOTE]

应用于近不可压缩情况下的应变不变量

\[ \begin{equation} \begin{aligned} \bar{\mathscr{I}}_1 &= \frac{\mathscr{I}_1}{J^{2/3}}, \\ \bar{\mathscr{I}}_2 &= \frac{\mathscr{I}_2}{J^{4/3}}, \\ J &= \sqrt{\mathscr{I}_3} \end{aligned} \label{eq:invariants_3} \end{equation} \]