跳转至

内力虚功的离散化

增量分析框架中给出的时刻 \(t + \Delta t\) 的虚功方程，根据参考构形的选择分为 Updated Lagrange 法和 Total Lagrange 法两种形式。本章使用形函数与有限元近似和形函数的空间导数中引入的有限元近似，对两种公式中的内力虚功进行空间离散化，得到单元内力向量 \(\boldsymbol{q}^e\)（UL 法）和 \(\boldsymbol{Q}^e\)（TL 法）。

对于单元 \(e\)，设构成节点为 \(\alpha = 1, \ldots, n_e\)，其位移为 \(\boldsymbol{u}^e_\alpha\)，并将单元节点位移向量排列为 \(\boldsymbol{u}^e = (\boldsymbol{u}^{eT}_1, \ldots, \boldsymbol{u}^{eT}_{n_e})^T\)；虚位移 \(\delta \boldsymbol{u}^e\) 也按相同顺序定义。单元内位移由形函数插值为 \(\boldsymbol{u} = \sum_\alpha N_\alpha^e \boldsymbol{u}^e_\alpha\)。

Updated Lagrange 法的内力虚功

在 Updated Lagrange 法中，以时刻 \(t\) 的当前构形 \({}^{t}\Omega\) 为参考构形，并用 Cauchy 应力 \(\boldsymbol{\sigma}\) 和 Almansi 应变线性部分 \(\boldsymbol{A}_{(L)}\) 将内力虚功写为

\[ \delta W^{\mathrm{int}} = \sum_e \int_{\Omega^e} \boldsymbol{\sigma} : \delta \boldsymbol{A}_{(L)}\, dv \]

。\(\delta \boldsymbol{A}_{(L)}\) 的各分量可以表示为关于当前构形坐标 \(\boldsymbol{x}\) 的形函数导数 \(\partial N_\alpha^e/\partial x_i\) 与节点虚位移 \(\delta u^e_{i\alpha}\) 的线性组合。在 Voigt 表示中可整理为

\[ \delta \boldsymbol{A}_{(L)} = \boldsymbol{B}_L\, \delta \boldsymbol{u}^e, \qquad \boldsymbol{B}_L = [\boldsymbol{B}_{L1}, \ldots, \boldsymbol{B}_{Ln_e}] \]

。节点块 \(\boldsymbol{B}_{L\alpha}\) 是按照 Voigt 规则排列 \(\partial N_\alpha^e/\partial x_i\) 得到的 \(6 \times 3\) 矩阵，\(\boldsymbol{B}_L\) 即 UL 法的应变-位移矩阵。将其代入内力虚功并分离 \(\delta \boldsymbol{u}^e\)，得到

\[ \delta W^{\mathrm{int}} = \sum_e \delta \boldsymbol{u}^{eT}\, \boldsymbol{q}^e, \qquad \boldsymbol{q}^e = \int_{\Omega^e} \boldsymbol{B}_L^T \boldsymbol{\sigma}\, dv \]

。\(\boldsymbol{q}^e\) 的节点块 \(\boldsymbol{q}^e_\alpha = \int_{\Omega^e} \boldsymbol{B}_{L\alpha}^T \boldsymbol{\sigma}\, dv\) 是单元 \(\Omega^e\) 作用于构成节点 \(\alpha\) 的内力。

Total Lagrange 法的内力虚功

在 Total Lagrange 法中，以参考构形 \(\Omega_0\) 为参考，并用第 2 Piola-Kirchhoff 应力 \(\boldsymbol{S}\) 和 Green-Lagrange 应变 \(\boldsymbol{E}\) 写为

\[ \delta W^{\mathrm{int}} = \sum_e \int_{\Omega^e_0} \boldsymbol{S} : \delta \boldsymbol{E}\, dV \]

。\(\delta \boldsymbol{E}\) 可分为关于虚位移的线性项，以及包含与当前位移梯度 \(\partial u_k/\partial X_j\) 乘积的非线性项：

\[ \delta E_{(L)ij} = \frac{1}{2}\left( \frac{\partial \delta u_i}{\partial X_j} + \frac{\partial \delta u_j}{\partial X_i} \right), \quad \delta E_{(NL)ij} = \frac{1}{2}\left( \frac{\partial \delta u_k}{\partial X_i}\, \frac{\partial u_k}{\partial X_j} + \frac{\partial u_k}{\partial X_i}\, \frac{\partial \delta u_k}{\partial X_j} \right). \]

