邱小伙的博客

Vim 的四种模式

• Normal mode: 普通模式, 移动光标或使用命令修改文本
• Visual mode: 浏览模式, 选择文本
• Insert mode: 插入模式, 编辑文本
• Command mode: 命令模式, 执行命令
  阅读全文 »

极限的定义, 以及否定形式

关于序列 \(\{a_n\}\) 的极限命题为:

对任意的 \(\varepsilon > 0\), 存在 \(N \in \mathbb{N}\), 当 \(n>N\) 时, \(|a_n - a| < \varepsilon\)

我们用含参数的陈述表示上述命题的后半段:

\[ S(N, \varepsilon): n > N \Rightarrow |a_n - a| < \varepsilon \]

定义某一空间 \(V\) 上的范数, 就是确定该空间上某一个非负实值映射 \(\| \cdot \|\), 且该映射满足如下条件:

1. 对任意的 \(v \in V\), \(\| v \| \geq 0\) 成立, 当且仅当 \(v = 0\) 时 \(\| v \| = 0\);
2. 对任意的 \(c \in \mathbb{R}\) 和 \(v \in V\), \(\|cv\| = |c| \|v\|\);
3. 对任意的 \(u, v \in V\), 不等式 \(\|u+v\| \leq \|u\| + \|v\|\) 成立.

我现在想处理的是, 当空间的概念从有限维的向量空间, 拓展到无穷维的函数空间上时, 条件 1 中 \(v = 0\) 的定义也需要加以推广, 以适应这个条件中 \(v=0\) 与 \(\|v\| = 0\) 的等价性.

首先我们先理清一下这三者的关系：

• Hexo 本身是自动生成博客网页的框架代码，它不需要我们自己编写网页（例如 .css，.js ），而是修改 _config.yml 文件提供的参量，自动生成静态网页，可以在本机浏览；
• NexT 是 Hexo 的主题插件，本身依赖于 Hexo，提供额外的主题。它同样也提供 _config.yml 文件供我们配置；
• 使用 Github 将生成的静态网页推送到互联网当中。
  阅读全文 »

一般的有限元平衡方程

平衡态的方程可以表示为: \[ 0 = \mathbf{f}^{\rm int}(\mathbf{d}^{n+1}, t^{n+1}) - \mathbf{f}^{\rm ext}(\mathbf{d}^{n+1}, t^{n+1}) \triangleq \mathbf{r}^{n+1}(\mathbf{d}^{n+1}) \] 式中, \(\mathbf{f}^{\rm int}\) 和 \(\mathbf{f}^{\rm ext}\) 分别表示 \(t^{n+1}\) 时刻的内力与外力. 内力与外力是关于节点位移 \(\mathbf{d}^{n+1}\) 和时刻 \(t^{n+1}\) 的函数. \(\mathbf{r}^{n+1}: \mathbb{R}^n \mapsto \mathbb{R}^n\) 表示在时刻 \(t^{n+1}\) 关于节点 \(\mathbf{d}\) 的非线性函数, 方程 \[ \mathbf{r}^{n+1}(\mathbf{d})=0 \tag{1} \] 的零点即为使得系统在 \(t^{n+1}\) 平衡的节点位移.

壳体插值函数

三维空间中壳单元的插值函数为 \[ \begin{aligned} \boldsymbol{x}(\xi, \eta, \zeta) &= \overline{\boldsymbol{x}}(\xi, \eta) + \boldsymbol{X}(\xi, \eta, \zeta) \\ \overline{\boldsymbol{x}}(\xi, \eta) &= \sum_{a}^{\mathrm{nen}}N^{a}(\xi, \eta) \overline{\boldsymbol{x}}^a \\ \boldsymbol{X}(\xi, \eta, \zeta) &= \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \boldsymbol{X}^a(\zeta) \\ \boldsymbol{X}^a(\zeta) &= z^{a}(\zeta) \hat{\boldsymbol{X}}^a \quad(\text{不求和})\\ z^a(\zeta) &= N_{+}(\zeta) z_{+}^{a} + N_{-}(\zeta) z_{-}^{a} \\ N_{+}(\zeta) &= \frac{1}{2}(1 + \zeta), \quad N_{-}(\zeta) = \frac{1}{2}(1 - \zeta) \end{aligned} \tag{1} \] 式中, \(\overline{\boldsymbol{x}}\) 为壳单元中参考面中的点; \(\boldsymbol{X}\) 是参考平面中指向 “纤维” 方向的向量, \(\hat{\boldsymbol{X}}^a\) 是节点处纤维法向量 (注意, 纤维方向并不一定垂直于参考平面); \(z_a\) 是厚度方程. 壳单元需要输入的参数为 \[ \boxed{ \overline{\boldsymbol{x}}^a, \quad \hat{\boldsymbol{X}}^a, \quad z_{+}^{a}, \quad z_{-}^{a} } \tag{2} \]

• 纤维: 固定坐标 \((\xi, \eta)\) 得到的关于 \(\zeta\) 的参数曲线;
• lamina: 固定坐标 \(\zeta\) 得到的关于 \((\xi, \eta)\) 的参数曲面

