三节点膜单元列式

发表于 更新于 分类于 阅读次数：

参考文献：[1]TAYLOR R L, OÑATE E, UBACH P A. Finite Element Analysis of Membrane Structures[M/OL]//OÑATE E, KRÖPLIN B. Textile Composites and Inflatable Structures: 卷 3. Berlin/Heidelberg: Springer-Verlag, 2005: 47-68[2024-05-24]. http://link.springer.com/10.1007/1-4020-3317-6_4. DOI:10.1007/1-4020-3317-6_4.

单元的插值函数

考虑三节点的膜单元, 如图所示

image-20240524182715120

在大变形的框架下推导膜单元的运动学表述, 参考构型中的物质点使用 \(\boldsymbol{X}\) 表示, 当前构型中的物质点使用 \(\boldsymbol{x}\) 表示. 物质点 \(\boldsymbol{X}\) 从参考构型中经过位移 \(\boldsymbol{u}\) 运动到点 \(\boldsymbol{x}\): \[ \boldsymbol{x} = \boldsymbol{X} + \boldsymbol{u} \] 单元内部的点通过节点线性插值得到, 也即 \[ \boldsymbol{x} = N^{a} \tilde {\boldsymbol{x}}^a, \quad \boldsymbol{X} = N^{a} \tilde{\boldsymbol{X}}^a \] 式中, \(N^a\) 为三角形单元节点处的形函数, \(\tilde \square\) 表示节点处的物理量. 上式以及在接下来的陈述中默认对成对的节点指标进行求和, 求和范围为 \(1,2,\ldots,\mathrm{nen}\). \(\mathrm{nen}\) 为单元节点数, 本文中三节点膜单元 \(\mathrm{nen}=3\). 三角形三节点单元, 在自然坐标系下的形函数为 \[ N^{1}(\xi_1, \xi_2) = 1 - \xi_1 - \xi_2, \quad N^{2}(\xi_1, \xi_2) = \xi_1, \quad N^{3}(\xi_1, \xi_2) = \xi_2 \] 待求解的未知量为变形映射函数 \(\boldsymbol{\varphi}\): \[ \boldsymbol{x} = \boldsymbol{\varphi}(\boldsymbol{X}, t) \] > [!NOTE] > > 在本文中, \(x\) 和 \(\varphi\) 是同一个函数. 一般使用 \(\delta \varphi\) 表示对 \(x\) 的变分.

膜单元的面内坐标系

上述列式是在一般的三维空间中确定了膜单元的运动过程. 但我们更关心膜单元面内的运动, 所以定义单元的面内坐标系 \(y\) 和 \(Y\)​​. 在引入面内坐标系之前, 需要约定单元的编号顺序与单元法向方向. 使用 ABAQUS 的惯例如下: 编号顺序与单元正法向符合右手螺旋方式:

img

在此基础之上, 定义面内坐标1 方向单位向量 \(\boldsymbol{m}_1\) 与三角形单元边 \(12\) 重合, 并指向节点 \(2\): \[ \boldsymbol{m}_1 \triangleq \frac{ \tilde{\boldsymbol{x}}^{21}}{\| \tilde{\boldsymbol{x}}^{21} \|}, \quad \tilde{\boldsymbol{x}}^{ij} \triangleq \tilde{\boldsymbol{x}}^i - \tilde{\boldsymbol{x}}^j, \quad i,j=1,2,3 \] 单元的法向 \(\boldsymbol{m}_3\) 由边 \(12\), \(13\) 叉乘得到: \[ \boldsymbol{m}_3 \triangleq \frac{ \tilde{\boldsymbol{x}}^{21} \times \tilde{\boldsymbol{x}}^{31}} {\| \tilde{\boldsymbol{x}}^{21} \times \tilde{\boldsymbol{x}}^{31}\|} \] 边 \(12\), \(13\) 叉乘得到的向量模长为膜单元面积 \(A\) 的 2 倍: \[ 2A \triangleq \| \tilde{\boldsymbol{x}}^{21} \times \tilde{\boldsymbol{x}}^{31} \| \] 面内坐标的另一方向 \(\boldsymbol{m}_2\) 可通过 \(\boldsymbol{m}_3\) 与 \(\boldsymbol{m}_1\) 之间叉乘得到 \[ \boldsymbol{m}_2 \triangleq \boldsymbol{m}_3 \times \boldsymbol{m}_1 \] 同样的过程可以定义参考构型下的面内坐标系 \(\{\boldsymbol{M}_1, ~\boldsymbol{M}_2, ~ \boldsymbol{M}_3 \}\).

面内坐标与全局坐标的变换关系

当面内坐标系的原点放在节点 \(1\) 处时, 全局坐标 \(x\), \(X\) 到单元面内坐标 \(y\), \(Y\) 的变换关系为 \[ Y_i = (\boldsymbol{X} - \tilde{\boldsymbol{X}}^1) \cdot \boldsymbol{M}_i , \quad y_i = (\boldsymbol{x} - \tilde{\boldsymbol{x}}^1) \cdot \boldsymbol{m}_i \] 如果定义坐标变换矩阵 \([\mathrm{Q}]\), \([\mathrm{q}]\) 为 \[ [\mathrm{Q}] \triangleq [\mathbf{M}_1, \mathbf{M}_2, \mathbf{M}_3], \quad [\mathrm{q}] \triangleq [\mathbf{m}_1, \mathbf{m}_2, \mathbf{m}_3] \] 式中, \(\mathbf{M}_i\), \(\mathbf{m}_i\) 是在全局坐标系下表示的列向量. 式 (3) 可以改写成如下形式 \[ Y_{i} = Q_{mi}(X_m - \tilde{X}_m^1), \quad y_{i} = q_{mi}(x_m - \tilde{x}_m^1) \] 因此得到变换公式为

关于坐标的导数变换 \[ \frac{ \partial ~ }{ \partial Y_i } = Q_{mi} \frac{ \partial ~ }{ \partial X_m }, \quad \frac{ \partial ~ }{ \partial y_i } = q_{mi} \frac{ \partial ~ }{ \partial x_m } \] 节点坐标之间的变换 \[ \begin{equation} \tilde Y_i^a = \tilde{\boldsymbol{X}}^{a1} , \quad \tilde y_i^a = \tilde{\boldsymbol{x}}^{a1} \cdot \boldsymbol{m}_i \quad a = 1,2,3 \quad i = 1,2 \label{eq:coord_transform} \end{equation} \] 在新坐标系 \(y\) 下, 单元内的点可以表示为 \[ \boldsymbol{Y} = N^a \tilde{\boldsymbol{Y}}^a, \quad \boldsymbol{y} = N^a \tilde{\boldsymbol{y}}^a \] 注意到, 此时节点 \(1\) 在坐标系 \(y\) 中的坐标为 \((0,0)\), 所以上式展开为 \[ \boldsymbol{Y} = N^2 \tilde{\boldsymbol{Y}}^2 + N^3 \tilde{\boldsymbol{Y}}^3, \quad \boldsymbol{y} = N^2 \tilde{\boldsymbol{y}}^2 + N^3 \tilde{\boldsymbol{y}}^3 \]

