[Peridynamics] 各向同性均匀材料的耦合有限元方法

• 方程与符号说明
  • 几何方程
  • 力平衡方程
  • 本构方程
• 虚功原理
• 计算流程

各向异性的推导与各向同性的类似，这里只讨论各向同性均匀材料的耦合有限元方法的推导。

方程与符号说明

对于耦合方法，我们现在集中说明所有的方程和所用到的所有符号。

几何方程

在连续介质力学，我们有几何方程

$$\underline{\underline{\epsilon}}(\underline{x}) = \frac12 ( \underline{\underline{\nabla}} \cdot \underline{u}(\underline{x}) + {}^t\underline{\underline{\nabla}} \cdot \underline{u}(\underline{x}) ) \quad \forall \underline{x} \in \Omega$$

连续介质力学的几何方程描述的是位移 $\underline{u}(\underline{x})$ 和应变 $\underline{\underline{\epsilon}}(\underline{x})$ 之间的关系。

而在近场动力学模型下的几何方程，实际上描述的是一对物质点在键力（或力态）的作用下所发生的的位移。即

$$\eta_{\underline{\xi}}(\underline{p} - \underline{x}) = u_{\underline{\xi}}(\underline{p}) - u_{\underline{\xi}}(\underline{x}) \quad \forall (\underline{x}, \underline{p}) \in \Omega$$

在这个公式里，$\underline{\xi}$ 表示一对物质点的键向量 $\underline{\xi} = \underline{p} - \underline{x}$。而 $u_{\underline{\xi}}(\underline{x})$ 表示点 $\underline{x}$ 发生的位移 $\underline{u}(\underline{x})$ 在这个键上的投影，即 $u_{\underline{\xi}}(\underline{p}) = \underline{u}(\underline{p}) \cdot \underline{e}_{\underline{\xi}}$，其中 ${e}_{\underline{\xi}}$ 表示单位方向向量，即 ${e}_{\underline{\xi}} = \underline{\xi} / \|\underline{\xi}\|$。我们用 $\eta_{\underline{\xi}}$ 描述一对物质点的位移差在键上的投影。同样也可以写做 $\underline{\eta}(\underline{p} - \underline{x}) = \underline{u}(\underline{p}) - \underline{u}(\underline{x})$ 表示位移差。

最后在边界上需要保持一致性，所以额外有位移边界条件

$$\underline{u} = \overline{\underline u} \quad \forall \underline{x} \in S_{\overline{\underline u}}$$

力平衡方程

对于整个系统，需要保持内力与外力的静力平衡，否则由牛顿第二定律可知系统表现为有加速度项。在这里我们只考虑静力平衡方程。

$$\underline{div} \cdot \underline{\underline{\sigma}} + \int_{H_\delta(\underline{x})} \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} V_{\underline{p}} + \underline{b} = \underline{0} \quad \forall \underline{x} \in \Omega$$

在这里，可以这样理解：第一项表示连续介质力学所描述的内力，第二项表示近场动力学模型所描述的物质点之间的内力，最后一项表示体力。三项均是一阶张量，单位为 $\mathrm N/\mathrm m^3$。

公式中，$\underline{\underline{\sigma}}$ 表示应变张量，可以通过连续介质力学的本构计算得到。$\underline{\hat{f}}(\underline{p}, \underline{x})$ 表示近场动力学模型所描述的在物质点 $\underline{x}$ 的近场邻域 $H_\delta(\underline{x})$ 内的另一物质点 $\underline{p}$ 对 $\underline{x}$ 的作用。两者的差刻画了这一对物质点的相互作用。

在边界上，我们只考虑连续介质力学模型，所以有边界条件

$$\underline{\underline{\sigma}} \cdot \underline{n} = \overline{\underline T} \quad \forall \underline{x} \in \Omega_{\overline{\underline T}}$$

这里的 $\underline{n}$ 表示外法向方向。

本构方程

由于耦合方法使用的是两套模型，所以我们同时需要连续介质力学的本构方程及近场动力学模型的本构方程。

$$\underline{\underline{\sigma}} = \underline{\underline{\underline{\underline{E}}}}(\underline{x}) : \underline{\underline{\epsilon}} \quad \forall \underline{x} \in \Omega$$

连续介质力学的本构方程描述的是应变与应力之间的关系，即刚度张量，由一个四阶张量 $\underline{\underline{\underline{\underline{E}}}}(\underline{x})$ 来表示。

近场动力学模型的本构描述的是物质点的相互作用与位移差的关系，在局部等应变假设的前提下，我们将相互作用做一阶泰勒展开（线性化处理），从而得到近场动力学模型的本构方程