线性项通过对 \(\partial N_\alpha^e/\partial X_i\) 应用与 UL 法相同的排列规则所得的节点块 \(\boldsymbol{B}_{L\alpha}\) 写为

\[ \delta \boldsymbol{E}_{(L)} = \boldsymbol{B}_L\, \delta \boldsymbol{u}^e \]

（由于参考构形不同，构成项仅由 \(\partial N_\alpha^e/\partial x_i\) 替换为 \(\partial N_\alpha^e/\partial X_i\)，符号与 UL 法共用）。非线性项通过按照 Voigt 规则排列当前位移梯度 \(\partial u_k/\partial X_j\) 与 \(\partial N_\alpha^e/\partial X_i\) 的乘积所得的节点块 \(\boldsymbol{B}_{NL\alpha}\) 写为

\[ \delta \boldsymbol{E}_{(NL)} = \boldsymbol{B}_{NL}\, \delta \boldsymbol{u}^e, \qquad \boldsymbol{B}_{NL} = [\boldsymbol{B}_{NL1}, \ldots, \boldsymbol{B}_{NLn_e}] \]

。因此 \(\delta \boldsymbol{E} = (\boldsymbol{B}_L + \boldsymbol{B}_{NL})\, \delta \boldsymbol{u}^e\)，而 \(\boldsymbol{B}_L + \boldsymbol{B}_{NL}\) 即 TL 法的应变-位移矩阵。将其代入内力虚功，得到

\[ \delta W^{\mathrm{int}} = \sum_e \delta \boldsymbol{u}^{eT}\, \boldsymbol{Q}^e, \qquad \boldsymbol{Q}^e = \int_{\Omega^e_0} (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T \boldsymbol{S}\, dV \]

。节点块 \(\boldsymbol{Q}^e_\alpha\) 是单元 \(\Omega^e_0\) 作用于构成节点 \(\alpha\) 的内力。

UL/TL 的对应关系与计算流程

Updated Lagrange 法与 Total Lagrange 法的单元内力向量对应关系如下。

项目 Updated Lagrange 法 Total Lagrange 法
参考构形 当前构形 \({}^{t}\Omega^e\) 参考构形 \(\Omega^e_0\)
应力张量 Cauchy 应力 \(\boldsymbol{\sigma}\) 第 2 PK 应力 \(\boldsymbol{S}\)
应变变分 \(\delta \boldsymbol{A}_{(L)}\) \(\delta \boldsymbol{E} = \delta \boldsymbol{E}_{(L)} + \delta \boldsymbol{E}_{(NL)}\)
B 矩阵 \(\boldsymbol{B}_L\)（\(\partial N/\partial x\)） \(\boldsymbol{B}_L + \boldsymbol{B}_{NL}\)（\(\partial N/\partial X\)，\(\partial u/\partial X\)）
单元内力 \(\boldsymbol{q}^e = \int_{\Omega^e} \boldsymbol{B}_L^T \boldsymbol{\sigma}\, dv\) \(\boldsymbol{Q}^e = \int_{\Omega^e_0} (\boldsymbol{B}_L + \boldsymbol{B}_{NL})^T \boldsymbol{S}\, dV\)

两者可采用相同的处理流程：“由形函数的空间导数组成 \(\boldsymbol{B}_L\)”；“TL 法中由当前位移梯度组成 \(\boldsymbol{B}_{NL}\) 并相加”；“按照本构关系更新应力（\(\boldsymbol{\sigma}\) 或 \(\boldsymbol{S}\)）”；“在积分点上对 \(\boldsymbol{B}^T \boldsymbol{\sigma}\) 或 \((\boldsymbol{B}_L+\boldsymbol{B}_{NL})^T \boldsymbol{S}\) 在单元区域内进行数值积分（数值积分）”。除切换参考构形（节点坐标及 \(\boldsymbol{B}\) 矩阵的构成）和替换应力张量之外，其余处理均相同，因此 FrontISTR 中两种方法的内力计算采用共用子程序实现。从单元内力向量 \(\boldsymbol{q}^e\)、\(\boldsymbol{Q}^e\) 向全局内力向量的组装见外力虚功与全局方程的组装。

相关项目