应变度量

面内坐标系下的变形梯度定义为 \[ F_{iI} \triangleq \frac{\partial y_{i}}{\partial Y_{I}} \] 注意到在面内坐标系下, \(y_3, Y_3 \equiv 0\), 因此上述指标的范围为 \(i,I = 1,2\).

[!NOTE]

在全局坐标系下, 变形梯度为 \[ F_{jJ} \triangleq \frac{\partial x_j}{\partial X_J}, \quad j,J = 1,2,3 \] 在矩阵记法下, \([\mathrm{F}]\) 是 \(3 \times 3\) 的矩阵. 可以通过坐标变换得到面内坐标系下的变形梯度转换关系 \[ F_{jJ} = \frac{\partial x_j}{\partial y_i} \frac{\partial y_i}{\partial Y_I} \frac{\partial Y_I}{\partial X_J} = q_{ji} Q_{JI} \frac{\partial y_i}{\partial Y_I}, \quad j,J = 1,2,3, \quad i,I = 1,2 \]

在自然坐标系下, 变形梯度可以通过链式法则得到 \[ F_{iI} = \frac{\partial y_i}{\partial \xi_k} \frac{\partial \xi_k}{\partial Y_I} \] 如果记不同构型下对自然坐标的 Jacobian 变换 \(\partial \square / \partial \xi\) 分别为 \[ j \triangleq \frac{\partial y}{\partial \xi}, \quad J \triangleq \frac{\partial Y}{\partial \xi} \] 那么, 变形梯度可写作 \[ F = j \cdot J^{-1} \] 对应的右 Cauchy-Green 变形张量为 \[ \begin{equation} C \triangleq F^{\top} \cdot F = J^{-\top} \cdot j^{\top} \cdot j \cdot J^{-1} \label{eq:C_rep_as_jacobian} \end{equation} \] 将上式写作矩阵形式, 并注意到, 对于线性插值的三角形三节点单元, Jacobian 为常数, 为面内坐标系的节点坐标: \[ [\mathrm{j}] = \begin{bmatrix} \tilde y_1^2 & \tilde y_1^3 \\ \tilde y_2^2 & \tilde y_2^3 \end{bmatrix}, \quad [\mathrm{J}] = \begin{bmatrix} \tilde Y_1^2 & \tilde Y_1^3 \\ \tilde Y_2^2 & \tilde Y_2^3 \end{bmatrix} \] 再根据式 \(\ref{eq:coord_transform}\) 给出的节点坐标关系, 可以得到 \[ [\mathrm{j}] = \begin{bmatrix} \tilde{\boldsymbol{x}}^{21} \cdot \boldsymbol{m}_1 & \tilde{\boldsymbol{x}}^{31} \cdot \boldsymbol{m}_1 \\ \tilde{\boldsymbol{x}}^{21} \cdot \boldsymbol{m}_2 & \tilde{\boldsymbol{x}}^{31} \cdot \boldsymbol{m}_2 \end{bmatrix} \] 在面内坐标系下, 节点 2 的 \(2\) 方向分量恒等于 0, 因此 Jacobian 矩阵可以进一步简化为 \[ [\mathrm{j}] = \begin{bmatrix} \| \tilde{\boldsymbol{x}}^{21} \| & \tilde{\boldsymbol{x}}^{31} \cdot \boldsymbol{m}_1 \\ 0 & \tilde{\boldsymbol{x}}^{31} \cdot \boldsymbol{m}_2 \end{bmatrix} \] 类似的, 参考构型中的 Jacobian 矩阵可以写作 \[ [\mathrm{J}] = \begin{bmatrix} \| \tilde{\boldsymbol{X}}^{21} \| & \tilde{\boldsymbol{X}}^{31} \cdot \boldsymbol{M}_1 \\ 0 & \tilde{\boldsymbol{X}}^{31} \cdot \boldsymbol{M}_2 \end{bmatrix} \] 在给出参考构型与当前构型中的 Jacobian 矩阵之后, 就可以将式 \(\ref{eq:C_rep_as_jacobian}\) 中的右 Cauchy-Green 变形张量表示为 \[ C = G^\top g ~G \] 式中, \(G\) 为参考构型 Jacobian 变换的逆, \(G \triangleq J^{-1}\), \(g \triangleq j^\top j\). 将它们写成矩阵形式为 \[ \begin{equation} [\mathrm{G}] = \frac{1}{J_{11} J_{22}} \begin{bmatrix} J_{22} & -J_{12} \\ 0 & J_{11} \end{bmatrix}, \quad [\mathrm{g}] = \begin{bmatrix} j_{11}^2 & j_{11}j_{12} \\ j_{11}j_{12} & j_{12}^2 + j_{22}^2 \end{bmatrix} = \begin{bmatrix} \tilde{\boldsymbol{x}}^{21} \cdot \tilde{\boldsymbol{x}}^{21} & \tilde{\boldsymbol{x}}^{21} \cdot \tilde{\boldsymbol{x}}^{31} \\ \tilde{\boldsymbol{x}}^{21} \cdot \tilde{\boldsymbol{x}}^{31} & \tilde{\boldsymbol{x}}^{31} \cdot \tilde{\boldsymbol{x}}^{31} \end{bmatrix} \label{eq:G_and_g} \end{equation} \] 所以, 张量 \(C\) 最终可以表示为如下矩阵形式 \[ [\mathrm{C}] = [\mathrm{G}]^\top [\mathrm{g}] [\mathrm{G}] \]

Green 应变定义为 \[ E \triangleq \frac{1}{2}(C - I) \]

本构方程

St.Venant-Kirchhoff 材料