壳单元的另一类参数: 如果给出节点纤维的上下端点 \(\boldsymbol{x}_{+}^a\), \(\boldsymbol{x}_{-}^a\), 以及参考面对应的参数 \(\overline{\zeta}\), 那么由此得到式 (2) 所需要的参数为 \[ \overline{\boldsymbol{x}}^a = \frac{1}{2}(1 + \overline{\zeta}) \boldsymbol{x}_{+}^a + \frac{1}{2}(1 - \overline{\zeta}) \boldsymbol{x}_{-}^a \\ \hat{\boldsymbol{X}}^a = \frac{\boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a} {\| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \|} \\ \quad z_{+}^{a} = \frac{1}{2}(1 - \overline{\zeta}) \| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \| \\ \quad z_{-}^{a} = -\frac{1}{2}(1 + \overline{\zeta}) \| \boldsymbol{x}_{+}^a - \boldsymbol{x}_{-}^a \| \]

壳单元坐标系

单元积分点处 Lamina 坐标系

定义在单元积分点上的笛卡尔坐标系 \(\{ \boldsymbol{e}_1^l, ~\boldsymbol{e}_2^l, ~\boldsymbol{e}_3^l \}\), 其中 \(\boldsymbol{e}_3^l\) 垂直于过该点的 Lamina. 该坐标系需要进行构造. 首先根据单元的插值函数得到: \[ \boldsymbol{e}_{\xi} = \frac{\boldsymbol{x}_{,\xi}}{\|\boldsymbol{x}_{,\xi}\|}, \quad \boldsymbol{e}_{\eta} = \frac{\boldsymbol{x}_{,\eta}}{\|\boldsymbol{x}_{,\eta}\|} \] 这两个单位向量不一定相互垂直, 但是切于 Lamina. 因此 \(e_{3}^l\) 可通过这两个向量叉乘得到 \[ \boldsymbol{e}_{3}^{l}\triangleq \frac{\boldsymbol{e}_{\xi} \times \boldsymbol{e}_{\eta}} {\| \boldsymbol{e}_{\xi} \times \boldsymbol{e}_{\eta} \|} \] \(\boldsymbol{e}_\xi\), \(\boldsymbol{e}_\eta\) 的角平分单位向量为 \[ \boldsymbol{e}_{\alpha} = \frac{\boldsymbol{e}_\xi + \boldsymbol{e}_\eta} {\|\boldsymbol{e}_\xi + \boldsymbol{e}_\eta\|}, \quad \boldsymbol{e}_{\beta} = \frac{\boldsymbol{e}_3^l \times \boldsymbol{e}_\alpha} {\|\boldsymbol{e}_3^l \times \boldsymbol{e}_\alpha\|} \] 绕 \(\boldsymbol{e}_3^l\) 顺时针旋转 \(45\deg\) 得到 \[ \boldsymbol{e}_1^l = \frac{\sqrt{2}}{2} (\boldsymbol{e}_{\alpha} - \boldsymbol{e}_{\beta}), \quad \boldsymbol{e}_2^l = \frac{\sqrt{2}}{2} (\boldsymbol{e}_{\alpha} + \boldsymbol{e}_{\beta}) \]

节点处 fiber direction 坐标系

节点处坐标系要求 \(\boldsymbol{e}_3^f\) 与 \(\hat{\boldsymbol{X}}^a\) 重合. 其它两个方向不定, Hughes 给出一个算法, 能够在 \(\hat{\boldsymbol{X}}^a\) 接近 \(e_3\) 时 保持与全局坐标系近乎一致. 不过这类算法, 应该不是唯一的, 只要能避免 \(\hat{\boldsymbol{X}}^a\) 与一个接近平行的向量叉乘即可.

坐标变换

宏观坐标系与上述定义的坐标系之间的坐标变换关系为 \[ \boldsymbol{e}_{i}^{\cdot} = Q_{\alpha i}^{\cdot} \boldsymbol{e}_{\alpha} \] 式中, \(\cdot\) 表示积分点上的 lamina 坐标系或节点处的 fiber 坐标系. 因此, 向量 \(\boldsymbol{e}_{i}^{\cdot}\) 在宏观坐标系中的坐标列向量排布得到旋转矩阵 \([\mathrm{Q}]\): \[ [\mathrm{Q}^{\cdot}] = [\boldsymbol{e}_{1}^{\cdot}, ~\boldsymbol{e}_{2}^{\cdot}, ~\boldsymbol{e}_{3}^{\cdot}] \] 由此得到坐标变换矩阵为 \[ [\mathrm{x}^{\cdot}] = [\mathrm{Q}^{\cdot}]^{\top} [\mathrm{x}] \] lamina 和 fiber 坐标系之间的关系为 \[ \boldsymbol{e}_{i}^{f} = Q_{\alpha i}^{f} \boldsymbol{e}_{\alpha} \\ \boldsymbol{e}_{i}^{l} = Q_{\alpha i}^{l} \boldsymbol{e}_{\alpha} \\ \]

单元的运动学描述

\[ \boldsymbol{u}(\xi, \eta, \zeta) = \overline{\boldsymbol{u}}(\xi, \eta) + \boldsymbol{U}(\xi, \eta, \zeta) \\ \overline{\boldsymbol{u}}(\xi, \eta) = \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \overline{\boldsymbol{u}}^a \\ \boldsymbol{U}(\xi, \eta, \zeta) = \sum_{a}^{\mathrm{nen}} N^{a}(\xi, \eta) \boldsymbol{U}^a(\zeta) \\ \boldsymbol{U}^a(\zeta) = z^{a}(\zeta) \hat{\boldsymbol{U}}^a \quad(\text{不求和}) \]