$$\begin{aligned} \underline{\hat{f}}(\underline{p}, \underline{x}) &= \frac12 \cdot \underline{\underline{C}}(\underline{x}, \underline{\xi}) \cdot \underline{\eta}(\underline{p}, \underline{x}) \\ &= \frac12 \cdot c(\underline{x}, \underline{\xi}) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) \\ \end{aligned}$$

在此基础上，我们引入耦合函数 $\alpha(x)$，即 $\underline{\underline{C}}(\underline{x}, \underline{\xi}) = \alpha(x) \cdot c^0(\|\underline{\xi}\|) \cdot \underline{\xi} \otimes \underline{\xi}$，整理近场动力学模型的本构

$$\begin{aligned} \underline{f}(\underline{p} \rightarrow \underline{x}) &= \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \\ &= \frac12 \cdot c(\underline{x}, \underline{\xi}) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) -\frac12 \cdot c(\underline{p}, \underline{\xi}) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{x}, \underline{p}) \\ &= \frac12 \cdot (c(\underline{x}, \underline{\xi}) + c(\underline{p}, \underline{\xi})) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) \\ &= \frac12 \cdot c^0(\|\underline{\xi}\|) \cdot (\alpha(\underline{x}) + \alpha(\underline{p})) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) \end{aligned}$$

连续介质力学的本构和近场动力学模型的本构所联立的方程组是不可解的，因为连续介质力学中的刚度张量 $\underline{\underline{\underline{\underline{E}}}}(\underline{x})$ 并不完全等效于近场动力学本构方程中的 $c^0(\|\underline{\xi}\|)$，即我们仍然需要一个方程来建立两者的联系。在 Lubineau et al. (2012) 这篇文章中，提到了可以通过应变能密度的一致性来推导两个系数之间的联系，即

$$\begin{aligned} \frac12 \underline{\underline{\epsilon}}(\underline{x}) : \underline{\underline{\underline{\underline{E}}}}^0 : \underline{\underline{\epsilon}}(\underline{x}) &= \frac12 \underline{\underline{\epsilon}}(\underline{x}) : \underline{\underline{\underline{\underline{E}}}}(\underline{x}) : \underline{\underline{\epsilon}}(\underline{x}) \\ &\quad+ \frac14 \int_{H_\delta(\underline{x})}c^0(\|\underline{\xi}\|) \frac{ \alpha(\underline{x}) + \alpha(\underline{p}) }{2} \{ \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) \}^2 {\mathrm{d}\!} V_{\underline{p}} \end{aligned}$$

化简得到

$$\begin{aligned} \underline{\underline{\underline{\underline{E}}}}(\underline{x}) = (1 - \alpha(\underline{x}))\underline{\underline{\underline{\underline{E}}}}^0 + \int_{H_\delta(\underline{x})}c^0(\|\underline{\xi}\|) \frac{ \alpha(\underline{x}) - \alpha(\underline{p}) }{2} \cdot \frac{\underline{\xi} \otimes \underline{\xi} \otimes \underline{\xi} \otimes \underline{\xi}}{2\|\underline{\xi}\|^2} {\mathrm{d}\!} V_{\underline{p}} \end{aligned}$$

从而得到全局意义下的等效刚度张量 $\underline{\underline{\underline{\underline{E}}}}(\underline{x})$，即只要提供一个合理的 $\alpha(x)$ 就可以直接计算出整个求解域内任意点的等效刚度张量。

虚功原理

我们把整个区域进行网格剖分成 $n$ 个单元，同时依据 Lubineau et al. (2012) 文中所说，将整个求解域分割成三个大块，纯连续介质力学区域 $\Omega_1$，纯近场动力学区域 $\Omega_2$ 和耦合区域 $\Omega_m$。

回到我们的平衡方程

$$\underline{div} \cdot \underline{\underline{\sigma}} + \int_{H_\delta(\underline{x})} \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} V_{\underline{p}} + \underline{b} = \underline{0} \quad \forall \underline{x} \in \Omega$$

和

$$\underline{\underline{\sigma}} \cdot \underline{n} = \overline{\underline T} \quad \forall \underline{x} \in \Omega_{\overline{\underline T}}$$

如果从虚功原理进行推导的话，相当于系统的所有力和其对应所产生的虚位移的总做功为 0。我们令 $\underline{\delta u}(\underline{x})$ 表示虚位移场，同时要求在边界上，虚位移场 $\underline{\delta u}(\underline{x}) = \underline{0}$，$\forall \underline{x} \in \Gamma_u$。则此时我们的平衡方程在区域 $\Omega$ 及其边界 $\Omega_{\overline{\underline T}}$ 上所做的虚功为：