使用 St.Venant-Kirchhoff 材料模型, 膜单元上的应力 \(S\) 为 \[ \begin{equation} S = \mathbb{L} : E \label{eq:constitutive} \end{equation} \] 式中, \(\mathbb{L}\) 为弹性模量.

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\). 其它应变不变量分别为 \[ \begin{equation} \begin{aligned} \bar{\mathscr{I}}_1 &= \frac{\mathscr{I}_1}{J^{2/3}}, \quad \mathscr{I}_1 = \mathrm{tr} B \\ \bar{\mathscr{I}}_2 &= \frac{\mathscr{I}_2}{J^{4/3}}, \quad \mathscr{I}_2 = \frac{1}{2} \left( (\mathrm{tr}B)^2 - \mathrm{tr}B^2 \right)\\ J &= \det F \end{aligned} \label{eq:cons_rep_as_B_inva_3} \end{equation} \] 将 \(W\) 对不变量的导数带入到式 \(\ref{eq:cons_rep_as_B_inva_3}\) 中得到 \[ \begin{equation} \boxed{ \boldsymbol{\sigma} = \frac{\mu}{J^{5/3}} ~\mathrm{dev}\boldsymbol{B} + K(J-1)\boldsymbol{I}, \quad \mathrm{dev} \boldsymbol{B} = \boldsymbol{B} - \frac{1}{3} \mathscr{I}_1 \boldsymbol{I} } \label{eq:neo_hookean_cttt} \end{equation} \]

控制方程弱形式

膜单元的控制方程的弱形式由下式给出: \[ \begin{equation} \boxed{ \begin{aligned} \delta \Pi & = \int_{S_0} \delta \varphi_i \rho_0 \ddot{\varphi}_i h \mathrm{~d} X + \int_{S_0} \delta \varphi_i c_0 \dot{\varphi}_i h \mathrm{~d} X + \int_{S_0} \delta E_{I J} S_{I J} h \mathrm{~d} X \\ &- \int_{S} \delta \varphi_i n_i p \mathrm{~d} x - \int_{\partial S_{t}} \delta \varphi_i \bar{t}_i \mathrm{~d} \sigma(x) = 0 \end{aligned} } \label{eq:weak_form} \end{equation} \] 式中, 积分区域 \(S\) 和 \(S_0\) 分别指代当前构型和参考构型的膜单元区域. \(h\) 表示膜的厚度, 在本问题中考虑为常数. \(n_i\) 为当前构型下的面法向, \(p\) 为内压力. \(\partial S_{t}\) 表示力边界, \(\bar t_i\) 表示施加在力边界上的作用力. 前三项积分在参考构型中得到, 后两项积分在当前构型中得到.

上式中内力虚功使用的是功共轭的 PK2 应力 \(S_{IJ}\) 和 Green 应变 \(E_{IJ}\). 对于 Neo-Hookean 材料, 本构关系由 Cauchy 应力和变形梯度给出. 因此使用功共轭的 Kirchhoff 应力 $ _{ij}$ 和 变形率张量 \(D_{ij}\) (在参考构型下积分) 更方便一些: \[ \begin{equation} \delta W^{\mathrm{int}} = \int_{S_0} \delta D_{ij} \tau_{ij} ~\mathrm{d} X \label{eq:tau-D} \end{equation} \]

内力虚功

\(E\) - \(S\) 共轭对

考察积分项 \(\delta E_{IJ}S_{IJ}\), \[ \delta E_{IJ} S_{IJ} = \delta \tilde{\boldsymbol{\varphi}}^a \cdot \frac{ \partial E_{IJ}}{ \partial \tilde{\boldsymbol{\varphi}}^a} S_{IJ} = [\delta \varphi]^{\top} [\mathrm{B}]^{\top} [\mathrm{S}] \] 将单元上的自由度变量 \(\tilde{\boldsymbol{\varphi}}^a\) 排列为列向量, 对称的 PK2 应力按照 Voigt 记法记作列向量 \[ [\delta \varphi] = \begin{bmatrix} \delta \tilde{\boldsymbol{\varphi}}^{1} \\ \delta \tilde{\boldsymbol{\varphi}}^{2} \\ \delta \tilde{\boldsymbol{\varphi}}^{3} \end{bmatrix}_{9 \times 1}, \quad [\mathrm{S}] = \begin{bmatrix} S_{11}\\ S_{22}\\ S_{12} \end{bmatrix}_{3 \times 1} \] 在定义矩阵 \([\mathrm{S}]\) 和 \([\delta \varphi]\) 的形式之后, \([\mathrm{B}]\) 矩阵的形式也就唯一确定了: \[ [\mathrm{B}] = \begin{bmatrix} \dfrac{ \partial E_{11} }{ \partial \tilde{\boldsymbol{\varphi}}^1 } & \dfrac{ \partial E_{11} }{ \partial \tilde{\boldsymbol{\varphi}}^2 } & \dfrac{ \partial E_{11} }{ \partial \tilde{\boldsymbol{\varphi}}^3 } \\ \dfrac{ \partial E_{22} }{ \partial \tilde{\boldsymbol{\varphi}}^1 } & \dfrac{ \partial E_{22} }{ \partial \tilde{\boldsymbol{\varphi}}^2 } & \dfrac{ \partial E_{22} }{ \partial \tilde{\boldsymbol{\varphi}}^3 } \\ 2\dfrac{ \partial E_{12} }{ \partial \tilde{\boldsymbol{\varphi}}^1 } & 2\dfrac{ \partial E_{12} }{ \partial \tilde{\boldsymbol{\varphi}}^2 } & 2\dfrac{ \partial E_{12} }{ \partial \tilde{\boldsymbol{\varphi}}^3 } \end{bmatrix}_{3 \times 9} \] 继续写出 \([\mathrm{B}]\) 矩阵分量的具体形式. 因为变分符号只作用在与函数变量 \(\varphi\) 相关的变量上, 所以 \[ \delta E_{IJ} S_{IJ} = \frac{1}{2} \delta C_{IJ} S_{IJ} = \frac{1}{2} \delta g_{ij} G_{iI} G_{jJ} S_{IJ} \] 若定义应力 \(s\) 为 \[ s_{ij} \triangleq G_{iI} G_{jJ} S_{IJ} \] 并将 \(s\) 写成如下矩阵形式 \[ [\mathrm{s}] = \begin{bmatrix} s_{11}\\ s_{22}\\ s_{12} \end{bmatrix} \] 那么, 对应的映射矩阵 \([\mathrm{M}]\) 以及分量为 \[ \boxed{ [\mathrm{s}] = [\mathrm{M}]^{\top} [\mathrm{S}], \quad [\mathrm{M}] = \begin{bmatrix} 1/J_{11}^2 & 0 & 0 \\ 1/J_{22}^2 & 1/J_{22}^2 & -J_{12}/J_{11}J_{22}^2 \\ -2J_{12}/J_{11}^2J_{22} & 0 & 1/J_{11}J_{22} \\ \end{bmatrix} } \] 将 \(\delta g_{ij}\) 排列为列向量的形式: \[ [\delta g] = \begin{bmatrix} \delta g_{11} \\ \delta g_{22} \\ 2\delta g_{12} \end{bmatrix} \] 因此, 功共轭项可以写成如下矩阵形式 \[ \delta E_{IJ} S_{IJ} = \frac{1}{2} [\delta g]^{\top} [\mathrm{s}] \] 对 \(g_{ij}\) 的变分转换为对节点位移的变分. 根据式 \(\ref{eq:G_and_g}\), \(g\) 的分量可以表示为 \[ g_{11} = \tilde{\boldsymbol{x}}^{21} \cdot \tilde{\boldsymbol{x}}^{21}, \quad g_{22} = \tilde{\boldsymbol{x}}^{31} \cdot \tilde{\boldsymbol{x}}^{31}, \quad g_{12} = g_{21} = \tilde{\boldsymbol{x}}^{21} \cdot \tilde{\boldsymbol{x}}^{31} \] 因此, \[ \begin{equation} \begin{aligned} \delta g_{11} &= 2(\delta \tilde{\boldsymbol{\varphi}}^2 - \delta \tilde{\boldsymbol{\varphi}}^1) \cdot \tilde{\boldsymbol{x}}^{21} &= -2\tilde{\boldsymbol{x}}^{21} \cdot \delta \tilde{\boldsymbol{\varphi}}^1 &+ 2\tilde{\boldsymbol{x}}^{21} \cdot \delta \tilde{\boldsymbol{\varphi}}^2 & \\ \delta g_{22} &= 2(\delta \tilde{\boldsymbol{\varphi}}^3 - \delta \tilde{\boldsymbol{\varphi}}^1) \cdot \tilde{\boldsymbol{x}}^{31} &= -2\tilde{\boldsymbol{x}}^{31} \cdot \delta \tilde{\boldsymbol{\varphi}}^1 &&+ 2\tilde{\boldsymbol{x}}^{31} \cdot \delta \tilde{\boldsymbol{\varphi}}^3 \\ 2\delta g_{12} &= 2(\delta \tilde{\boldsymbol{\varphi}}^2 - \delta\tilde{\boldsymbol{\varphi}}^1) \cdot \tilde{\boldsymbol{x}}^{31} + 2(\delta \tilde{\boldsymbol{\varphi}}^3 - \delta \tilde{\boldsymbol{\varphi}}^1) \cdot \tilde{\boldsymbol{x}}^{21} &= -2(\tilde{\boldsymbol{x}}^{21} + \tilde{\boldsymbol{x}}^{31}) \cdot \delta \tilde{\boldsymbol{\varphi}}^1 &+ 2\tilde{\boldsymbol{x}}^{31} \cdot \delta \tilde{\boldsymbol{\varphi}}^2 &+ 2\tilde{\boldsymbol{x}}^{21} \cdot \delta \tilde{\boldsymbol{\varphi}}^3 \end{aligned} \label{eq:variation_g} \end{equation} \] 写成矩阵形式 \[ \boxed{ [\delta g] = 2 \begin{bmatrix} -\tilde{\boldsymbol{x}}^{21} & \tilde{\boldsymbol{x}}^{21} & 0 \\ -\tilde{\boldsymbol{x}}^{31} & 0 & \tilde{\boldsymbol{x}}^{31} \\ -\tilde{\boldsymbol{x}}^{21} - \tilde{\boldsymbol{x}}^{31} & \tilde{\boldsymbol{x}}^{31} & \tilde{\boldsymbol{x}}^{21} \end{bmatrix}_{3\times 9} \begin{bmatrix} \delta \tilde{\boldsymbol{\varphi}}^{1} \\ \delta \tilde{\boldsymbol{\varphi}}^{2} \\ \delta \tilde{\boldsymbol{\varphi}}^{3} \end{bmatrix}_{9 \times 1} \triangleq 2 [\mathrm{b}] [\delta \varphi] } \] 式中, \(\tilde{\boldsymbol{x}}^{ab}\) 是行向量. 积分项 \(\delta E_{IJ}S_{IJ}\) 写成矩阵形式为 \[ \delta E_{IJ}S_{IJ} = [\delta \varphi]^{\top} [\mathrm{b}]^\top [\mathrm{M}]^{\top} [\mathrm{S}] \]

