完全拉格朗日形式的矩阵列式
完全拉格朗日形式的弱形式可以解释为虚功原理, 也即内力虚功等于惯性力和外力虚功: \[ \delta w^{\mathrm{int}} = \delta w^{\mathrm{kin}} + \delta w^{\mathrm{ext}} \tag{1} \] 使用形函数将待求解的变量 \(\varphi(X,t)\) 的空间维度进行离散得到 \[ \varphi_i^{h}(X,t) = N^a(X) \varphi_{i}^{a}(t) \] 式中, 指标 \(i\) 表示拉格朗日描述下的空间指标, 指标 \(a\) 表示节点编号指标, 并默会 Einstein 求和约定. 在完全拉格朗日形式下, 节点内力, 节点外力, 和单元质量矩阵分别在空间维度上离散为 \[ M_{ij}^{ab} \triangleq \delta_{ij}\int_{B_0} \rho_0 N^{a} N^{b} ~\mathrm{d} X, \\ f_{i}^{a,\mathrm{int}} \triangleq \int_{B_0} P_{iJ} \frac{ \partial N^a}{ \partial X_J} ~\mathrm{d} X, \\ f_{i}^{a,\mathrm{ext}} \triangleq \int_{B_0} \rho_0 b_{i} N^a ~\mathrm{d} X + \int_{\partial B_0^t} \overline{T}_i N^a ~\mathrm{d} \Sigma(X) \tag{2} \] 式中, \(\rho_0(X)\) 为参考构型的材料密度; \(P_{IJ}\) 为 PK1 应力, \(P \triangleq J \sigma \cdot F^{-\top}\); \(b_i(X)\) 为体积力; \(\overline{T}_{i}(X)\) 为力边界条件上指定的面力; \(N^a\) 为形函数. 上式中出现的积分是对物质坐标 \(X\) 的积分, 因此积分域为参考构型或参考构型的边界. 在式 (2) 的定义下, 虚功原理 (1) 可以表示为 \[ \delta \varphi_{i}^{a} (M_{ij}^{ab} \ddot{\varphi}_{j}^{b} + f_{i}^{a,\mathrm{int}} - f_{i}^{a,\mathrm{ext}}) = 0 \] 对任意可允许的节点虚位移 \(\delta \varphi_{i}^{a}\) 成立, 所以得到在 \((a,i)\) 自由的方程组 \[ M_{ij}^{ab} \ddot{\varphi}_{j}^{b} + f_{i}^{a,\mathrm{int}} - f_{i}^{a,\mathrm{ext}} = 0 \]
式 (1) 中的虚功项可以拆解成各个单元上的积分, 以内力虚功为例: \[ \delta w^{\mathrm{int}} = \delta \varphi_{i}^{a} \sum_{e}\int_{B_0^e} P_{iJ} \frac{ \partial N^a}{ \partial X_J} ~\mathrm{d} X \tag{3} \] 如果定义单元内力为 \[ (f_{i}^{a,\mathrm{int}})^{e} \triangleq \int_{B_0^e} P_{iJ} \frac{ \partial N^a}{ \partial X_J} ~\mathrm{d} X \] 那么式 (3) 可以重新写作 \[ \delta w^{\mathrm{int}} = \delta \varphi_{i}^{a} \sum_{e}(f_{i}^{a,\mathrm{int}})^{e} \] 由于有限元插值函数的局部性, 单元内力只会和所属于该单元的节点处做功, 在其它节点处的虚位移不会对虚功做出贡献. 因此, 内力虚功可以分解为单元内力虚功之和: \[ \delta w^{\mathrm{int}} = \sum_{e}\delta \varphi_{i}^{a,e} (f_{i}^{a,\mathrm{int}})^{e} = \sum_{e} \delta w^{\mathrm{int},e} \] 上式说明, 单元上的虚功列式在形式上与系统的虚功列式完全相同. 因此总是可以在单元上计算得到单元的内力,外力, 质量矩阵, 再通过 “gather” - “scatter” 数据操作组装得到系统的内力, 外力和质量矩阵. 在之后的叙述中, 我们总是默认式 (2) 的运算在单元上进行, 不再加上标 \(e\).
定义的标量值
| 变量名 | 变量说明 |
|---|---|
nedim |
参数坐标维数 |
nsdim |
空间维数 |
nndof |
节点自由度数 |
nen |
单元节点数 |
nedof |
单元自由度数, nedof = nndof * nen |
在实现中, \(\varphi_I^a\)
经常被排列成列向量的形式: \[
[\varphi] \triangleq [\varphi^{1}; \varphi^{2}; \ldots;
\varphi^{\mathrm{nen}}]
\] 上式中的列向量长度为 nedof*nn. 引入自由度指标
\(\alpha: (a,i) \mapsto \alpha\),
并定义形函数矩阵 \([\mathrm{N}]_{I\overline{a}}\):
nndof x nndof \[
[\mathrm{N}]_{I \alpha} = [N^{1} \mathrm{I}, N^{2} \mathrm{I}, \ldots,
N^{\mathrm{nn}} \mathrm{I}]
\] 式中, \(\mathrm{I}\) 是维度为
nndof 的单位矩阵, 因此, 被插值的变量可以写成如下矩阵形式
\[
[\varphi^{h}] = [\mathrm{N}] [\varphi]
\] 或者按照指标记法 \[
\varphi_{I}^{h} = N_{I \alpha} \varphi_{\alpha}
\] 使用这种指标记法, 式 (2) 可以改写为 \[
M_{\alpha \beta} = \int_{B_0} \rho_0 N_{i \alpha} N_{i
\beta} ~\mathrm{d} X, \\
f_{\alpha}^{\mathrm{int}} = \int_{B_0} P_{iJ} \frac{ \partial N_{i
\alpha} }{ \partial X_J}
~\mathrm{d} X, \\
f_{\alpha}^{\mathrm{ext}} = \int_{B_0} \rho_0 b_{i} N_{i
\alpha} ~\mathrm{d} X
+ \int_{\partial B_0^t} \overline{T}_i N_{i \alpha} ~\mathrm{d}
\Sigma(X)
\tag{2'}
\]
B 矩阵
内力 \(F_{\alpha}^{\mathrm{int}}\) 在写成向量形式时, 因为式中 PK1 应力 \(P_{iJ}\) 是不对称的, 所以没有办法将应力张量 \(P\) 按照 Voigt 记法压缩成一列. 在这里我们将 PK1 转换为对称的 PK2 应力, 并且用式 (2) 的形式表示形函数 \[ f_{i}^{a,\mathrm{int}} \triangleq \int_{B_0} F_{iI} S_{IJ} \frac{ \partial N^a}{ \partial X_J} ~\mathrm{d} X \] 定义 B 矩阵为 \[ B_{(I,J)(a,i)} \triangleq \frac{ \partial x_i}{ \partial X_I} \frac{ \partial N^a}{ \partial X_J} \] 式中, \((a,i)\) 通过自由度编号映射压缩为 B 矩阵列的索引 \(\alpha\), \((I,J)\) 通过 Voigt 映射压缩为 B 矩阵行的索引 \(I'\). 固定节点 \(a\), 得到 B 矩阵的子矩阵 \([\mathrm{B}^{a}]\), 由此 B 矩阵可以写作 \[ [\mathrm{B}] = [[\mathrm{B}^{1}], [\mathrm{B}^{2}], \ldots, [\mathrm{B}^{\mathrm{nnd}}]] \] 单元内力写成矩阵形式为 \[ [f^{\mathrm{int}}] = \int_{B_0} [\mathrm{B}]^{\top} [\mathrm{S}] \] \([S]\) 是按照 Voigt 排布的列向量.
二维情景下, 矩阵 \([\mathrm{B}^{a}]\) 的维度为 3
x 2, 三维情景下的维度为 6 x 3.
对应 B 矩阵的维度分别为 3 x nedof 和
6 x nedof. B 矩阵的展开式为 \[
[\mathrm{B}^{a}] = \left[\begin{array}{cc}
\frac{\partial N^a}{\partial X} \frac{\partial x}{\partial X} &
\frac{\partial N^a}{\partial X} \frac{\partial y}{\partial X} \\
\frac{\partial N^a}{\partial Y} \frac{\partial x}{\partial Y} &
\frac{\partial N^a}{\partial Y} \frac{\partial y}{\partial Y} \\
\frac{\partial N^a}{\partial X} \frac{\partial x}{\partial
Y}+\frac{\partial N^a}{\partial Y} \frac{\partial x}{\partial X} &
\frac{\partial N^a}{\partial X} \frac{\partial y}{\partial
Y}+\frac{\partial N^a}{\partial Y} \frac{\partial y}{\partial X}
\end{array}\right]
\]
以及三维情景: \[ [\mathrm{B}^{a}] = \left[\begin{array}{ccc} \frac{\partial N^a}{\partial X} \frac{\partial x}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial y}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial z}{\partial X} \\ \frac{\partial N^a}{\partial Y} \frac{\partial x}{\partial Y} & \frac{\partial N^a}{\partial Y} \frac{\partial y}{\partial Y} & \frac{\partial N^a}{\partial Y} \frac{\partial z}{\partial Y} \\ \frac{\partial N^a}{\partial Z} \frac{\partial x}{\partial Z} & \frac{\partial N^a}{\partial Z} \frac{\partial y}{\partial Z} & \frac{\partial N^a}{\partial Z} \frac{\partial z}{\partial Z} \\ \frac{\partial N^a}{\partial Y} \frac{\partial x}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial x}{\partial Y} & \frac{\partial N^a}{\partial Y} \frac{\partial y}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial y}{\partial Y} & \frac{\partial N^a}{\partial Y} \frac{\partial z}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial z}{\partial Y} \\ \frac{\partial N^a}{\partial X} \frac{\partial x}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial x}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial y}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial y}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial z}{\partial Z}+\frac{\partial N^a}{\partial Z} \frac{\partial z}{\partial X} \\ \frac{\partial N^a}{\partial X} \frac{\partial x}{\partial Y}+\frac{\partial N^a}{\partial Y} \frac{\partial x}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial y}{\partial Y}+\frac{\partial N^a}{\partial Y} \frac{\partial y}{\partial X} & \frac{\partial N^a}{\partial X} \frac{\partial z}{\partial Y}+\frac{\partial N^a}{\partial Y} \frac{\partial z}{\partial X} \end{array}\right] \]
系统中的格林应变可通过 B 矩阵得到 \[ [\dot{\mathrm{E}}] = [\mathrm{B}] [\dot{\varphi}] \]
单元坐标系
式 (2) 的积分是在参考构型的单元上进行, 可以通过积分换元, 将单元上的积分转换到同一个标准形状的单元上进行. 如果用 \(\xi\) 表示单元坐标, 那么微元之间的映射关系为 \[ ~\mathrm{d} X = \frac{ \partial X }{ \partial \xi } \cdot \mathrm{d} \xi \triangleq F_{\xi}^{0} \cdot \mathrm{d} \xi \] 因此式 (2) 转换为 \[ M_{\alpha \beta} = \int_{\square} \rho_0 N_{I \alpha} N_{I \beta} F_{\xi}^{0} ~\mathrm{d} \xi, \\ F_{\alpha}^{\mathrm{int}} = \int_{\square} S_{IJ} \frac{ \partial N_{I \alpha} }{ \partial \xi_K} \left( \frac{ \partial X_K }{ \partial \xi_J } \right)^{-1} ~\mathrm{d} \xi, \\ F_{\alpha}^{\mathrm{ext}} = \int_{\square} \rho_0 f_{I} N_{I \alpha} F_{\xi}^{0} ~\mathrm{d} \xi + \int_{\partial \Omega_0^t} \overline{T}_I N_{I \alpha} ~\mathrm{d} S_0 \tag{2'} \]