应用主从节点法在有限元方程中施加多点约束
有限元离散后, 要求节点处内力与外力平衡, \(\boldsymbol{\mathsf{f}}^{\mathrm{int}}=\boldsymbol{\mathsf{f}}^{\mathrm{ext}}\). 其中, 内力可由刚度矩阵 \(\boldsymbol{\mathsf{K}}\) 和节点自由度列阵 \(\boldsymbol{\mathsf{d}}\) 给出, 最终得到离散后的平衡方程: \[ \begin{equation} \boldsymbol{\mathsf{K}} \boldsymbol{\mathsf{d}} = \boldsymbol{\mathsf{f}}^{\mathrm{ext}}. \end{equation} \] 除平衡方程之外, 节点之间一般还需满足若干运动学约束: \[ F_{k}( \boldsymbol{d}_{a}, \ldots ) = \boldsymbol{0}, \quad k=1,\ldots \] 当上述运动学约束只与单个节点 \(\boldsymbol{d}_{a}\) 的自由度相关时, 此时等价于在该节点处施加位移边界条件. 当运动学约束与多个节点相关时, 此时称该约束为多点约束. 由于这些运动学约束的存在, 刚度矩阵 \(\boldsymbol{\mathsf{K}}\) 一般不是满秩矩阵, 需要在施加约束之后才能唯一定解节点位移. 施加位移边界条件 (单点约束) 可直接消除刚度矩阵对应自由度编号的行和列, 而施加多点约束一般有如下三种方式: (1) 主从节点法 (2) 罚函数法 和 (3) 拉格朗日乘子法.
主从节点法是将多点约束直接代入到线性方程组中, 消除从节点自由度. 虽然该方法思路简单, 但由于
- 确定主从节点对的方案并不唯一, 并且很难找到一个通用的算法, 在任意给定的约束方程组中确定正确约束数量的主从节点对;
- 施加约束需要对刚度矩阵的元素进行调整, 移动或删除操作, 这可能会在刚度矩阵按照稀疏矩阵存储时遇到困难;
- 该方法只能施加线性约束条件,
所以主从节点法在通用性方面远不如 罚函数法 或 拉格朗日乘子法. 不过在针对特定问题时 (例如, 求解在周期性边界条件约束下单胞内的线弹性问题), 多点约束的形式较单一, 所以仍可以找到一套系统的算法, 确定问题中的主从节点对和高效调整刚度矩阵.
主从节点法
主从节点法本质上是代入方法, 对于原始未知量 \(\boldsymbol{\mathsf{d}}\), 维度为
NDOF x 1, 在施加 NcDOF 约束之后, 可通过数量为
NiDOF=NDOF-NcDOF 的独立未知量 \(\hat{\boldsymbol{\mathsf{d}}}\) 表示为
\[
\begin{equation}
\underbrace{\boldsymbol{\mathsf{d}}}_{\scriptsize \mathrm{NDOF}\times1}
= \underbrace{\boldsymbol{\mathsf{T}}}_{\scriptsize
\mathrm{NDOF}\times\mathrm{NiDOF}}\
\underbrace{\hat{\boldsymbol{\mathsf{d}}}}_{\scriptsize
\mathrm{NiDOF}\times1}.
\end{equation}
\] 将该方程代入平衡方程之后, 再左乘 \(\boldsymbol{\mathsf{T}}^{\top}\),
就得到满秩的线性代数方程组: \[
\begin{equation}
\underbrace{\hat{\boldsymbol{\mathsf{K}}}}_{\scriptsize
\mathrm{NiDOF}\times\mathrm{NiDOF}}\
\underbrace{\hat{\boldsymbol{\mathsf{d}}}}_{\scriptsize
\mathrm{NiDOF}\times1} =
\underbrace{\hat{\boldsymbol{\mathsf{f}}}}_{\scriptsize
\mathrm{NiDOF}\times1}, \quad \text{式中}\quad
\hat{\boldsymbol{\mathsf{K}}} = \boldsymbol{\mathsf{T}}^{\top}
\boldsymbol{\mathsf{K}} \boldsymbol{\mathsf{T}}, \quad
\hat{\boldsymbol{\mathsf{f}}} = \boldsymbol{\mathsf{T}}^{\top}
\boldsymbol{\mathsf{f}}.
\end{equation}
\]
约束自由度又可以分成两类: 施加位移边界条件的自由度 和
多点约束中的从节点自由度, 其数量分别为 NdcDOF 和
NscDOF. 因此, 最终方程中独立未知数的个数为 \[
\mathrm{NiDOF} = \mathrm{NDOF} - \mathrm{NdcDOF} - \mathrm{NscDOF}
\]
使用主从节点法的主要难点在于:
- 如何确保约束矩阵 \(\boldsymbol{\mathsf{T}}\) 是满秩矩阵?
- 如何在代码实现中避免构造 \(\boldsymbol{\mathsf{T}}\) 获取修正后的刚度矩阵 \(\hat{\boldsymbol{\mathsf{K}}}\)?
以下将针对线下阶段影响张量的求解, 分别给出对上述两个问题的具体实现.
应用: 线下阶段求解影响张量
在对方形单胞相对面 \(S^{-}\), \(S^{+}\) 施加周期性边界条件时, 方形单胞的相对面节点位移满足如下约束: \[ \boldsymbol{d}_{a}^{+}-\boldsymbol{d}_{a}^{-} = \boldsymbol{d}_{b}^{+} - \boldsymbol{d}_{b}^{-}, \quad \forall ~\text{node}~ a,b \in S \]
单胞问题主从节点对的确定
在求解影响张量时, 施加在单胞上的宏观应变等于零, 这意味着对于方形单胞,
8 个角点的位移全部等于零: \[
\boldsymbol{d}_{a}=0, \quad \text{node}\ a \text{ at corner}
\] 正是因为宏观应变等于零, 所以对任意相对面上的节点, 位移均相等:
\[
\begin{equation}
\boldsymbol{d}_{a}^{+} = \boldsymbol{d}_{a}^{-}, \quad \forall
~\text{node}~ a,b \in S
\end{equation}
\] 因此, 确定约束矩阵 \(\boldsymbol{\mathsf{T}}\),
就转变为寻找所有在 (1) 边界面 和 (2) 棱线 处的节点对 \(\{ a^{+}, a^{-} \}\), \(\{a^{++}, a^{+-}, a^{-+}, a^{--}\}\),
并在节点对中确定一个节点作为主节点, 其它所有节点作为从节点,
将从节点的三个自由度全部设置为属性 slave,
将在之后调整线性方程组时消去该自由度.
[!NOTE]
\(\text{FE}^{2}\) 问题中单胞的宏观应变不等于零, 因此两个节点之间的约束关系并不是齐次的, 在确定约束关系时与此处给出的方法不同.
线性方程组的调整
首先, 对于属于位移边界条件的自由度, 在组装全局刚度矩阵时直接跳过,
由此得到数量为 NDOF-NdcDOF 的线性方程组: \[
\boldsymbol{\mathsf{K}}_{uu} \boldsymbol{\mathsf{d}}_{u} =
\boldsymbol{\mathsf{f}}_{u}
\] 将未知量 \(\boldsymbol{\mathsf{d}}_{u}\) 分解为 \((\hat{\boldsymbol{\mathsf{d}}};
\boldsymbol{\mathsf{d}}_{s})\), 其中 \(\hat{\boldsymbol{\mathsf{d}}}\)
表示独立的自由度, 维度为 NiDOF x 1, \(\boldsymbol{\mathsf{d}}_{s}\)
表示从属自由度, 维度为 NscDOF x 1. 对于式 (4) 中的约束关系,
可以定义维度为 NscDOF x NiDOF 矩阵 \(\boldsymbol{\mathsf{A}}\), 在第
sidx 行处只在列索引为 s2mMap(sidx) 处等于1,
其它位置均等于零, s2mMap
将从属自由度索引映射到对应的主自由度索引. 由此, 式 (4) 约束关系,
以及约束矩阵 \(\boldsymbol{\mathsf{T}}\) 可以表示为 \[
\boldsymbol{\mathsf{d}}_{s} = \boldsymbol{\mathsf{A}}
\hat{\boldsymbol{\mathsf{d}}}, \quad
\boldsymbol{\mathsf{T}} = \begin{bmatrix}\boldsymbol{\mathsf{I}} \\
\boldsymbol{\mathsf{A}} \end{bmatrix}
\] 所以, 添加约束关系后的刚度矩阵和外力列阵等于 \[
\hat{\boldsymbol{\mathsf{K}}} = \begin{bmatrix}\boldsymbol{\mathsf{I}}
& \boldsymbol{\mathsf{A}} ^{\top}\end{bmatrix}
\boldsymbol{\mathsf{K}} \begin{bmatrix}\boldsymbol{\mathsf{I}} \\
\boldsymbol{\mathsf{A}} \end{bmatrix}, \quad
\hat{\boldsymbol{\mathsf{f}}} = \begin{bmatrix}\boldsymbol{\mathsf{I}}
& \boldsymbol{\mathsf{A}} ^{\top}\end{bmatrix}
\boldsymbol{\mathsf{f}}.
\] 此时, \(\boldsymbol{\mathsf{A}}^{\top}\)
的作用相当于将索引为 sidx 的从属自由度组装在
s2mMap(sidx) 位置.