定义单元上的 \([\mathrm{B}]\) 矩阵: \[ [\mathrm{B}] = [\mathrm{M}][\mathrm{b}] \] 最终得到 \[ \delta E_{IJ}S_{IJ} = [\delta \varphi]^{\top} ~[\mathrm{B}]^\top [\mathrm{S}] \] 对于线性插值得到的位移场函数, 膜单元上的应变为常值, 又因为材料本构方程 \(\ref{eq:constitutive}\) 同样是常数, 所以单元上的应力也是常数, 所以, 弱形式 \(\ref{eq:weak_form}\) 中内力虚功项为 \[ \begin{equation} \boxed{ [\mathrm{f}^{\mathrm{int}}] \triangleq h A_0 [\mathrm{B}]^\top [\mathrm{S}], \quad \delta W^{\mathrm{int}} = \int_{S_0} \delta E_{I J} S_{I J} h \mathrm{~d} X = [\delta \varphi]^{\top} [\mathrm{f}^{\mathrm{int}}] } \label{eq:work_int} \end{equation} \] 式中, \(A_0\) 为参考构型中单元的面积.

\(D\) - \(\tau\) 共轭对

将插值函数代入到式 \(\ref{eq:tau-D}\) 中得到内力项: \[ \begin{equation} f_{i}^{\mathrm{int},a} = \int_{S_0} \tau_{ij} N_{,I}^a F_{Ij}^{-1} ~\mathrm{d} X \label{eq:f_int_tau_D} \end{equation} \]

惯性力和粘力虚功

式 \(\ref{eq:weak_form}\) 中的惯性力虚功为 \[ \delta W^{\mathrm{kin}} = \int_{S_0} \delta \varphi_i \rho_0 h \ddot{\varphi}_i ~\mathrm{d} X \] 将插值形函数代入得到 \[ \delta W^{\mathrm{kin}} = \delta \varphi_i^a \int_{S_0} \rho_0 h N^a N^b ~\mathrm{d} X ~\ddot{\varphi}_i^b \triangleq \delta \varphi_i^a M^{ab} \ddot{\varphi}_i^b \] 因此惯性力虚功写成矩阵形式为 \[ \begin{equation} \boxed{ \begin{gathered} \underbrace{[\mathrm{M}]}_{\rm 3 nen \times 3 nen} = [\mathrm{M}^{ab}], \quad \underbrace{[\mathrm{M}^{ab}]}_{3 \times 3} = 2 h A_0 \int_{\triangle} \rho_0 N^a N^b ~\mathrm{d} \xi \underbrace{[\mathrm{I}]}_{3 \times 3} \\ \delta W^{\mathrm{kin}} = [\delta \varphi]^{\top} [\mathrm{M}] [\ddot{\varphi}] \end{gathered} } \label{eq:work_kin} \end{equation} \] > [!TIP] > > 如果材料的密度为常数, 那么积分得到单元的质量矩阵为 > \[ > [\mathrm{M}]_{9 \times 9} = \rho_0 A_0 h \begin{bmatrix} > \frac{1}{6} \mathrm{I} & \frac{1}{12} \mathrm{I} & \frac{1}{12} \mathrm{I} \\ > \frac{1}{12} \mathrm{I} & \frac{1}{6} \mathrm{I} & \frac{1}{12} \mathrm{I} \\ > \frac{1}{12} \mathrm{I} & \frac{1}{12} \mathrm{I} & \frac{1}{6} \mathrm{I} > \end{bmatrix} > \] > 一般使用 lump matrix, 矩阵 \([\mathrm{M}]\) 对角线元素为每一行之和, 其它元素置为 0: > \[ > [\tilde{\mathrm{M}}]_{9 \times 9} = \frac{\rho_0 A_0 h}{3} \begin{bmatrix} > \mathrm{I} && \\ > & \mathrm{I} & \\ > && \mathrm{I} > \end{bmatrix} > \]

