二维三角形单元

我们现将有限元方法的原理扩展到二维域。离散二维连续体最简单的单元是三角形单元。考虑三角形的每个角点处有三个节点，节点按逆时针顺序编号。每个节点有两个自由度 $u_{i}$ 和 $v_{i}$ 。

1. 二维弹性力学公式

在推导该单元的特性之前，我们必须建立二维弹性力学的基本关系。

任何点的应力和应变状态都可以用向量表示。

Voigt应力向量定义为： ${σ} = {\begin{matrix} σ_{x x} \\ σ_{y y} \\ σ_{x y} \end{matrix}}$

相应的Voigt应变向量为： ${ϵ} = {\begin{matrix} ϵ_{x x} \\ ϵ_{y y} \\ γ_{x y} \end{matrix}}$

应力与应变之间的关系是本构定律，对于线弹性材料，可以用矩阵形式表示为： ${σ} = [E] {ϵ}$ 这里，E是弹性矩阵，包含材料属性。（该矩阵不要与其定义中使用的标量杨氏模量E相混淆）。

对于二维问题，通常采用两种简化假设之一：平面应力或平面应变。

平面应力：假设用于在其自身平面内受载的薄结构（ $σ_{z z} = 0$ ）。弹性矩阵为： $𝐄 = \frac{E}{1 - ν^{2}} [\begin{matrix} 1 & ν & 0 \\ ν & 1 & 0 \\ 0 & 0 & \frac{1 - ν}{2} \end{matrix}]$

平面应变：假设用于具有均匀横截面的长结构，其面外应变为零（ $ϵ_{z z} = 0$ ）。弹性矩阵为： $𝐄 = \frac{E (1 - ν)}{(1 + ν) (1 - 2 ν)} [\begin{matrix} 1 & \frac{ν}{1 - ν} & 0 \\ \frac{ν}{1 - ν} & 1 & 0 \\ 0 & 0 & \frac{1 - 2 ν}{2 (1 - ν)} \end{matrix}]$

2 应变-位移矩阵 [B]

假设：单元内的位移场，x方向的u(x, y)和y方向的v(x, y)，通过形函数N_i从节点位移值插值得到。假设给定方向的位移仅取决于同方向的节点位移。即，

应变是位移场的导数。对于二维问题，关系如下：

$ϵ_{x x} = \frac{\partial u}{\partial x}, ϵ_{y y} = \frac{\partial v}{\partial y}, γ_{x y} = \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x}$

通过将插值后的位移场代入这些应变定义，我们可以建立Voigt应变向量 ${ϵ}$ 与节点位移向量 $𝐪$ 之间的关系。此关系定义了应变-位移矩阵B。

${ϵ} = [B] {q}$

对于3节点三角形单元，节点位移向量排列为：

$𝐪 = {\begin{matrix} u_{1} \\ v_{1} \\ u_{2} \\ v_{2} \\ u_{3} \\ v_{3} \end{matrix}}$

因此，应变-位移矩阵B是一个 $3 \times 6$ 矩阵，由形函数的导数组成：

$𝐁 = [\begin{matrix} \frac{\partial N_{1}}{\partial x} & 0 & \frac{\partial N_{2}}{\partial x} & 0 & \frac{\partial N_{3}}{\partial x} & 0 \\ 0 & \frac{\partial N_{1}}{\partial y} & 0 & \frac{\partial N_{2}}{\partial y} & 0 & \frac{\partial N_{3}}{\partial y} \\ \frac{\partial N_{1}}{\partial y} & \frac{\partial N_{1}}{\partial x} & \frac{\partial N_{2}}{\partial y} & \frac{\partial N_{2}}{\partial x} & \frac{\partial N_{3}}{\partial y} & \frac{\partial N_{3}}{\partial x} \end{matrix}]$

3. 形函数的推导

为了确定B矩阵，我们必须首先找到形函数N_i的显式形式。对于3节点三角形单元，我们假设最简单的可能位移场：一个线性多项式。

系数A₁、A₂和A₃是常数。为了将它们与物理节点位移联系起来，我们在三个节点处分别施加该方程：

${\begin{matrix} u_{1} \\ u_{2} \\ u_{3} \end{matrix}} = [\begin{matrix} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix}] {\begin{matrix} A_{1} \\ A_{2} \\ A_{3} \end{matrix}}$

可以求解该方程组，用节点位移u_i和节点坐标表示系数A_i。

其中 $A = \frac{1}{2} | \begin{matrix} 1 & x_{1} & y_{1} \\ 1 & x_{2} & y_{2} \\ 1 & x_{3} & y_{3} \end{matrix} |$ 为三角形的面积。

将其代回线性多项式并重新排列，可得到形函数的最终形式：

由线性多项式场推导出的形函数形式为： $N_{i} (x, y) = \frac{1}{2 A} (a_{i} + b_{i} x + c_{i} y)$ 其中A是单元面积，而a_i、b_i和c_i是由节点坐标组成的系数。

注意，N_i(x, y) 是一个平面方程，在节点i处值为1，在另外两个节点处值为零。

4. 常应变B矩阵与单元刚度矩阵

这些形函数的导数为： $\frac{\partial N_{i}}{\partial x} = \frac{b_{i}}{2 A}, \frac{\partial N_{i}}{\partial y} = \frac{c_{i}}{2 A}$

线性形函数的一个重要结果是其导数为常数。由于 [B] 矩阵完全由这些导数组成，因此该单元的 [B] 矩阵也是常数。由关系式 {ϵ} = [B]{q} 可知，这意味着应变 {ϵ} 在整个单元内是均匀的，这就是它被称为常应变三角形（CST）的原因。

单元刚度矩阵K通过对单元域Ω进行积分求得： $𝐊 = \int_{Ω} 𝐁^{𝖳} 𝐄 𝐁 d V$

对于厚度恒为t的单元，该体积分变为面积分： $𝐊 = \int_{A} 𝐁^{𝖳} 𝐄 𝐁 t d A$

由于B、E和t都是常数，可以提到积分号外面，积分简化为面积A。由此得到单元刚度矩阵的最终表达式：

$𝐊 = (𝐁^{𝖳} 𝐄 𝐁) t A$

5. 单元间连续性与高阶单元

关于位移的连续性存在一个问题。我们知道位移在节点处是连续的，但它在节点之间的边上是否连续？

答案是肯定的。由于位移场是线性的，沿任何边的变化是两个角节点值之间的线性插值。共享该边的相邻单元，其位移场将由完全相同的线性插值描述。这确保了单元之间的位移兼容性，该条件称为 $C^{0}$ 连续性。

完备性条件确保单元能够表示刚体运动，它得到满足，如下性质所示：

$\sum N_{i} = 1$

CST单元的精度受限于其常应变假设。为了更准确地模拟具有变化应变场的问题，使用了高阶单元。这是通过向单元添加节点并对形函数使用更高阶多项式来实现的。例如，添加边中点节点可得到二次位移场（6节点三角形），从而产生线性变化的应变场。

这些高阶单元仍然只具有 $C^{0}$ 连续性。要实现更高的连续性（例如 $C^{1}$ ，导数也连续），则需要更复杂的单元列式。