壳单元假设纤维方向只有旋转没有拉伸, 因此在小转动假设下, 节点处的向量为 \(\hat{\boldsymbol{U}}^a\) \[ \hat{\boldsymbol{U}}^a = \theta_2^a \boldsymbol{e}_{a1}^f - \theta_1^a \boldsymbol{e}_{a2}^f \]

本构方程

在壳体 Lamina 坐标系下的线弹性本构方程为 \[ \sigma^{l} = \mathbb{L}^{l} : \varepsilon^{l} \] 表示成 Voigt 记法后的矩阵形式为 \[ [\sigma^{l}] = [\mathrm{L}^{l}] : [\varepsilon^{l}] \] 壳体的 zero normal stress 假设要求应力分量 \(\sigma_{33}^{l}\) 等于 0, 由此得到关于 \(\varepsilon_{3}^{l}\) 的方程: \[ \varepsilon_{3}^{l} = -\frac{1}{L_{33}^{l}} (L_{3i'}^{l} \varepsilon_{i'}^{l}) \]

本构方程 \[ \sigma_{i'}^{l} = L_{i'j}^{l} \varepsilon_{j}^{l} = L_{i'j'}^{l} \varepsilon_{j'}^{l} + L_{i'3}^{l} \varepsilon_{3}^{l} = \left( L_{i'j'}^{l} - \frac{1}{L_{33}}^{l} L_{i'3}L_{3j'}^{l} \right) \varepsilon_{j'}^{l} \] 得到各向同性线弹性本构矩阵 \[ [\mathrm{L}^{l}] = \frac{E}{(1-\nu^2)} \begin{bmatrix} 1 & \nu & & & \\ & 1 & & & \\ && \dfrac{1-\nu}{2} && \\ &&& \kappa\dfrac{1-\nu}{2} & \\ &&&& \kappa \dfrac{1-\nu}{2} \end{bmatrix} \] 再加上 bending 修正系数 \(\kappa = 5/6\)

B 矩阵的构造

节点的自由度为 3 个位移自由度和 2 个转角自由度. 由此得到 B 矩阵的列式为 \[ [\varepsilon^{l}] = \sum_{a}^{\mathrm{nen}} [\mathrm{B}^a] \begin{bmatrix} \boldsymbol{u}^a \\ \theta_{a1} \\ \theta_{a2} \end{bmatrix} \] 式中, 矩阵 \([\mathrm{B}^a]\) 定义为 \[ [\mathrm{B}^a] = \left[\begin{array}{ll} b_I^u & b_I^\theta \end{array}\right]=\left[\begin{array}{cc} b_1^u & b_1^\theta \\ b_2^u & b_2^\theta \\ b_3^u & b_3^\theta \\ \hdashline b_4^u & b_4^\theta \\ b_5^u & b_5^\theta \end{array}\right] \] 式中的 \(b_{I}^u\) 和 \(b_{I}^{\theta}\) 是行向量, \[ \begin{aligned} & b_I^u=\left\langle b_{I m}^u\right\rangle=\left\langle\begin{array}{lll} b_{I1}^u & b_{I2}^u & b_{I3}^u \end{array}\right\rangle \\ & b_I^\theta=\left\langle b_{I a}^\theta\right\rangle=\left\langle\begin{array}{lll} b_{I 1}^\theta & b_{I2}^\theta \end{array}\right\rangle \end{aligned} \] 注意到, 节点上的位移自由度是在全局坐标系下的分量, 应变是 Lamina 坐标系下的表达式, 因此在求位移导数时需要进行坐标变换. 根据壳单元的运动学方程, 得到 Lamina 坐标系下的位移导数为 \[ \frac{ \partial u_{i}^{l} }{ \partial x_{j}^{l}} = Q_{im}^{l} \frac{ \partial u_{m} }{ \partial x_{j}^{l}} = Q_{im}^{l} \sum_{a}^{\mathrm{nen}}\left( \frac{ \partial N^{a} }{ \partial x_{j}^{l}} \overline{u}_m^a + \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}} (\theta_2^a \boldsymbol{e}_{a1}^f - \theta_1^a \boldsymbol{e}_{a2}^f) \cdot \boldsymbol{e}_{m} \right)\\ = Q_{im}^{l} \sum_{a}^{\mathrm{nen}}\left( \frac{ \partial N^{a} }{ \partial x_{j}^{l}} \overline{u}_m^a - \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}} \boldsymbol{e}_{am2}^f ~\theta_1^a + \frac{ \partial (N^{a} z^a)}{ \partial x_{j}^{l}} \boldsymbol{e}_{am1}^f ~\theta_2^a \right) \]

式中出现的关于 lamina 坐标的导数, 可以通过坐标变换转换为宏观坐标导数的计算: \[ \frac{ \partial ~}{ \partial x_{i}^{l}} = Q_{im} \frac{ \partial ~}{ \partial x_{m}} \]

刚度矩阵