类似的, 粘力虚功和粘性矩阵为 \[ \begin{equation} \boxed{ \begin{gathered} \underbrace{[\mathrm{C}]}_{\rm 3 nen \times 3 nen} = [\mathrm{C}^{ab}], \quad \underbrace{[\mathrm{C}^{ab}]}_{3 \times 3} = 2 h A_0 \int_{\triangle} c_0 N^a N^b ~\mathrm{d} \xi \underbrace{[\mathrm{I}]}_{3 \times 3}\\ \delta W^{\mathrm{vic}} = [\delta \varphi]^{\top} [\mathrm{C}] [\dot{\varphi}] \end{gathered} } \end{equation} \]

外力虚功

式 \(\ref{eq:weak_form}\) 中的外力虚功项为 \[ \delta W^{\mathrm{ext}} = \int_{S} \delta \varphi_i n_i p \mathrm{~d} x + \int_{\partial S_{t}} \delta \varphi_i \bar{t}_i \mathrm{~d} \sigma(x) \] 第一项是作用在膜单元的 SNEG 面压强造成的外力: \[ \int_{S} \delta \varphi_i n_i p \mathrm{~d} x = \delta \boldsymbol{\varphi}^{a} \cdot \int_{S} N^a p \boldsymbol{n} ~\mathrm{d}x \triangleq [\delta \varphi^{a}]^{\top} [\mathrm{f}_{p}^{a}] \] 对于三节点膜单元, 单元的法向 \(\boldsymbol{n}\) 在单元内是常数, 因此单元节点上的外力指向与面法向一致. 同时如果假设 \(p\) 也是常数, 那么单元节点上由压力造成的外力相等: \[ \boldsymbol{f}_{p}^{a} = \frac{1}{3} p A \boldsymbol{n} \] 式中, \(A\) 为变形后的单元面积. 注意到 \(A \boldsymbol{n}\) 可以通过变形后单元边向量叉乘得到: \[ \boldsymbol{f}_{p}^{a} = \frac{1}{6} p \tilde{\boldsymbol{x}}^{21} \times \tilde{\boldsymbol{x}}^{31} \] 上式中的叉乘也可以写成矩阵的形式, 构造 \(3\times 3\)矩阵 \([\widehat{\tilde{x}^{21}}]\) \[ [\widehat{\tilde{x}^{21}}] = \begin{bmatrix} 0 & -\tilde{x}^{21}_3 & \tilde{x}^{21}_2 \\ \tilde{x}^{21}_3 & 0 & -\tilde{x}^{21}_1 \\ \tilde{x}^{21}_3 & -\tilde{x}^{21}_2 & 0 \end{bmatrix} \] 并将 \(\tilde{\boldsymbol{x}}^{31}\) 写作列向量的形式 \([\tilde{x}^{31}]\), 因此有 \[ \begin{equation} \boxed{ \begin{gathered} \underbrace{[\mathrm{f}^{\mathrm{ext}}]}_{\rm 3 nen \times 1} = [\mathrm{f}_{p}^{a}], \quad \underbrace{[\mathrm{f}_{p}^{a}]}_{3 \times 1} = \frac{1}{6} p [\widehat{\tilde{x}^{21}}] [\tilde{x}^{31}] \\ \delta W^{\mathrm{ext}} = [\delta \varphi]^{\top} [\mathrm{f}^{\mathrm{ext}}] \end{gathered} } \label{eq:work_ext} \end{equation} \]

空间离散方程组

在使用形函数对空间离散之后, 结合式 \(\ref{eq:work_int}\) ~ \(\ref{eq:work_ext}\), 就得到 \[ \delta W^{\mathrm{kin}} + \delta W^{\mathrm{vic}} + \delta W^{\mathrm{int}} - \delta W^{\mathrm{ext}} = 0, \quad \forall~ \delta \varphi \\ \Rightarrow [\delta \varphi]^{\top} \left( [\mathrm{M}] [\ddot{\varphi}] + [\mathrm{C}] [\dot{\varphi}] + h A_0 [\mathrm{B}]^\top [\mathrm{S}] - [\mathrm{f}_{p}] \right) = 0, \quad \forall ~ [\delta \varphi] \] 由 \([\delta \varphi]\) 在未约束自由度的任意性, 得到关于节点自由自由度处的时域的常微分方程组: \[ \begin{equation} \boxed{ [\mathrm{M}] [\ddot{\varphi}] + [\mathrm{C}] [\dot{\varphi}] + h A_0 [\mathrm{B}]^\top [\mathrm{S}] - [\mathrm{f}_{p}] = 0 } \label{eq:ode} \end{equation} \]

隐式求解方法

Newmark 方法 – 瞬态问题

