跳转至

外力虚功与总体方程的组装

在内力虚功的离散化中，已将弱形式左端归结为单元内力向量 \(\boldsymbol{q}^e\)（UL 法）或 \(\boldsymbol{Q}^e\)（TL 法）。本章从外力虚功引入单元节点外力向量 \(\boldsymbol{F}^e\)，并经过将单元节点物理量按全局节点编号重新排列和汇总的组装操作，最终得到 FrontISTR 非线性结构分析所需求解的关于节点位移的非线性方程组。

外力虚功的单元分解

虚功原理右端可按单元分解为由体积力（物体力）和力学边界上的规定表面力组成的外力虚功。为以矩阵形式表示形函数与有限元近似中引入的位移插值，使用将节点 \(\alpha\) 的形函数 \(N_\alpha^e\) 排列在对角线上的 \(d \times d\) 块 \(\boldsymbol{N}_\alpha\)，以及将这些块横向排列得到的 \(\boldsymbol{N} = [\boldsymbol{N}_1, \ldots, \boldsymbol{N}_{n_e}]\)，令 \(\delta\boldsymbol{u} = \boldsymbol{N}\, \delta\boldsymbol{u}^e\)。将其代入以参考构形表示的外力虚功，可得

\[ \delta W^{\mathrm{ext}} = \sum_e \delta\boldsymbol{u}^{eT} \boldsymbol{F}^e, \qquad \boldsymbol{F}^e_\alpha = \int_{\Omega^e_0} \boldsymbol{N}_\alpha^T \rho_0 \boldsymbol{g}\, dV + \int_{\Gamma^e_{0t}} \boldsymbol{N}_\alpha^T \bar{\boldsymbol{t}}_0\, d\Gamma_0 \]

由此得到上式。这里，将单元节点外力向量排列为 \(\boldsymbol{F}^e = (\boldsymbol{F}^{eT}_1, \ldots, \boldsymbol{F}^{eT}_{n_e})^T\)。这样，外力虚功就可以与内力一侧相同，整理为“单元节点向量 × 测试函数”的形式（即使用当前构形表示，通过 \(dV \to dv\)、\(\rho_0 \to \rho\)、\(\bar{\boldsymbol{t}}_0 \to \bar{\boldsymbol{t}}\) 的替换也可得到相同形式）。

单元节点物理量的组装

将各单元得到的节点物理量 \(\boldsymbol{Q}^e_\alpha, \boldsymbol{F}^e_\alpha\) 汇总到按全局节点编号排列的全局向量中。将单元 \(\Omega^e\) 的局部节点编号 \(\alpha\) 所对应的全局节点编号表示为

\[ \mathrm{gdx}(e, \alpha) = i_g \]

则单元节点物理量与全局节点物理量的相应分量一致（例如 \(\boldsymbol{u}^e_\alpha = \boldsymbol{u}_{i_g}\)）。一般而言，节点 \(i_g\) 会被多个单元共享，因此定义由全局节点编号为 \(i_g\) 的 \((e, \alpha)\) 组合所构成的集合

\[ \mathcal{E}(i_g) = \{ (e, \alpha) \mid \mathrm{gdx}(e, \alpha) = i_g \} \]

利用该集合将求和改写为 \(\sum_e \sum_\alpha = \sum_{i_g} \sum_{(e,\alpha) \in \mathcal{E}(i_g)}\)，即可得到遍及全部 \(n_g\) 个节点的节点内力以及全局内力向量

\[ \boldsymbol{Q}_{i_g} = \sum_{(e,\alpha) \in \mathcal{E}(i_g)} \boldsymbol{Q}^e_\alpha, \qquad \boldsymbol{Q} = (\boldsymbol{Q}^T_1, \ldots, \boldsymbol{Q}^T_{n_g})^T \]

由此得到上式。\(\boldsymbol{Q}_{i_g}\) 相当于作用于节点 \(i_g\) 的各单元节点内力的合力；若无外力作用且处于平衡状态，则该量为 \(\boldsymbol{0}\)。UL 法也可通过相同步骤得到 \(\boldsymbol{q}_{i_g}, \boldsymbol{q}\)，并且数值上 \(\boldsymbol{q} = \boldsymbol{Q}\)，因此下文除需要区分之处外统一记为 \(\boldsymbol{Q}\)。全局外力向量 \(\boldsymbol{F}\) 也通过相同的汇总得到。

在实现中并不显式构造集合 \(\mathcal{E}(i_g)\)，而是在单元循环中累加到相应分量。

将全局内力向量 Q 初始化为 0：Q_{i_g} = 0  (i_g = 1, ..., n_g)
for e = 1 to (单元数)
    for α = 1 to n_e
        i_g = gdx(e, α)
        Q_{i_g} += Q^e_α
    end for
end for

全局外力向量 \(\boldsymbol{F}\) 也按相同步骤构造。将单元节点物理量加到按全局节点编号编号的向量或矩阵中并存储的操作称为组装（assemble）。对于与两个节点编号相关的二阶张量（如刚度矩阵），使用集合 \(\mathcal{E}^2(i_g, i_h) = \{ (e, \alpha, \beta) \mid \mathrm{gdx}(e, \alpha) = i_g\ \mathrm{and}\ \mathrm{gdx}(e, \beta) = i_h \}\) 可得到同类的组装操作（具体构造参见切线刚度矩阵）。

待求解的非线性方程

将内力和外力的组装结果代入虚功原理，并利用其对满足几何边界条件的任意测试函数 \(\delta\boldsymbol{u}^n\) 均成立这一事实，可得

\[ \boldsymbol{Q}(\boldsymbol{u}^n) - \boldsymbol{F}(\boldsymbol{u}^n) = \boldsymbol{0} \]

由此得到上式。在增量分析（增量分析框架）的语境下恢复时间下标 \(_{n+1}\)，并省略表示全局节点向量的上标 \(^n\) 后，待求方程为

\[ \boldsymbol{Q}(\boldsymbol{u}_{n+1}) - \boldsymbol{F}(\boldsymbol{u}_{n+1}) = \boldsymbol{0} \]

即为上式。求时刻 \(t_{n+1}\) 的节点位移 \(\boldsymbol{u}_{n+1}\) 的离散边值问题，归结为将这个关于位移的非线性方程与几何边界条件联立求解。方程的线性化及切线刚度矩阵的构造参见切线刚度矩阵，迭代求解方法参见 Newton-Raphson 法。

相关内容