\[ \underbrace{[\mathrm{k}]}_{\rm 5 nen \times 5 nen} = [[\mathrm{k}^{ab}]], \\ \underbrace{[\mathrm{k}^{ab}]}_{\rm 5 \times 5} = \int_{\square} \int_{-1}^{+1} [\mathrm{B}^a]^{\top} [\mathrm{L}^{l}] [\mathrm{B}] J ~\mathrm{d} \zeta ~\mathrm{d} \square \]

外力向量

体积力

\[ \underbrace{[f^{\mathrm{body}}]}_{\rm 5 nen \times 1} = [[f_a^{\mathrm{body}}]] \\ \underbrace{[f_a^{\mathrm{body}}]}_{\rm 5 \times 1} = \int_{\square} \int_{-1}^{+1} [\mathrm{N}^a]^{\top} [\mathrm{b}] J ~\mathrm{d} \zeta ~\mathrm{d} \square \]

式中, 矩阵 \([\mathrm{N}^a]\) 为 \[ [\mathrm{N}^a] = \left[\begin{array}{ccc:cc} N^a & 0 & 0 & -N^a z^a e_{a 12}^f & N^a z^a e_{a 11}^f \\ 0 & N^a & 0 & -N^a z^a e_{a 22}^f & N^a z^a e_{a 21}^f \\ 0 & 0 & N^a & -N^a z^a e_{a 32}^f & N^a z^a e_{a 31}^f \end{array}\right] \]

面力

\[ \underbrace{[f^{\mathrm{surf}}]}_{\rm 5 nen \times 1} = [[f_a^{\mathrm{surf}}]] \\ \underbrace{[f_a^{\mathrm{surf}}]}_{\rm 5 \times 1} = \int_{\square} [\mathrm{N}^a]^{\top} [\mathrm{t}] J_s ~\mathrm{d} \square, \quad \zeta = \begin{cases} +1 & \text{top surface} \\ -1 & \text{bottom surface} \end{cases} \]

式中, \(J_s = \| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta} \|\).

当面力为压力时, 对应的面力为 \[ t = -\zeta p \boldsymbol{n}, \quad \boldsymbol{n} = \frac{\boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta}} {\| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\eta} \|} \] 当面力为壳体面内的剪切力时, 假设已知沿参数坐标线的面力分量, 由此得到 \[ t = h_{\xi} \boldsymbol{e}_{\xi} + h_{\eta} \boldsymbol{e}_{\eta} \]

边力

\[ [f^{\mathrm{edge}}] = [[f_a^{\mathrm{edge}}]] \\ [f_a^{\mathrm{edge}}] = \int_{-1}^{+1} \int_{-1}^{+1} [\mathrm{N}^a]^{\top} [\mathrm{t}] J_e ~\mathrm{d} \xi ~\mathrm{d} \zeta, \quad \eta = \{+1, -1\} \]

式中, \(J_e = \| \boldsymbol{x}_{,\xi} \times \boldsymbol{x}_{,\zeta} \|\).

应力结果

弯矩, 膜力以及横向剪力可以通过对应力分量 在 \(\zeta\) 方向积分得到

完全拉格朗日形式的弱形式可以解释为虚功原理, 也即内力虚功等于惯性力和外力虚功: \[ \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\).

定义的标量值