对于瞬态问题, 一般使用 Newmark 方法进行求解. 设在 \(t_n\) 时刻的变量 \(\{ \varphi_{n},~ \dot{\varphi}_n, ~ \ddot{\varphi}_n\}\) 已知, 根据中值定理: \[ \frac{\dot{\varphi}_{n+1} - \dot{\varphi}_n}{\Delta t} = \ddot{\varphi}_{n+\xi}, \quad \xi \in [0,1] \] 假设 \(\ddot{\varphi}_{n+\xi}\) 在 \(\ddot{\varphi}_{n}\) 和 \(\ddot{\varphi}_{n+1}\) 之间变化, 引入参数 \(\gamma\), 对 \(\ddot{\varphi}_{n}\) 和 \(\ddot{\varphi}_{n+1}\) 插值得到 \(\ddot{\varphi}_{n+\xi}\): \[ \dot{\varphi}_{n+1} = \dot{\varphi}_n + \Delta t \ddot{\varphi}_{\gamma}, \quad \ddot{\varphi}_{\gamma} = (1 - \gamma) \ddot{\varphi}_{n} + \gamma \ddot{\varphi}_{n+1}, \quad \gamma \in [0,1] \] 确定 \(\gamma\) 后就确定了下一时刻的速度更新方程. 下一时刻的位移项通过修正的加速度项得到: \[ \varphi_{n+1} = \varphi_n + \Delta t \dot{\varphi}_{n} + \frac{1}{2} \Delta t^2 \ddot{\varphi}_{\beta}, \quad \ddot{\varphi}_{\beta} = (1 - 2 \beta) \ddot{\varphi}_{n} + 2 \beta \ddot{\varphi}_{n+1}, \quad 2\beta \in [0,1] \] 再根据微分方程 \(\ref{eq:ode}\) 在时刻 \(t_{n+1}\) 得到的等式, 就得到关于 \(\{ \varphi_{n+1},~ \dot{\varphi}_{n+1}, ~ \ddot{\varphi}_{n+1}\}\) 的非线性代数方程组: \[ \boxed{ \begin{aligned} M \ddot{\varphi}_{n+1} &= f^{\mathrm{ext}}(\varphi_{n+1}, t_{n+1}) - f^{\mathrm{int}}(\varphi_{n+1}, t_{n+1}) -C \dot{\varphi}_{n+1} \\ \dot{\varphi}_{n+1} &= \dot{\varphi}_n + (1 - \gamma) \Delta t \ddot{\varphi}_{n} + \gamma \Delta t \ddot{\varphi}_{n+1} \\ \varphi_{n+1} &= \varphi_n + \Delta t \dot{\varphi}_{n} + (\frac{1}{2} - \beta) \Delta t^2 \ddot{\varphi}_{n} + \beta \Delta t^2 \ddot{\varphi}_{n+1} \end{aligned} } \] 整理得到关于 \(\varphi_{n+1}\) 的非线性方程组: \[ \boxed{ \begin{aligned} K^{\mathrm{nmk}} \varphi_{n+1} &+ f^{\mathrm{int}}(\varphi_{n+1}, t_{n+1}) - f^{\mathrm{ext}}(\varphi_{n+1}, t_{n+1}) \\ &= \underbrace{ K^{\mathrm{nmk}} \left( \varphi_{n} + \Delta t \dot{\varphi}_n + (\frac{1}{2} - \beta) \Delta t^2 \ddot{\varphi}_{n} \right) -C \left( \dot{\varphi}_n + (1 - \gamma) \Delta t \ddot{\varphi}_{n} \right) }_{\triangleq f_n^{\mathrm{nmk}}} \end{aligned} } \] 式中, \(K^{\mathrm{nmk}}\) 是根据 Newmark 算法产生的刚度矩阵, 维度和质量矩阵相同; \(f_n^{\mathrm{nmk}}\) 通过时刻 \(t_n\) 的变量 \(\{ \varphi_{n},~ \dot{\varphi}_n, ~ \ddot{\varphi}_n\}\) 计算得到. \[ \begin{equation} \boxed{ K^{\mathrm{nmk}} \triangleq \frac{1}{\beta \Delta t } (\frac{1}{\Delta t} M + \gamma C) } \label{eq:stiff_nmk} \end{equation} \] 定义 \(t_{n+1}\) 时刻的残差函数 \(R_{n+1}: \mathbb{R}^{\mathrm{ndof}} \mapsto \mathbb{R}^{\mathrm{ndof}}\): \[ \begin{equation} R_{n+1}(\varphi_{n+1} ) \triangleq f^{\mathrm{int}}(\varphi_{n+1}, t_{n+1}) - f^{\mathrm{ext}}(\varphi_{n+1}, t_{n+1}) + K^{\mathrm{nmk}} \varphi_{n+1} - f_{n}^{\mathrm{nmk}} \label{eq:residual} \end{equation} \] 式中, \(f^{\mathrm{int}}\) 和 \(f^{\mathrm{ext}}\) 是方程非线性的来源, 使用 Newton 迭代法进行求解. 设第 \(i\) 步迭代初始值为 \(\varphi_{n+1}^{(i)}\), 在该位置对残差函数进行 Taylor 展开至一阶项: \[ R_{n+1}(\varphi_{n+1}) = R_{n+1}(\varphi_{n+1}^{(i)}) + \left. \frac{\partial R_{n+1}}{\partial \varphi_{n+1}} \right|_{\varphi_{n+1}^{(i)}} \cdot (\varphi_{n+1} - \varphi_{n+1}^{(i)}) \] 记 \[ \begin{equation} K_{n+1}^{(i)} \triangleq \left. \frac{\partial R_{n+1}}{\partial \varphi_{n+1}} \right|_{\Delta \varphi_{n+1}^{(i)}}, \quad \Delta \varphi_{n+1}^{(i+1)} \triangleq \varphi_{n+1} - \varphi_{n+1}^{(i)} \label{eq:derive_of_residual} \end{equation} \] 令 Taylor 展开式等于 0, 就得到关于 \(\Delta \varphi_{n+1}^{(i+1)}\) 的线性方程组: \[ K_{n+1}^{(i)} \Delta \varphi_{n+1}^{(i+1)} = -R_{n+1}(\varphi_{n+1}^{(i)}) \] 求解上述方程组后, 下一迭代步更新为 \[ \varphi_{n+1}^{(i+1)} := \varphi_{n+1}^{(i)} + \Delta \varphi_{n+1}^{(i+1)} \]

Consistent tangent matrix

将残差函数表达式 \(\ref{eq:residual}\) 代入到式 \(\ref{eq:derive_of_residual}\) 中得到: \[ \frac{\partial R_{n+1}}{\partial \varphi_{n+1}} = \frac{\partial f^{\mathrm{int}}(\varphi_{n+1}, t_{n+1})}{\partial \varphi_{n+1}} - \frac{\partial f^{\mathrm{ext}}(\varphi_{n+1}, t_{n+1})}{\partial \varphi_{n+1}} + K^{\mathrm{nmk}} \] 由上式可见, 刚度矩阵由 内力 外力 Newmark 算法三部分构成. 其中 \(K^\mathrm{nmk}\) 由式 \(\ref{eq:stiff_nmk}\) 给出, 接下来将讨论前两项对刚度矩阵的贡献.

内力刚度

St.Venant-Kirchhoff 材料

需要刻画函数变量 \(\varphi\) 的增量对单元内力的影响. 内力虚功可以写作节点虚位移和内力的点积形式: \[ \delta W^{\mathrm{int}} = \int_{S_0} \delta E_{I J} S_{I J} h \mathrm{~d} S_0 = [\delta \varphi]^{\top} [\mathrm{f}^{\mathrm{int}}] = \delta \varphi_{i}^{a} f_i^{\mathrm{int},a} \]