$$\int_\Omega \underline{\delta u}(\underline{x}) \cdot ( \underline{div} \cdot \underline{\underline{\sigma}} + \underline{b}(\underline{x}) + \underline{f}(\underline{p} \rightarrow \underline{x}) ) {\mathrm{d}\!} \Omega - \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot ( \underline{\underline{\sigma}} \cdot \underline{n} - \overline{\underline T} ) {\mathrm{d}\!} S = 0$$

其中 $\underline{f}(\underline{p} \rightarrow \underline{x})$ 表示 $\underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p})$。

我们对上式的第一项逐项分析。$\underline{div} \cdot \underline{\underline{\sigma}}$ 表示应力在各个方向的分量，那么由格林公式我们可以进一步得到应力的虚功和应变能的关系，即区域内部的做功总和应当等于区域边界上做功总和，而这里区域内部的做功总和具体表现为虚位移产生的做功，和表现为应力应变的做功（即应变能）。所以我们有

$$\int_\Omega \underline{\delta u}(\underline{x}) \cdot ( \underline{div} \cdot \underline{\underline{\sigma}} ) {\mathrm{d}\!} \Omega = -\int_\Omega \underline{\underline{\delta \epsilon}}(\underline{x}) : \underline{\underline{\sigma}}(\underline{x}) {\mathrm{d}\!} \Omega + \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot \underline{\underline{\sigma}}(\underline{x}) \cdot \underline{n} {\mathrm{d}\!} S$$

把上式带入化简，得到

$$\begin{aligned} &-\int_\Omega \underline{\underline{\delta \epsilon}}(\underline{x}) : \underline{\underline{\sigma}}(\underline{x}) {\mathrm{d}\!} \Omega + \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot \underline{\underline{\sigma}}(\underline{x}) \cdot \underline{n} {\mathrm{d}\!} S + \int_\Omega \underline{\delta u}(\underline{x}) \cdot \underline{b}(\underline{x}) {\mathrm{d}\!} \Omega + \int_\Omega \underline{\delta u}(\underline{x}) \cdot \underline{f}(\underline{p} \rightarrow \underline{x}) {\mathrm{d}\!} \Omega \\ &- \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot\underline{\underline{\sigma}} \cdot \underline{n} {\mathrm{d}\!} S + \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot \overline{\underline T} {\mathrm{d}\!} S = 0 \end{aligned}$$

在这里，第 2 项与第 5 项相互抵消了，所以边界上实际只剩下了外部力与对应的虚位移的做功。

下面分析第四项，即近场动力学模型中，两个物质点的作用力所做虚功。

