基函数
设有限元解空间 \(\mathcal{U}_{h}\) 的基函数为
\[
\mathbf{v}_{1},\mathbf{v}_{2},\cdots,\mathbf{v}_{N},
\]
\(\mathbf{v}_{i} = \left[v_{x,i},v_{y,i},v_{z,i}\right]\),其中,\(v_{x,i},v_{y,i},v_{z,i}\) 表示位移在 \(x,y,z\) 方向的插值基函数,通常 \(v_{x,i}=v_{y,i}=v_{z,i}=:v_{i}\),于是位移分布函数为
(3)\[
\mathbf{u} = \sum_{i=1}^{N}\mathbf{u}_{i}v_{i}=\sum_{i=1}^{N}\left(u_{x,i}\mathbf{v}_{x,i}+u_{y,i}\mathbf{v}_{y,i}+u_{z,i}\mathbf{v}_{z,i}\right),
\]
其中,\(\mathbf{u}_{i} = \left[u_{x,i},u_{y,i},u_{z,i}\right]^{T}\),且
\[\begin{split}
\mathbf{v}_{x,i}=\begin{bmatrix}
v_{i} \\ 0 \\ 0
\end{bmatrix},\quad
\mathbf{v}_{y,i}=\begin{bmatrix}
0 \\ v_{i} \\ 0
\end{bmatrix},\quad
\mathbf{v}_{z,i}=\begin{bmatrix}
0 \\ 0 \\ v_{i}
\end{bmatrix},
\end{split}\]
此处,测试函数空间 \(\mathcal{V}_{h} = \mathcal{U}_{h}\)
于是,共计 \(3N\) 个位移自由度
\[
u_{x,1},u_{y,1},u_{z,1},u_{x,2},u_{y,2},u_{z,2},\cdots,u_{x,N},u_{y,N},u_{z,N},
\]
和 \(3N\) 个求解方程
(4)\[\begin{split}
\begin{equation}
\begin{aligned}
&\int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{u}) \, \mathrm{d}V + \int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S \\
=& \int_{\Omega}\mathbf{f}\cdot \mathbf{v}_{*,i}\,\mathrm{d}V + \int_{\Gamma_{N}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S + \int_{\Gamma_{R}} \mathbf{f}_{R} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S,\quad *=x,y,z;\ i=1:N.
\end{aligned}
\end{equation}
\end{split}\]
或 Voigt 形式
(5)\[\begin{split}
\begin{equation}
\begin{aligned}
&\int_{\Omega} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{u}) \, \mathrm{d}V + \int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S \\
=& \int_{\Omega}\mathbf{f}\cdot \mathbf{v}_{*,i}\,\mathrm{d}V + \int_{\Gamma_{N}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i}\ \mathrm{d}S + \int_{\Gamma_{R}} \mathbf{f}_{R} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S,\quad *=x,y,z;\ i=1:N.
\end{aligned}
\end{equation}
\end{split}\]
式 (4) 和 (5) 是关于 \(u_{*,i}\) 的线性方程组,关键在于计算出求解矩阵,即 刚度矩阵,其主要贡献项是
\[
\int_{\Omega} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{u}) \, \mathrm{d}V\quad \text{或} \quad \int_{\Omega} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{u}) \, \mathrm{d}V,
\]
这些积分项在整个计算区域上进行的,将其转换为各个单元上的积分之和
\[
\sum_{E}\int_{E} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{u}) \, \mathrm{d}V\quad \text{或} \quad \sum_{E}\int_{E} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{u}) \, \mathrm{d}V,
\]
根据 (3),\(\mathbf{u}\) 可以表示为 \(\mathbf{v}_{*,i}\) 的线性组合,且起系数为待求解自由度,由于 \(a(\mathbf{u},\mathbf{v})\) 的双线性形式,故只需考虑
\[
\int_{E} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_{+,j}) \, \mathrm{d}V\quad \text{或} \quad \int_{E} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{+,j}) \, \mathrm{d}V,
\]
由于基函数具有局部支撑性,因此在单元上的积分仅与该单元有限元节点上的基函数有关,故在计算单元积分时,只需考虑在该单元上取值不为零的基函数,且在单元积分运算时,只需考虑其在单元上的限制 \(\mathbf{v}_{*,i} = \mathbf{v}_{*,i}|_{E}\),以下在单元积分时,\(\mathbf{v}_{*,i}\) 默认为 \(\mathbf{v}_{*,i}|_{E}\)
通过选取与单元相关的 \(\mathbf{v}_{*,i}\) 和 \(\mathbf{v}_{*,j}\) 进行计算,得到单元刚度矩阵。由于物理单元的复杂性和多样性,物理单元上的积分通常会转换到参考单元上进行,此时需要利用形函数进行坐标变换和积分计算
形函数
记参考单元上的形函数为
\[
N_{1},N_{2},\cdots,N_{n}.
\]
几何映射
在物理单元 \(E_{\text{}}\) 上(以下为局部编号)
\[
\begin{equation}
x = \sum_{i=1}^{n} N_{i}(\xi,\eta,\zeta)x_{i},\quad y = \sum_{i=1}^{n} N_{i}(\xi,\eta,\zeta)y_{i},\quad z = \sum_{i=1}^{n} N_{i}(\xi,\eta,\zeta)z_{i}.
\end{equation}
\]
等参变换
根据等参变换,基函数
(6)\[
v_{i}(x,y,z) = N_{i}(x(\xi,\eta,\zeta),y(\xi,\eta,\zeta),z(\xi,\eta,\zeta)),
\]
有
(7)\[
\nabla_{\mathbf{x}}v_{i} = \nabla_{\mathbf{x}}N_{i} = \mathbf{J}^{-1}\nabla_{\boldsymbol{\xi}}N_{i}.
\]
Jacobian 矩阵
在物理单元 \(E\) 上(以下为局部编号)
\[\begin{split}
\begin{equation}
\begin{aligned}
\mathbf{J} =
\begin{bmatrix}
\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi} \\
\frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta} \\
\frac{\partial x}{\partial \zeta} & \frac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta}
\end{bmatrix}
&=
\sum_{i=1}^{n}
\begin{bmatrix}
\frac{\partial N_i}{\partial \xi} x_i & \frac{\partial N_i}{\partial \xi} y_i & \frac{\partial N_i}{\partial \xi} z_i \\
\frac{\partial N_i}{\partial \eta} x_i & \frac{\partial N_i}{\partial \eta} y_i & \frac{\partial N_i}{\partial \eta} z_i \\
\frac{\partial N_i}{\partial \zeta} x_i & \frac{\partial N_i}{\partial \zeta} y_i & \frac{\partial N_i}{\partial \zeta} z_i
\end{bmatrix}\\
&=\begin{bmatrix}
\frac{\partial N_{1}}{\partial \xi} & \frac{\partial N_{2}}{\partial \xi} & \cdots & \frac{\partial {N_n}}{\partial \xi} \\
\frac{\partial N_{1}}{\partial \eta} & \frac{\partial N_{2}}{\partial \eta} & \cdots & \frac{\partial N_{n}}{\partial \eta} \\
\frac{\partial N_{1}}{\partial \zeta} & \frac{\partial N_{2}}{\partial \zeta} & \cdots & \frac{\partial N_{n}}{\partial \zeta}
\end{bmatrix}
\begin{bmatrix}
x_{1} & y_{1} & z_{1} \\
x_{2} & y_{2} & z_{2} \\
\vdots & \vdots & \vdots \\
x_{n} & y_{n} & z_{n}
\end{bmatrix}.
\end{aligned}
\end{equation}
\end{split}\]
场变量插值
在物理单元 \(E\) 上(以下为局部编号)
\[\begin{split}
\begin{equation}
\mathbf{u}
=
\begin{bmatrix}
u_{x} \\ u_{y} \\ u_{z}
\end{bmatrix}
=
\begin{bmatrix}
\sum_{i=1}^{n}N_{i}(\xi,\eta,\zeta)u_{x,i} \\
\sum_{i=1}^{n}N_{i}(\xi,\eta,\zeta)u_{y,i} \\
\sum_{i=1}^{n}N_{i}(\xi,\eta,\zeta)u_{z,i}
\end{bmatrix}
=
\sum_{i=1}^{n}N_{i}\mathbf{u}_{i}.
\end{equation}
\end{split}\]
梯度计算
在物理单元 \(E\) 上(以下为局部编号),对于场分布函数
\[\begin{split}
\begin{equation}
\begin{aligned}
\begin{bmatrix}
\frac{\partial u_{x}}{\partial x} & \frac{\partial u_{y}}{\partial x} & \frac{\partial u_{z}}{\partial x} \\
\frac{\partial u_{x}}{\partial y} & \frac{\partial u_{y}}{\partial y} & \frac{\partial u_{z}}{\partial y} \\
\frac{\partial u_{x}}{\partial z} & \frac{\partial u_{y}}{\partial z} & \frac{\partial u_{z}}{\partial z}
\end{bmatrix} &= \begin{bmatrix}
\frac{\partial x}{\partial \xi} & \frac{\partial y}{\partial \xi} & \frac{\partial z}{\partial \xi} \\
\frac{\partial x}{\partial \eta} & \frac{\partial y}{\partial \eta} & \frac{\partial z}{\partial \eta} \\
\frac{\partial x}{\partial \zeta} & \frac{\partial y}{\partial \zeta} & \frac{\partial z}{\partial \zeta}
\end{bmatrix}^{-1}
\begin{bmatrix}
\frac{\partial u_{x}}{\partial \xi} & \frac{\partial u_{y}}{\partial \xi} & \frac{\partial u_{z}}{\partial \xi} \\
\frac{\partial u_{x}}{\partial \eta} & \frac{\partial u_{y}}{\partial \eta} & \frac{\partial u_{z}}{\partial \eta} \\
\frac{\partial u_{x}}{\partial \zeta} & \frac{\partial u_{y}}{\partial \zeta} & \frac{\partial u_{z}}{\partial \zeta}
\end{bmatrix}\\
&=\mathbf{J}^{-1}
\begin{bmatrix}
\frac{\partial N_{1}}{\partial \xi} & \frac{\partial {N_2}}{\partial \xi} & \cdots & \frac{\partial {N_n}}{\partial \xi} \\
\frac{\partial N_{1}}{\partial \eta} & \frac{\partial N_{2}}{\partial \eta} & \cdots & \frac{\partial N_{n}}{\partial \eta} \\
\frac{\partial N_{1}}{\partial \zeta} & \frac{\partial N_{2}}{\partial \zeta} & \cdots & \frac{\partial N_{n}}{\partial \zeta}
\end{bmatrix}
\begin{bmatrix}
u_{x,1} & u_{y,1} & u_{z,1} \\
u_{x,2} & u_{y,2} & u_{z,2} \\
\vdots & \vdots & \vdots \\
u_{x,n} & u_{y,n} & u_{z,n}
\end{bmatrix}.
\end{aligned}
\end{equation}
\end{split}\]
对于形函数,例如
(8)\[\begin{split}
\begin{equation}
\begin{bmatrix}
\frac{\partial N_{1}(x,y,z)}{\partial x} \\
\frac{\partial N_{1}(x,y,z)}{\partial y} \\
\frac{\partial N_{1}(x,y,z)}{\partial z}
\end{bmatrix}
=
\mathbf{J}^{-1}
\begin{bmatrix}
\frac{\partial N_{1}(\xi,\eta,\zeta)}{\partial \xi} \\
\frac{\partial N_{1}(\xi,\eta,\zeta)}{\partial \eta} \\
\frac{\partial N_{1}(\xi,\eta,\zeta)}{\partial \zeta}
\end{bmatrix}.
\end{equation}
\end{split}\]
单元刚度矩阵
张量形式
将物理单元上的积分运算映射回参考单元
\[
\int_{E} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_{+,j}) \, \mathrm{d}V = \int_{E_{\text{参考}}} \boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_{+,j})\left|\det(\mathbf{J})\right| \, \mathrm{d}V_{\text{参考}}
\]
需要注意的是,左端式子的坐标是 \((x,y,z)\),而右端式子的坐标是 \((\xi,\eta,\zeta)\),转换到参考单元上的计算见 (6) 和 (7)
考虑
\[
\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_{+,j}),
\]
由于
\[
\mathbb{C} : \boldsymbol{\varepsilon} = \mathbb{C}:\boldsymbol{\varepsilon} = \lambda\text{tr}(\boldsymbol{\varepsilon})\mathbf{I} + 2\mu\boldsymbol{\varepsilon},
\]
故
\[
\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) : \mathbb{C} : \boldsymbol{\varepsilon}(\mathbf{v}_{+,j}) = \lambda\text{tr}(\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}))\text{tr}(\boldsymbol{\varepsilon}(\mathbf{v}_{+,j})) + 2\mu\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}):\boldsymbol{\varepsilon}(\mathbf{v}_{+,j}),
\]
其中
\[
\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}) = \frac{1}{2}\left(\nabla\mathbf{v}_{*,i} + \nabla\mathbf{v}_{*,i}^{T}\right).
\]
使用数值积分,得到
(9)\[
\mathbf{K}^{e}_{(*,i),(+,j)} \approx \sum_{q}\left(\lambda\text{tr}(\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}))\text{tr}(\boldsymbol{\varepsilon}(\mathbf{v}_{+,j})) + 2\mu\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}):\boldsymbol{\varepsilon}(\mathbf{v}_{+,j})\right)\cdot w_{q} \left|\det(\mathbf{J}_{q})\right|.
\]
通常 \(\mathbf{v}_{*,i}\) 的排序为
\[
\mathbf{v}_{x,1},\mathbf{v}_{y,1},\mathbf{v}_{z,1},\mathbf{v}_{x,2},\mathbf{v}_{y,2},\mathbf{v}_{z,2},\cdots
\]
由于
\[
\text{tr}(\boldsymbol{\varepsilon}(\mathbf{u})) = \frac{1}{2}\text{tr}(\nabla\mathbf{u} + \nabla\mathbf{u}^{T}) = \nabla\cdot\mathbf{u},
\]
故式 (9) 也可以写为
(10)\[
\mathbf{K}^{e}_{(*,i),(+,j)} \approx \sum_{q}\left(\lambda(\nabla\cdot\mathbf{v}_{*,i})(\nabla\cdot\mathbf{v}_{+,j}) + 2\mu\boldsymbol{\varepsilon}(\mathbf{v}_{*,i}):\boldsymbol{\varepsilon}(\mathbf{v}_{+,j})\right)\cdot w_{q} \left|\det(\mathbf{J}_{q})\right|.
\]
Vogit 形式
将物理单元上的积分运算映射回参考单元
(11)\[
\int_{E} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{+,j})\, \mathrm{d}V = \int_{E_{\text{参考}}} (\mathcal{B}\mathbf{v}_{*,i})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{+,j})\left|\det(\mathbf{J})\right|\, \mathrm{d}V_{\text{参考}},
\]
需要注意的是,左端式子的坐标是 \((x,y,z)\),而右端式子的坐标是 \((\xi,\eta,\zeta)\)
\(\mathcal{B}\mathbf{v}_{*,i}\) 转换到参考单元后的梯度计算公式见 (7)
接下来,我们将 (11) 转换为矩阵形式,即计算 \((\mathbf{v}_{*,i},\mathbf{v}_{+,j})\) 的所有组合,将它们依次选取为
\[
\mathbf{v}_{x,1},\mathbf{v}_{y,1},\mathbf{v}_{z,1},\mathbf{v}_{x,2},\mathbf{v}_{y,2},\mathbf{v}_{z,2},\cdots
\]
于是,式 (11) 写为
\[\begin{split}
\begin{equation}
\begin{bmatrix}
(\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{x,1}) & (\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{y,1}) & (\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{z,1}) & \cdots \\
(\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{x,1}) & (\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{y,1}) & (\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{z,1}) & \cdots \\
(\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{x,1}) & (\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{y,1}) & (\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}(\mathcal{B}\mathbf{v}_{z,1}) & \cdots \\
\vdots & \vdots & \vdots & \ddots
\end{bmatrix},
\end{equation}
\end{split}\]
进一步写为
\[\begin{split}
\begin{equation}
\begin{bmatrix}
(\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,1}\ \mathbf{v}_{y,1}\ \mathbf{v}_{y,1}\right]) & (\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,2}\ \mathbf{v}_{y,2}\ \mathbf{v}_{y,2}\right]) & \cdots \\
(\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,1}\ \mathbf{v}_{y,1}\ \mathbf{v}_{y,1}\right]) & (\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,2}\ \mathbf{v}_{y,2}\ \mathbf{v}_{y,2}\right]) & \cdots \\
(\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,1}\ \mathbf{v}_{y,1}\ \mathbf{v}_{y,1}\right]) & (\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}(\mathcal{B}\left[\mathbf{v}_{x,2}\ \mathbf{v}_{y,2}\ \mathbf{v}_{y,2}\right]) & \cdots \\
\vdots & \vdots & \ddots
\end{bmatrix},
\end{equation}
\end{split}\]
记
\[\begin{split}
\mathbf{B}_{i} = \mathcal{B}\left[\mathbf{v}_{x,i}\ \mathbf{v}_{y,i}\ \mathbf{v}_{y,i}\right] = \begin{bmatrix}
\frac{\partial N_{i}}{\partial x} & 0 & 0 \\
0 & \frac{\partial N_{i}}{\partial y} & 0 \\
0 & 0 & \frac{\partial N_{i}}{\partial z} \\
\frac{\partial N_{i}}{\partial y} & \frac{\partial N_{i}}{\partial x} & 0 \\
0 & \frac{\partial N_{i}}{\partial z} & \frac{\partial N_{i}}{\partial y} \\
\frac{\partial N_{i}}{\partial z} & 0 & \frac{\partial N_{i}}{\partial x}
\end{bmatrix},
\end{split}\]
且
\[
\mathbf{B}=
\begin{bmatrix}
\mathbf{B}_{1} & \mathbf{B}_{2} & \cdots & \mathbf{B}_{n}
\end{bmatrix},
\]
Note
记
\[\begin{split}
\mathbf{N} = \begin{bmatrix}
N_{1} & & & & N_{n} \\
& N_{1} & & \cdots & & N_{n} \\
& & N_{1} & & & & N_{n}
\end{bmatrix},
\end{split}\]
于是
\[
\mathbf{B} = \mathcal{B}\mathbf{N}.
\]
于是,上述矩阵进一步写为
\[\begin{split}
\begin{equation}
\begin{bmatrix}
(\mathcal{B}\mathbf{v}_{x,1})^{T}\mathbb{C}\mathbf{B}\\
(\mathcal{B}\mathbf{v}_{y,1})^{T}\mathbb{C}\mathbf{B}\\
(\mathcal{B}\mathbf{v}_{z,1})^{T}\mathbb{C}\mathbf{B}\\
\vdots
\end{bmatrix} = \begin{bmatrix}
(\mathcal{B}\mathbf{v}_{x,1})^{T}\\
(\mathcal{B}\mathbf{v}_{y,1})^{T}\\
(\mathcal{B}\mathbf{v}_{z,1})^{T}\\
\vdots
\end{bmatrix}\mathbb{C}\mathbf{B} = \mathbf{B}^{T}\mathbb{C}\mathbf{B}
\end{equation}
\end{split}\]
带入积分,得到单元刚度矩阵,记为
\[\begin{split}
\begin{equation}
\begin{aligned}
\mathbf{K}^{e} &= \int_{E_{\text{参考}}} \mathbf{B}^{T}\mathbb{C}\mathbf{B}\left|\det(\mathbf{J})\right|\, \mathrm{d}V\\
&\approx\sum_{q}\mathbf{B}_{q}^{T}\mathbb{C}\mathbf{B}_{q}\cdot w_{q}\cdot\left|\det(\mathbf{J}_{q})\right|
\end{aligned}
\end{equation}
\end{split}\]
其中,积分点是 \(\left\{(\xi_{q},\eta_{q},\zeta_{q})\right\}\),积分权重是 \(w_{q}\)
边界条件
在有限元计算中,边界条件通常在组装完成全局刚度矩阵和右端项后进行附加处理
Dirichlet 边界条件
Dirichlet 边界条件直接指定了边界节点的自由度值,因此只需对全局刚度矩阵和右端项进行相应的修正
\[\begin{split}
\begin{array}{c|ccc}
& \cdots & \cdots & u_{x,k} & \cdots & \cdots\\
\hline
\vdots & & & 0 & & & \\
\vdots & & & 0 & & & \\
u_{x,k} & 0 & 0 & 1 & 0 & 0 & \\
\vdots & & & 0 & & & \\
\vdots & & & 0 & & & \\
\end{array},\quad\quad\quad\quad
\begin{bmatrix}
\\
\vdots\\
\vdots\\
u_{D}\\
\vdots\\
\vdots\\
\end{bmatrix},
\end{split}\]
其中 \(u_{D}\) 是边界节点自由度 \(u_{x,k}\) 的值
Neumann 边界条件
边界积分通常通过划分至对应的边界单元表面进行计算
\[
\int_{\Gamma_{N}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i} \ \mathrm{d}S = \sum_{\Gamma_{N}^{E}}\int_{\Gamma_{N}^{E}} \tilde{\mathbf{p}} \cdot \mathbf{v}_{*,i} \ \mathrm{d}\Gamma_{N}^{E},\quad *=x,y,z;\ i=1:N.
\]
其中,\(\Gamma_{N}^{E}\) 表示单元 \(E\) 与边界 \(\Gamma_{N}\) 的交集面。类似地,可以将物理单元上的面积分转换为参考单元上的面积分,但需要事先确定物理单元上的边界面与参考单元上哪个面相对应(例如六面体单元的 \(\xi=1\) 面或 \(\eta=-1\) 面),此时,积分通过坐标变换转换到参考面上计算
\[\begin{split}
\begin{equation}
\begin{aligned}
\mathbf{F}_{\Gamma_{N}^{E}} &= \int_{\Gamma_{N}^{E}} \mathbf{N}^{T}\tilde{\mathbf{p}}\ \mathrm{d}\Gamma_{N}^{E}
=\int_{\partial E_{\text{参考}}} \mathbf{N}^{T}\tilde{\mathbf{p}}\ \left|\det(\mathbf{J}^{(\xi,\eta)})\right| \mathrm{d} S \\
&= \sum_{q}\mathbf{N}^{T}_{q}\tilde{\mathbf{p}}_{q}\cdot w_{q}\cdot\left|\det(\mathbf{J}^{(\xi,\eta)}_{q})\right|,
\end{aligned}
\end{equation}
\end{split}\]
其中,积分点是 \(\left\{(\xi_{q},\eta_{q})\right\}\),积分权重是 \(w_{q}\),\(\tilde{\mathbf{p}}\) 需使用 \(\xi,\eta\) 坐标表示
单元右端项到全局右端项的映射关系如下
\[
\mathbf{F}(g(E,i)) \ +\!\!= \ \mathbf{F}_{\Gamma_{N}^{E}}(i).
\]
Robin 边界条件
Robin 边界条件通常有两项需要处理
\[
\int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S
\quad \text{和} \quad
\int_{\Gamma_{R}} \mathbf{f}_{R} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S,\quad \forall \, \mathbf{v}_{*,i} \in \, \mathcal{V},\quad *=x,y,z;\ i=1:N.
\]
第二项的处理与 Neumann 边界条件类似。对于第一项
\[
\int_{\Gamma_{R}} \alpha\mathbf{u} \cdot \mathbf{v}_{*,i} \, \mathrm{d}S = \sum_{\Gamma_{R}^{E}}\int_{\Gamma_{R}^{E}} \alpha \mathbf{u} \cdot \mathbf{v}_{*,i} \ \mathrm{d}\Gamma_{R}^{E},\quad *=x,y,z;\ i=1:N.
\]
因为
\[\begin{split}
\mathbf{u} = \sum_{i=1}^{n}N_{i}\mathbf{u}_{i} =
\begin{bmatrix}
N_{1} & & & & N_{n} \\
& N_{1} & & \cdots & & N_{n} \\
& & N_{1} & & & & N_{n}
\end{bmatrix}\begin{bmatrix}
u_{x,1} \\ u_{y,1} \\ u_{z,1} \\ u_{x,2} \\ u_{y,2} \\ u_{z,2} \\ \vdots \\ u_{x,n} \\ u_{y,n} \\ u_{z,n}
\end{bmatrix}=\mathbf{N}\mathbf{u}_{E},
\end{split}\]
对于单元 \(E\) ,积分写为矩阵形式
\[\begin{split}
\begin{equation}
\begin{aligned}
\mathbf{F}_{\Gamma_{R}^{E}}
&=
\int_{\Gamma_{R}^{E}} \alpha\mathbf{N}^{T}\mathbf{N}\mathbf{u}_{E} \ \mathrm{d}\Gamma_{R}^{E}\\
&=
\int_{\Gamma_{R}^{E}} \alpha\mathbf{N}^{T}\mathbf{N} \ \mathrm{d}\Gamma_{R}^{E} \cdot \mathbf{u}_{E}\\
&=\int_{\partial E_{\text{参考}}} \alpha\mathbf{N}^{T}\mathbf{N}\ \left|\det(\mathbf{J}^{(\xi,\eta)})\right| \mathrm{d} S \ \mathbf{u}_{E} \\
&= \sum_{q}\alpha\mathbf{N}^{T}_{q}\mathbf{N}_{q}\cdot w_{q}\cdot\left|\det(\mathbf{J}^{(\xi,\eta)}_{q})\right|\mathbf{u}_{E}\\
&=\mathbf{K}_{r}^{e}\cdot \mathbf{u}_{E},
\end{aligned}
\end{equation}
\end{split}\]
其中,积分点是 \(\left\{(\xi_{q},\eta_{q})\right\}\),积分权重是 \(w_{q}\)
单元右端项到全局右端项的映射关系如下
\[
\mathbf{F}(g(E,i)) \ +\!\!= \ \mathbf{F}_{\Gamma_{R}^{E}}(i).
\]