对上式关于 \(\varphi_{j}^{b}\) 求偏导数得到内力的刚度矩阵 \[ \frac{ \partial (\delta \varphi_{i}^{a} f_i^{\mathrm{int},a}) }{ \partial \varphi_j^b } = 0 + \delta \varphi_{i}^{a} K_{ij}^{\mathrm{int},ab} \]

求导与积分交换次序得到, \[ \begin{equation} \delta \varphi_{i}^{a} K_{ij}^{\mathrm{int},ab} = \frac{ \partial ~}{ \partial \varphi_{j}^{b}} \int_{S_0} \delta E_{I J} S_{I J} h \mathrm{~d} X = \int_{S_0} \delta E_{I J} \frac{ \partial S_{I J}}{ \partial \varphi_{j}^{b}} h \mathrm{~d} X + \int_{S_0} \delta \frac{ \partial E_{I J}}{ \partial \varphi_{j}^{b}} S_{I J} h \mathrm{~d} X \label{eq:mem_tg_mat} \end{equation} \]

式中, 第一项是由材料本构引起的刚度项, 代入材料本构 \(\ref{eq:constitutive}\) \[ \int_{S_0} \delta E_{I J} \frac{ \partial S_{I J}}{ \partial \varphi_{j}^{b}} h \mathrm{~d} S_0 = \int_{S_0} \delta E_{I J} L_{IJKL} \frac{ \partial E_{KL}}{ \partial \varphi_{j}^{b}} h \mathrm{~d} S_0 = [\delta \varphi]^{\top} \cdot \int_{S_0} [\mathrm{B}]^{\top} [\mathrm{D}] [\mathrm{B}] h ~\mathrm{d} X \]

式中, \([\mathrm{D}]\) 是四阶张量 \(\mathbb{L}\) 在 Voigt 记法下的矩阵形式. 因此, \[ \boxed{ [\mathrm{K}^{\mathrm{mat}}] = \int_{S_0} [\mathrm{B}]^{\top} [\mathrm{D}] [\mathrm{B}] h ~\mathrm{d} X = h A_0 [\mathrm{B}]^{\top} [\mathrm{D}] [\mathrm{B}] } \] 式 \(\ref{eq:mem_tg_mat}\) 中的第二项的刚度由几何非线性导致, \[ \int_{S_0} \delta \frac{ \partial E_{I J}}{ \partial \varphi_{j}^{b}} S_{I J} h \mathrm{~d} X = \frac{1}{2} \int_{S_0} \delta \frac{ \partial g_{ij}}{ \partial \varphi_{j}^{b}} s_{ij} h \mathrm{~d} X \]

\[ \begin{aligned} & ~\mathrm{d}\left(\delta g_{11}\right)=2\left(\delta \tilde{\boldsymbol{\varphi}}^2-\delta \tilde{\boldsymbol{\varphi}}^1\right)^\top \cdot \left(~\mathrm{d} \tilde{\boldsymbol{\varphi}}^2-~\mathrm{d} \tilde{\boldsymbol{\varphi}}^1\right) \\ & ~\mathrm{d}\left(\delta g_{22}\right)=2\left(\delta \tilde{\boldsymbol{\varphi}}^3-\delta \tilde{\boldsymbol{\varphi}}^1\right)^\top \cdot \left(~\mathrm{d} \tilde{\boldsymbol{\varphi}}^3-~\mathrm{d} \tilde{\boldsymbol{\varphi}}^1\right) \\ & ~\mathrm{d}\left(2\delta g_{12}\right)=2\left(\delta \tilde{\boldsymbol{\varphi}}^2-\delta \tilde{\boldsymbol{\varphi}}^1\right)^\top \cdot \left(~\mathrm{d} \tilde{\boldsymbol{\varphi}}^3-~\mathrm{d} \tilde{\boldsymbol{\varphi}}^1\right) +2\left(\delta \tilde{\boldsymbol{\varphi}}^3-\delta \tilde{\boldsymbol{\varphi}}^1\right)^\top \cdot \left(~\mathrm{d} \tilde{\boldsymbol{\varphi}}^2-~\mathrm{d}\tilde{\boldsymbol{\varphi}}^1\right) \end{aligned} \]

因此 $$ {S_0} ~( g{11} ) s_{11} ~ X = []^{} _{[_{11}]} [~ ] \

{S_0} ~( g{22} ) s_{22} ~ X = []^{} _{[_{22}]} [~ ] \

{S_0} ~( 2g{12} ) s_{12} ~ X = []^{} _{[_{12}]} [~ ] \[ 最终得到 \] $$

Neo-Hookean 材料