$$\begin{aligned} \int_\Omega \underline{\delta u}(\underline{x}) \cdot \underline{f}(\underline{p} \rightarrow \underline{x}) {\mathrm{d}\!} \Omega &= \int_\Omega \int_\Omega \underline{\delta u}(\underline{x}) \cdot \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ &= \int_\Omega \int_\Omega \underline{\delta u}(\underline{p}) \cdot \left \{ \underline{\hat{f}}(\underline{x}, \underline{p}) - \underline{\hat{f}}(\underline{p}, \underline{x}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ \end{aligned}$$

把两个等式综合到一起，即得

$$\begin{aligned} \int_\Omega \underline{\delta u}(\underline{x}) \cdot \underline{f}(\underline{p} \rightarrow \underline{x}) {\mathrm{d}\!} \Omega &= -\frac12 \int_\Omega \int_\Omega \left [ \underline{\delta u}(\underline{p}) - \underline{\delta u}(\underline{x}) \right ] \cdot \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ \end{aligned}$$

所以原式可以整理为如下四项

$$\begin{aligned} &+ \frac12 \int_\Omega \int_\Omega \left [ \underline{\delta u}(\underline{p}) - \underline{\delta u}(\underline{x}) \right ] \cdot \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ &+\int_\Omega \underline{\underline{\delta \epsilon}}(\underline{x}) : \underline{\underline{\sigma}}(\underline{x}) {\mathrm{d}\!} \Omega - \int_\Omega \underline{\delta u}(\underline{x}) \cdot \underline{b}(\underline{x}) {\mathrm{d}\!} \Omega - \int_{\Omega_{\overline{\underline T}}} \underline{\delta u}(\underline{x}) \cdot \overline{\underline T} {\mathrm{d}\!} S = 0 \end{aligned}$$

虚功原理的推导结果与最小势能原理的推导结果是一致的。如果以最小势能原理为出发点进行分析，相当于考察泛函的最小值，而这个最小值在他的一阶导为 0 时取到，即最终结论与虚功原理等价。此时的泛函为

$$\begin{aligned} \Pi(\underline{u}) &= \frac12 \int_\Omega \underline{\underline{\epsilon}}(\underline{x}) : \underline{\underline{\sigma}}(\underline{x}) {\mathrm{d}\!} \Omega - \int_\Omega \underline{u}(\underline{x}) \cdot \underline{b}(\underline{x}) {\mathrm{d}\!} \Omega - \int_{\Omega_{\overline{\underline T}}} \underline{u}(\underline{x}) \cdot \overline{\underline T} {\mathrm{d}\!} S \\ &\quad+ \frac14 \int_\Omega \int_\Omega \left [ \underline{u}(\underline{p}) - \underline{u}(\underline{x}) \right ] \cdot \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ &= W_c(\underline{\underline{\epsilon}}) + W_s(\underline{\eta}) - F_{\Omega}(\underline{u}) - F_{\Gamma}(\underline{u}) \\ W_c(\underline{\underline{\epsilon}}) &= \frac12 \int_\Omega \underline{\underline{\epsilon}}(\underline{x}) : \underline{\underline{\sigma}}(\underline{x}) {\mathrm{d}\!} \Omega_{\underline{x}} = \frac12 \int_\Omega \underline{\underline{\epsilon}}(\underline{x}) : \underline{\underline{\underline{\underline{E}}}}(\underline{x}) : \underline{\underline{\epsilon}}(\underline{x}) {\mathrm{d}\!} \Omega_{\underline{x}} \\ F_{\Omega}(\underline{u}) &= \int_\Omega \underline{b}(\underline{x}) \cdot \underline{u}(\underline{x}){\mathrm{d}\!} \Omega \\ F_{\Gamma}(\underline{u}) &= \int_{\Omega_{\overline{\underline T}}} \overline{\underline T} \cdot \underline{u}(\underline{x}) {\mathrm{d}\!} S \\ \end{aligned}$$

注意，这里的 $\underline{\underline{\underline{\underline{E}}}}(\underline{x})$ 项是全局意义下的等效刚度张量，即

$$\begin{aligned} \underline{\underline{\underline{\underline{E}}}}(\underline{x}) = (1 - \alpha(\underline{x}))\underline{\underline{\underline{\underline{E}}}}^0 + \int_{H_\delta(\underline{x})}c^0(\|\underline{\xi}\|) \frac{ \alpha(\underline{x}) - \alpha(\underline{p}) }{2} \cdot \frac{\underline{\xi} \otimes \underline{\xi} \otimes \underline{\xi} \otimes \underline{\xi}}{2\|\underline{\xi}\|^2} {\mathrm{d}\!} V_{\underline{p}} \end{aligned}$$

我们单独考虑 $W_s(\underline{\eta})$ 这一项，化简得到

$$\begin{aligned} W_s(\underline{\eta}) &= \frac14 \int_\Omega \int_\Omega \underline{\eta}(\underline{p}, \underline{x}) \cdot \left \{ \underline{\hat{f}}(\underline{p}, \underline{x}) - \underline{\hat{f}}(\underline{x}, \underline{p}) \right \} {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ &= \frac14 \int_\Omega \int_\Omega \underline{\eta}(\underline{p}, \underline{x}) \cdot c^0(\|\underline{\xi}\|) \cdot \frac{\alpha(\underline{x}) + \alpha(\underline{p})}{2} \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ &= \frac14 \int_\Omega \int_\Omega c^0(\|\underline{\xi}\|) \cdot \frac{\alpha(\underline{x}) + \alpha(\underline{p})}{2} \cdot \underline{\eta}(\underline{p}, \underline{x}) \cdot \underline{\xi} \otimes \underline{\xi} \cdot \underline{\eta}(\underline{p}, \underline{x}) {\mathrm{d}\!} \Omega_{\underline{p}} {\mathrm{d}\!} \Omega_{\underline{x}} \\ \end{aligned}$$

整理完毕后，回到最初的能量泛函

$$\begin{aligned} \Pi(\underline{u}) &= W_c(\underline{\underline{\epsilon}}) + W_s(\underline{\eta}) - F_{\Omega}(\underline{u}) - F_{\Gamma}(\underline{u}) \end{aligned}$$

那么，将上式在 $\underline{u} + \underline{\delta u}$ 进行泰勒展开，即有

$$\begin{aligned} \Pi(\underline{u} + \underline{\delta u}) &= \Pi(\underline{u}) + \frac{\partial \Pi(\underline{u})}{\partial \underline{u}} \cdot \underline{\delta u} + \frac12 \underline{\delta u}^T \cdot \underline{\underline{H}}[\Pi]( u + \underline{\underline{\Psi}}\underline{\delta u} ) \cdot \underline{\delta u} \geq \Pi(\underline{u}) \\ \end{aligned}$$

其中 $\underline{\underline{\Psi}}$ 是一个对角矩阵 $\underline{\underline{\Psi}} = \text{diag}(\psi_1, \psi_2, \psi_3)$ 并且 $\psi_i \in [0, 1]$。 $\underline{\underline{H}}[\Pi](u + \underline{\underline{\Psi}}\underline{\delta u})$ 表示 $\Pi(\underline{u})$ 在 $u + \underline{\underline{\Psi}}\underline{\delta u}$ 处的黑塞矩阵。

如果我们能够证明 $\frac12 \underline{\delta u}^T \cdot \underline{\underline{H}}[\Pi](u + \underline{\underline{\Psi}}\underline{\delta u}) \cdot \underline{\delta u} \geq 0$，这样，最小势能原理所给出的解 $\underline{u}$ 的充分必要条件就可以表达为

$$\frac{\partial \Pi(\underline{u})}{\partial \underline{u}} \cdot \underline{\delta u} = 0, \quad \forall \underline{\delta u}$$

即可以通过计算刚度矩阵 $\underline{\underline{K}}$ 和载荷向量 $\underline{F}$ 来组成类似于如下形式的方程

$$\Pi(\underline{u}) = \frac12 \underline{u}^T \underline{\underline{K}} \underline{u} - \underline{u}^T \underline{F} = \text{Minimum}$$

从而可以直接得到 $\underline{\underline{K}}\underline{u} - \underline{F} = 0$ 这样形式的线性方程组，进一步可以得到 $\underline{u}$

计算流程

符号	说明
$\Omega$	求解区域
$V_i$, $i = 1, 2, \ldots, n$	$n$ 个单元
$S_i$, $i = 1, 2, \ldots, n$	单元边界
$\mathbf{d}_i$	单元 $i$ 顶点的位移
$\mathbf{N}_i(x)$	形函数矩阵
$\mathbf{u}_i(x)$	单元 $i$ 的位移
$\mathbf{H}$	微分算子
$\mathbf{R}_i$	单元与节点的关联矩阵

其中

$$\begin{aligned} \mathbf{d}_i &= \mathbf{R}_i\mathbf{d} \\ \mathbf{u}_i(x) &= \mathbf{N}_i(x)\mathbf{R}_i\mathbf{d} \\ \mathbf{R}_i &= \text{diag}(0, \ldots, 1, \ldots, 0) \\ \mathbf{R}_i[\text{nodeid}(V_i)] &= 1 \end{aligned}$$

借助以上的信息可以组建刚度矩阵 $\mathbf{K}$ 和载荷向量 $\mathbf{F}$

$$\begin{aligned} \mathbf{K} &= \sum_{i=1}^n \int_{V_i} [\mathbf{H} \mathbf{N}_i(x) \mathbf{R}_i]^T \cdot [\mathbf{E}(x)] \cdot [\mathbf{H} \mathbf{N}_i(x) \mathbf{R}_i] {\mathrm{d}\!} V_x \\ &\quad+ \frac12 \sum_{i=1}^n \sum_{j=1}^{h(x)} \int_{V_i} \int_{V_x^j} c^0(\|\underline{\xi}\|) \cdot \frac{\alpha(x) + \alpha(x')}{2} \cdot [\mathbf{N}_j(x')\mathbf{R}_j - \mathbf{N}_i(x)\mathbf{R}_i]^T \\ &\quad\quad \cdot [\underline{\xi} \otimes \underline{\xi}] \cdot [\mathbf{N}_j(x')\mathbf{R}_j - \mathbf{N}_i(x)\mathbf{R}_i] {\mathrm{d}\!} V_{x'} {\mathrm{d}\!} V_x \\ \mathbf{F} &= \sum_{i=1}^n \int_{V_i} [\mathbf{N}_i(x) \mathbf{R}_i]^T \cdot [\underline{b}(x)] {\mathrm{d}\!} V_x + \sum_{i=1}^n \int_{S_i} [\mathbf{N}_i(x) \mathbf{R}_i]^T \cdot [\overline{\underline T}(x)] {\mathrm{d}\!} S_x \end{aligned}$$

所以有 $\mathbf{K} \mathbf{d} = \mathbf{F}$。