在实现中, \(\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'} \]

Notation

在确定单元的几何类型, 并已知单元的节点坐标之后, 理应可以得到关于该单元的任何几何信息, 比如单元内部积分点的全局坐标, 形函数的取值, 形函数偏导数的取值等等. 这些在转化为代码中实际的矩阵计算时, 需要约定矩阵的行列分别对应的指标类型. 接下来将建立公式与实际计算中矩阵列式之间一一对应的关系.

在进行公式推导和变量定义之前, 有必要首先声明文中用到的 Notation:

• 使用 \(i,j,k,\ldots\) 表示空间和参数指标, 一般放在下标位置
• 使用 \(a, b,c,\ldots\) 表示单元节点的指标, 一般放在上标位置
• 如无特别说明, 均默认指标使用 Einstein 求和约定
• 公式字体 \(varible\) 表示一般的公式推导; 字体 \(\mathrm{matrix}\) 表示对应变量的矩阵形式, 并一般使用方括号 \([\cdot]\) 包围. \([\cdot]_{ij}\) 表示矩阵分量
• 按照 MATLAB 惯例, \([1 ~ 2 ~ 3]\) 或 \([1,2,3]\) 表示行向量; \([1;2;3]\) 表示列向量

定义的标量值

定义的向量与矩阵值

参数坐标得到全局坐标

给定单元的某一参数坐标 \(\xi \in \mathbb{R}^{\rm nedim}\) 后, 该点的全局坐标 \(x \in \mathbb{R}^{\rm nsdim}\) 为: \[ x(\xi) = x^a N^a(\xi) \tag{1} \] 其变量关系为 \[ \boxed{ \mathrm{ [x] = [N] [encrd]} } \]

Jacobian

Jacobian 可以看作是参数坐标系中的微元 \(\mathrm{d} \xi\) 到全局坐标系中的微元 \(\mathrm{d} x\) 的仿射变换 \[ \mathrm{d} x = J \cdot \mathrm{d} \xi \] 写成分量形式为 \[ \mathrm{d} x_i = \frac{\partial x_i}{\partial \xi_j} \mathrm{d} \xi_j \] 记 Jacobian 的矩阵形式为 \[ [\mathrm{J}]_{ij} = \frac{\partial x_i}{\partial \xi_j} \] 当使用式 (1) 表示的映射关系时, Jacobian 矩阵可以表示为形函数关于参数坐标导数的线性组合: \[ [\mathrm{J}]_{ij} = x_i^a \frac{\partial N^{a}}{\partial \xi_j} \] 将形函数关于参数坐标的导数组成的矩阵记作变量 dNdXi, Jacobian 矩阵就可以表示为 \[ \boxed{ [\mathrm{J}] = \rm [encrd]^{\top} [dNdxi] } \tag{2} \]

\(\partial N / \partial x\) 矩阵

对于某一点处形函数对 \(x\) 的偏导数, 可以通过链式法则获得: \[ \frac{\partial N^{a}}{\partial x_i} = \frac{\partial N^{a}}{\partial \xi_j} \frac{\partial \xi_j}{\partial x_i} \] 上式中, 当单元类型确定之后, 形函数关于参数坐标 \(\xi\) 的导数就是已知并有解析表达式的, 但是参数坐标关于全局坐标的函数表示, 并不总是那么方便地得到. 但庆幸地是, 我们只需要知道某一点上的取值, 而这可以通过反函数定理很快地得到: \[ \left.\frac{\partial x_i}{\partial \xi_k}\right|_{\xi} \left.\frac{\partial \xi_k}{\partial x_j}\right|_{x} = \delta_{ij} \] 当 nsdim = nedim 时, 可以通过 Jacobian 逆矩阵得到 \[ \boxed{ \mathrm{[dNdx] = [dNdxi]} [J]^{-1} } \] 当 nsdim > nedim 时 (空间维度总是大于或等于单元维度), 建立单元坐标系 \(\{e_1, e_2, e_3\}\) 为

1. 二维单元: \(e_3\) 垂直于单元所在平面, \(e_1\) 平行于边 \(1-2\) 并指向节点 \(2\), \(e_2 = e_3 \times e_1\) ;
2. 一维单元: \(e_1\) 从节点 \(1\) 指向节点 \(2\), 其它基矢量与 \(e_1\) 满足右手坐标系关系.

如果记 \(\{e_1, e_2, e_3\}\) 在旧坐标系下的坐标矩阵为 \[ [Q] \triangleq [e_1 ~ e_2 ~ e_3] \] 那么, 新旧坐标之间的变换关系为 \[ [x] = [x'] [Q]^{\top} , \quad [x'] = [x] [Q] \] 节点在新坐标系下坐标为 \(x'^{a}\), 并且单元内所有点的坐标分量 \(x'_i\) ( 二维单元 \(i=3\), 一维单元 \(i=2,3\) ) 为常数. 如果记变换后的单元的节点坐标变量为 encrd_ecsys, 那么与变量 encrd 之间的关系为 \[ \boxed{ \rm [encrd\_ecsys] = [encrd] [Q] } \tag{3} \] 在新坐标系下, 用形函数表示为参数坐标的映射关系为 \[ x' = x'^a N^a(\xi) \] 式中的形函数形式上与未变换之前的形函数是一致的, 因此微元之间的关系同样为 \[ \mathrm{d} x' = J' \cdot \mathrm{d} \xi \] 式中, \(J'\) 为新坐标系下的 Jacobian, 类似于式 (2) \[ [\mathrm{J}'] = \rm [encrd\_ecsys]^{\top} [dNdxi] \] 再根据坐标变换关系 (3), 就得到不同坐标系下 Jacobian 之间的关系 \[ \boxed{ [\mathrm{J}'] = \rm [Q]^{\top} [encrd]^{\top} [dNdxi] = [Q]^{\top} [J] } \]

矩阵 \([J']\) 大于 nedim 的行元素全部等于 0. 去除零行之后得到 \([J'']\): nedim x nedim, 得到在新坐标系下的形函数偏导数为: \[ \mathrm{[dNdx'] = [dNdxi]} [J'']^{-1} \] 如果希望得到关于旧坐标系的偏导数值, 那么根据坐标变换关系: \[ \frac{\partial ~~}{\partial x_{\beta}} = \frac{\partial ~~}{\partial x'_{i}} \frac{\partial x'_{i}}{\partial x_{\beta}} = Q_{\beta i} \frac{\partial ~~}{\partial x'_{i}} \] 需要强调的是上式中指标的范围. \(\beta\) 是旧坐标系下的指标, 因此范围为 \(1,2,3\). \(i\) 是新坐标系下的指标, 在新坐标系下, \(i\) 的范围是 \(1,2\) (二维单元) 或 \(1\) (一维单元). 写成矩阵的形式为 \[ \boxed{ \rm [dNdx] = [dNdx'] [Q(0,0,nsdim,nedim)]^{\top} } \]

六节点单元, 形函数没有问题, 积分点上每次计算得到的应变增量没有问题, 应力计算应该也是没有问题的, 但是在线性第二次迭代时, 得不到正确结果? 这是为什么?

问题应该出现在 src\PField\FElement\FSolid2D\FElemSld2D.cpp 中的迭代过程. 需要再想想

应该不是迭代方法的问题, 还是第一次迭代步中最后的位移计算出错了

B 矩阵

\(B\) 矩阵的形式与应变向量的排布相关. 在小位移假设下, 应变总是关于位移对 \(x\) 偏导数的线性组合, 也即关于 \(\partial N / \partial x\) 的线性组合: \[ \varepsilon_{ij} = \frac{1}{2}(u_{i,j} + u_{j,i}) \] 对位移使用形函数进行插值, 得到 \[ u_{i}(x) = u_{i}^{a} N^{a}(x) \tag{4} \] 式中, \(u_{i}^{a}\) 为节点上的位移分量, 为待求解的未知量. 将上式代入到小应变的表达式当中得到 \[ \varepsilon_{ij} = \frac{1}{2}(u_{i}^{a} \frac{\partial N^{a}}{\partial x_j} + u_{j}^{a} \frac{\partial N^{a}}{\partial x_i}) \tag{5} \] 将上式写成矩阵形式, 得到 \[ \boxed{ \rm [epsilon] = [B] [d] } \] 接下来我们将阐述上式中变量具体的排布方式. 对于变量 \(\rm [epsilon]\), 我们使用 Voigt 记法, 将应变存储为向量形式: \[ [\mathrm{epsilon}] = [\varepsilon_{11}; ~ \varepsilon_{22}; ~ 2\varepsilon_{12}] \quad \text{二维情景} \\ [\mathrm{epsilon}] = [\varepsilon_{11}; ~ \varepsilon_{22}; ~ \varepsilon_{33}; ~ 2\varepsilon_{23}; ~ 2\varepsilon_{13}; ~ 2\varepsilon_{12}] \quad \text{三维情景} \] 对于变量 \([d]\), 设节点 \(a\) 上自由度变量排布为 \[ [\mathrm{u}^a] = [u_1^a ; ~ u_2^a ; ~ \cdots u_{\mathrm{nndof}}^a] \] 单元自由度向量通过 \([u^{a}]\) 进行组装 \[ \rm [d] = [u^1 ; u^2 ; \cdots ; u^{nen}] \] 在定义单元自由度向量和应变向量之后, 式 (5) 中将之联系起来的矩阵 \(B\) 就唯一确定了. 我们首先使用定义的单元自由度向量表示式 (4) 中的 \(u(x)\) \[ [\mathrm{u}](x) = N^{a}(x) [\mathrm{u}^a] = [N^{1}(x)[\mathrm{I}] , N^2(x)[\mathrm{I}], \ldots , N^{\mathrm{nen}}(x)[\mathrm{I}]] [\mathrm{d}] \] 式中, \([\mathrm{I}]\) 为维度 nndof 的单位阵. 定义算子矩阵 \[ [\mathrm{D}] = \begin{bmatrix} \frac{\partial }{\partial x_1} & \\ & \frac{\partial }{\partial x_2} \\ \frac{\partial }{\partial x_2} & \frac{\partial }{\partial x_1} \end{bmatrix} (\text{二维情景}), \quad \begin{bmatrix} \frac{\partial }{\partial x_1} & & \\ & \frac{\partial }{\partial x_2} & \\ & & \frac{\partial }{\partial x_3} \\ & \frac{\partial }{\partial x_3} & \frac{\partial }{\partial x_2} \\ \frac{\partial }{\partial x_3} & & \frac{\partial }{\partial x_1} \\ \frac{\partial }{\partial x_2} & \frac{\partial }{\partial x_1} & \end{bmatrix} (\text{三维情景}) \] 由此应变可以表示为 \[ [\mathrm{epsilon}] = [\mathrm{D}] [\mathrm{u}](x) = [\mathrm{D}] [N^{1}(x)[\mathrm{I}] , N^2(x)[\mathrm{I}], \ldots , N^{\mathrm{nen}}(x)[\mathrm{I}]] [\mathrm{d}] \] 如果定义 \[ [\mathrm{D}] N^{a}(x) = [\mathrm{B}^{a}] \] 那么 \[ [\mathrm{B}] = [\mathrm{D}] [N^{1}(x)[I] , N^2(x)[I], \ldots , N^{\mathrm{nen}}(x)[I]] = [\mathrm{B}^1, \mathrm{B}^2, \ldots, \mathrm{B}^{\mathrm{nen}}] \]

使用下标 \(i,j,k\) 表示空间维度的指标; 使用上标 \(a\) 表示单元节点的指标; 使用上标 \(i\) 表示单元积分点的指标

三角形单元

高阶三角形单元的形函数

一阶的三角形单元形函数为 \[ N^1 = L_1, \quad N^2 = L_2, \quad N^3 = L_3 \] 三角形单元随着阶次增加, 其所覆盖的完备的多项式阶数也恰好逐次增加, 不会像四边形单元, 增加节点数也会增加不完备的多项式项. 更高阶的三角形单元形函数可以通过 Lagrange 函数进行构造, 其构造方法为 \[ N^{a} = l_{I}^{I}(L_1) l_{J}^{J}(L_2) l_{K}^{K}(L_3) \] 式中, \(I+J+K = M\), \(M\) 为三角形单元的阶数. \(l_{k}^{n}(\xi)\) 为 Lagrange 函数, 定义为 \[ l_{k}^{n}(\xi) \triangleq \frac{(\xi - \xi_0)(\xi - \xi_1) \cdots (\xi - \xi_{k-1})(\xi - \xi_{k+1}) \cdots (\xi - \xi_n)} {(\xi_k - \xi_0)(\xi_k - \xi_1) \cdots (\xi_k - \xi_{k-1})(\xi_k - \xi_{k+1}) \cdots (\xi_k - \xi_n)} \] 如此构造的 \(N^a\) 是满足性质 \(N^a(L^{b}) = \delta_{ab}\) 的. 对于节点 \(a\), \[ l_{I}^{I}(L_1^a)=1,\quad l_{J}^{J}(L_2^a)=1,\quad l_{K}^{K}(L_3^a)=1 \] 对于其它节点 \(a'\), 其索引为 \((I', J', K')\), 总存在一个索引分量 (不妨选取 \(I'\)), 使得 \(I' < I\), 这样就得到 \[ l_{I}^{I}(L_1^{a'}) = 0 \Rightarrow N^a(L^{a'}) = 0 \] 对于二阶单元, 备选的 Lagrange 函数为 \[ l_0^0(\xi) = 1, \quad l_1^1(\xi) = 2\xi, ~\{0,~\frac{1}{2}\}, \quad l_2^2(\xi) = \xi(2\xi-1), ~\{0,~\frac{1}{2},~1\} \] 花括号中是 Lagrange 函数的插值点. 因此, 二阶三角形单元的形函数为

角点处的节点: \[ N^a = L_a (2L_a - 1), \quad a=1,2,3 \] 边上的节点: \[ N^4 = 4 L_1 L_2, \quad N^5 = 4 L_2 L_3, \quad N^S6 = 4 L_3 L_1 \]

三角形单元内的坐标可以使用面积坐标进行表示 \[ x(L_1, L_2, L_3) = x_1 L_1 + x_2 L_2 + x_3 L_3 \\ y(L_1, L_2, L_3) = y_1 L_1 + y_2 L_2 + y_3 L_3 \] 式中, \((x_i, y_i)\) 是三角形节点在笛卡尔全局坐标系下的坐标, 面积坐标 \(L_i\) 在 \([0,1]\) 内变化, 并满足如下关系: \[ L_{i}(x^a, y^a) = \delta_{ia}, \quad \sum_{i} L_i = 1 \] 在有限元中, 需要计算单元 \(\Omega^e\) 上形如 \[ I = \iint_{\Omega^e} f(x,y) ~\mathrm{d}x \mathrm{d}y \] 的积分, 而面积坐标提供了一种坐标变换的方式, 使得在任意形状三角形区域的积分可以转化到标准区域进行. 这对数值积分非常友好, 因为 Gauss 型数值积分的精度通常是在标准区域上进行衡量的. 进行积分变换, 第一步是要计算 Jacobian 矩阵, 衡量微元之间的仿射变换关系: \[ \begin{pmatrix} \mathrm{d}x \\ \mathrm{d}y \end{pmatrix} = \frac{\partial (x,y)}{\partial(L_1,L_2)} \begin{pmatrix} \mathrm{d}L_1 \\ \mathrm{d}L_2 \end{pmatrix} \] 对于实体单元, 微元之间的系数为 Jacobian 矩阵的行列式值: \[ \mathrm{d}x \mathrm{d}y = \det \frac{\partial (x,y)}{\partial(L_1,L_2)} \mathrm{d}L_1 \mathrm{d}L_2 \] 对于三节点三角形单元. 行列式值为常值, 等于 2 倍的单元面积: \[ J \triangleq \det \frac{\partial (x,y)}{\partial(L_1,L_2)} = 2S_{\Omega^e} \] 经过坐标变换, 就得到原积分 \(I\) 为 \[ I = J \iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2 \tag{1} \] 式中, \(\triangle\) 表示标准单元区域, 在文中设置的面积坐标下, \(\triangle\) 为 \((0,0)\) \((1,0)\) \((0,1)\) 三节点连线围成的直角三角形区域. 在标准区域 \(\triangle\) 中, 考虑如下积分:

\[ \iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2 \] 数值积分将积分近似为单元区域某些点的取值和权重的乘积, 并进行求和: \[ \iint_{\triangle} \tilde f(L_1, L_2) ~\mathrm{d}L_1 \mathrm{d}L_2 \approx \sum_i \tilde f(\xi^i) w^i \] 为了方便验证求积公式的代数精度, 给出多项式函数在标准区域的积分解析公式: \[ \iint_{\triangle} L_1^a L_2^b L_3^c ~\mathrm{d}L_1 \mathrm{d}L_2 = \frac{a!b!c!}{(a+b+c+2)!} \] 一种常用的三角形区域的数值积分方式为 Hammer 积分, 这里列举选用一个积分点和三个积分点的情景. 只选取一个积分点时 \[ \xi^{1} = (\frac{1}{3}, \frac{1}{3}), \quad w^{1} = \frac{1}{2} \] 可以验证, 这种求积公式是满足一阶代数精度的. 选取三个积分点时, 有两种可能: 第一种求积公式为 \[ \xi^{1} = (0, \frac{1}{2}), \quad \xi^{2} = (\frac{1}{2}, 0), \quad \xi^{3} = (\frac{1}{2}, \frac{1}{2}) \\ w^{1} = w^{2} = w^{3} = \frac{1}{6} \]

第二种求积公式为

现在将标准区域的求积公式代入到原积分公式 (1) 当中得到 \[ I \approx J \sum_i \tilde f(\xi^i) w^i \] 令 \(\widetilde{w}^i \triangleq J w^i\), 就得到 \[ I \approx \sum_i \tilde f(\xi^i) \widetilde{w}^i \]

[!NOTE]

一些教材 (Zienkiewicz) 在处理三角形单元的数值积分时, 将 Jacobian 行列式写成 2 倍的三角形面积, 并将系数 2 并入权重当中: \[ I \approx S_{\Omega^e} \sum_i \tilde f(\xi^i) \hat{w}^i, \quad \hat{w}^i = 2 w^i \] 由此得到的一阶 Hammer 积分: \[ \xi^{1} = (\frac{1}{3}, \frac{1}{3}), \quad w^{1} = 1 \] 和二阶 Hammer 积分 \[ \xi^{1} = (0, \frac{1}{2}), \quad \xi^{2} = (\frac{1}{2}, 0), \quad \xi^{3} = (\frac{1}{2}, \frac{1}{2}) \\ w^{1} = w^{2} = w^{3} = \frac{1}{3} \]

四边形单元

对于四边形四节点单元, 使用的形函数为 \[ N^{I}(\xi, \eta) = \frac{1}{4} (1 + \xi^I \xi) (1 + \eta^I \eta) \] 式中, \((\xi^I, \eta^I)\) 为节点 \(I\) 的参数坐标. 等参变换为 \[ x(\xi, \eta) = \sum_{I=1}^{4} x^I N^{I}(\xi, \eta), \quad y(\xi, \eta) = \sum_{I=1}^{4} y^I N^{I}(\xi, \eta) \] 将积分转换到参数坐标区域中, 对应的 Jacobian 为 \[ J(\xi, \eta) = \frac{1}{4} \begin{bmatrix} \sum\limits_{I=1}^{4} x^I \xi^I (1 + \eta^I \eta) & \sum\limits_{I=1}^{4} x^I \eta^I (1 + \xi^I \xi) \\ \sum\limits_{I=1}^{4} y^I \xi^I (1 + \eta^I \eta) & \sum\limits_{I=1}^{4} y^I \eta^I (1 + \xi^I \xi) \end{bmatrix} \] 因此, 对于原积分为 \[ I = \iint_{\Omega^e} f(x,y) ~\mathrm{d}x \mathrm{d}y \] 将积分区域转换到参数坐标区域, 得到 \[ I = \int_{ -1}^{1}\int_{-1}^{ 1} \tilde f(\xi, \eta) J(\xi, \eta) ~\mathrm{d}\xi \mathrm{d}\eta \] 对上述乘积区域的积分 (\((-1,1) \times (-1,1)\)), 可以使用 Gauss 积分方式进行计算. 首先, 一维情景下, 在区间 \((-1,1)\) 上的 Gauss 积分列式为 \[ \int_{-1}^{1} f(\xi) ~\mathrm{d}\xi \approx \sum_{i=1}^{N_{\rm ipt}} w^{i} f(\xi^{i}) \] 式中, \(N_{\rm ipt}\) 是积分点的数量. 常用的 Gauss 积分点和权重如下表所示

在给出一维情景的 Gauss 积分列式之后, 二维或更高维的情景可以通过一维的例子方便地得到. 对于维度为 \(N_{\rm d}\) 的情景, 在乘积空间 \(\Pi(-1,1)\) 中, Gauss 积分的列式为 \[ \int_{\scriptsize \Pi(-1,1)} f(\xi_1, \xi_{2}, \ldots, \xi_{N^{\rm d}}) \mathrm{d} \xi \approx \underbrace{\sum_{i=1}^{N_{\rm ipt}} \sum_{j=1}^{N_{\rm ipt}} \cdots \sum_{k=1}^{N_{\rm ipt}}}_{N^{\rm d}} f(\xi_1^i, \xi_{2}^j, \ldots, \xi_{N^{\rm d}}^k) \underbrace{w^{i} w^{j} \cdots w^{k}}_{N^{\rm d}} \] 例如, 二维情景的 Gauss 积分为 \[ I = \int_{ -1}^{1}\int_{-1}^{ 1} \tilde f(\xi, \eta) J(\xi, \eta) ~\mathrm{d}\xi \mathrm{d}\eta \approx \sum_{i=1}^{N_{\rm ipt}} \sum_{j=1}^{N_{\rm ipt}} \tilde f(\xi^{i}, \eta^{j}) J(\xi^{i}, \eta^{j}) w^i w^{j} \] 当 \(N_{\rm ipt} = 1\) 时, 积分公式为 \[ I \approx 4\tilde f(0, 0) J(0, 0) \]

[!NOTE]

四边形单元 Jacobian 在参数坐标点 \((0,0)\) 取值为 \[ J(0,0) = \frac{1}{4}S_{\Omega^e} \]

参考文献:

M.A.Crisfield, Non-linear Finite Element Analysis of Solids and Structures.

算例简述

本算例是希望通过一个简单的一维系统，给出非线性计算过程中增量步的概念，这是线性与非线性有限元计算之间一个很大的区别。

我们准备分析的系统如图所示