式 \(\ref{eq:f_int_tau_D}\) 对节点位移微分, 得到刚度矩阵: \[ ~\mathrm{d} f_{i}^{\mathrm{int},a} = K_{ik}^{\mathrm{int},ab} ~\mathrm{d} \varphi_{k}^{b} \] 微分与积分号交换顺序, 得到 \[ \begin{equation} ~\mathrm{d} f_{i}^{\mathrm{int},a} = \underbrace{\int_{S_0} ~\mathrm{d} \tau_{ij} N_{,I}^a F_{Ij}^{-1} ~\mathrm{d} X}_{\triangleq I_1} + \underbrace{\int_{S_0} \tau_{ij} N_{,I}^a ~\mathrm{d} F_{Ij}^{-1} ~\mathrm{d} X}_{\triangleq I_2} \label{eq:neo_hookean_stiff} \end{equation} \] 等号右边第一项积分中, 对 Kirchhoff 应力的微分为 \[ \mathrm{d} \tau_{ij} = \frac{ \partial \tau_{ij} }{ \partial F_{lL} } \frac{ \partial F_{lL} }{ \partial \varphi_{k}^{b} } ~\mathrm{d} \varphi_{k}^{b} \] 其中, \[ \frac{ \partial F_{lL} }{ \partial \varphi_{k}^{b} } = \frac{ \partial \varphi_l^a }{ \partial \varphi_{k}^{b} } N_{,L}^{a} = \delta_{kl} N_{,L}^{b} \] 代入积分项中得到 \[ I_1 = \int_{S_0} \frac{ \partial \tau_{ij} }{ \partial F_{kL} } N_{,L}^{b} N_{,I}^a F_{Ik}^{-1} ~\mathrm{d} X ~\mathrm{d} \varphi_{k}^{b} \] 如果取形函数的微分为关于空间坐标的微分, 那么有 \[ N_{,j}^a = N_{,I}^a F_{Ij}^{-1} \] 定义 \[ \begin{equation} C_{ijkl} \triangleq \frac{ \partial \tau_{ij} }{ \partial F_{kL} } F_{lL} = J \delta_{kl} \sigma_{ij} + J \frac{ \partial \sigma_{ij} }{ \partial F_{kL} } F_{lL} \label{eq:def_C} \end{equation} \] 上两式代入 \(I_1\) 中得到 \[ \boxed{ I_1 = \int_{S_0} C_{ijkl} N_{,l}^{b} N_{,j}^a ~\mathrm{d} X ~\mathrm{d} \varphi_{k}^{b} } \] 接下来要考虑 \(C_{ijkl}\) 定义式 \(\ref{eq:def_C}\) 中出现的 \({ \partial \sigma_{ij} }/{ \partial F_{kL} }\) 项. Neo-Hookean 材料应力应变关系式为 \[ \sigma_{ij} = \frac{\mu}{J^{5/3}} \left( B_{ij} - \frac{1}{3} B_{kk} \delta_{ij} \right) + K ( J - 1 ) \delta_{ij} \] 根据如下导数公式 \[ \frac{ \partial (F_{iI} F_{jI}) }{ \partial F_{kL} } = \delta_{ik} F_{jL} + \delta_{jk} F_{iL} \\ \frac{ \partial B_{ii} }{ \partial F_{kL} } = 2 F_{kL} \\ \frac{ \partial J }{ \partial F_{kL} } = J F_{Lk}^{-1} \] 就得到 \[ \begin{aligned} \frac{ \partial \sigma_{ij} }{ \partial F_{kL} } &= \mu (-\frac{5}{3}) J^{-5/3} F_{Lk}^{-1} \left( B_{ij} - \frac{1}{3} B_{kk} \delta_{ij} \right) \\ &+ \mu J^{-5/3} \left( \delta_{ik} F_{jL} + \delta_{jk} F_{iL} - \frac{2}{3} F_{kL} \delta_{ij} \right) \\ &+ K J F_{Lk}^{-1} \delta_{ij} \end{aligned} \] 再代入到式 \(\ref{eq:def_C}\) 得到 \[ \begin{equation} \boxed{ C_{ijkl} = \mu J^{-2/3} \left( \delta_{ik} B_{jl} + \delta_{jk} B_{il} - \frac{2}{3} \delta_{ij} B_{kl} - \frac{2}{3} B_{ij} \delta_{kl} + \frac{2}{9} B_{kk} \delta_{ij} \delta_{kl} \right) + K J ( 2J - 1 ) \delta_{ij} \delta_{kl} } \label{eq:fml_C} \end{equation} \] 现在来考察积分式 \(I_2\). \(I_2\) 中出现的 \(~\mathrm{d} F_{Ij}^{-1}\) 的微分为 \[ ~\mathrm{d} \boldsymbol{F}^{-1} = (\boldsymbol{I} + \boldsymbol{F}^{-1} \cdot \mathrm{d}\boldsymbol{F} )^{-1} \cdot \boldsymbol{F}^{-1} -\boldsymbol{F}^{-1} \sim -\boldsymbol{F}^{-1} \cdot \mathrm{d}\boldsymbol{F} \cdot \boldsymbol{F}^{-1}\quad \scriptsize{\left( (I + A)^{-1} \sim I - A \right)} \] 因此有 \[ ~\mathrm{d} F_{Ij}^{-1} = -F_{Ik}^{-1} F_{Jj}^{-1} N_{,J}^b ~\mathrm{d} \varphi_{k}^{b} \] 代入到式 \(\ref{eq:neo_hookean_stiff}\) 中积分项 \(I_2\) 有 \[ \boxed{ I_2 = -\int_{S_0} \tau_{ij} N_{,I}^a N_{,J}^b F_{Ik}^{-1} F_{Jj}^{-1} ~\mathrm{d} X, \quad \text{或者} \quad I_2 = -\int_{S_0} \tau_{ij} N_{,k}^a N_{,j}^b ~\mathrm{d} X } \]

载荷刚度

\[ \begin{align*} ~\mathrm{d} \boldsymbol{f}_{p}^{a} &= \frac{1}{6} p \left( ~\mathrm{d}\tilde{\boldsymbol{x}}^{21} \times \tilde{\boldsymbol{x}}^{31} + \tilde{\boldsymbol{x}}^{21} \times ~\mathrm{d}\tilde{\boldsymbol{x}}^{31} \right) \\ &= \frac{1}{6} p \left( ( \tilde{\boldsymbol{x}}^{31} - \tilde{\boldsymbol{x}}^{21}) \times \mathrm{d}\tilde{\boldsymbol{x}}^{1} - \tilde{\boldsymbol{x}}^{31} \times ~\mathrm{d}\tilde{\boldsymbol{x}}^{2} + \tilde{\boldsymbol{x}}^{21} \times ~\mathrm{d}\tilde{\boldsymbol{x}}^{3} \right) \\ &= \frac{1}{6}p \begin{bmatrix} [\widehat{\tilde{x}^{31}}] - [\widehat{\tilde{x}^{21}}] & - [\widehat{\tilde{x}^{31}}] & [\widehat{\tilde{x}^{21}}] \end{bmatrix}_{3 \times 9} [~\mathrm{d} \varphi] \end{align*} \]

所以载荷刚度矩阵为 \[ \boxed{ \mathrm{K}^{\mathrm{ext}} = \begin{bmatrix} [\widehat{\tilde{x}^{31}}] - [\widehat{\tilde{x}^{21}}] & - [\widehat{\tilde{x}^{31}}] & [\widehat{\tilde{x}^{21}}] \\ [\widehat{\tilde{x}^{31}}] - [\widehat{\tilde{x}^{21}}] & - [\widehat{\tilde{x}^{31}}] & [\widehat{\tilde{x}^{21}}] \\ [\widehat{\tilde{x}^{31}}] - [\widehat{\tilde{x}^{21}}] & - [\widehat{\tilde{x}^{31}}] & [\widehat{\tilde{x}^{21}}] \\ \end{bmatrix} } \]

准稳态问题

膜结构一般在初始阶段没有刚度, 因此即使对于准静态问题, 也应添加人为的阻尼项稳定计算, 在速率项足够小的情况下再切换到稳态计算. 稳态计算需要保留方程 \(\ref{eq:ode}\) 中的阻尼项: \[ [\mathrm{C}] [\dot{\varphi}] + h A_0 [\mathrm{B}]^\top [\mathrm{S}] - [\mathrm{f}] = 0 \] 并使用向后欧拉方法进行求解: \[ \varphi_{n+1} = \varphi_{n} + \Delta t \dot{\varphi}_{n+1} \] 这对应 Newmark 方法中 \(\beta = \gamma\).