1. 基本概念
1.1 四面体单元特点
- 最简单的三维实体单元
- 4个节点，每个节点3个自由度（位移u, v, w）
- 线性位移模式
- 常应变/常应力单元
2. 形函数和坐标变换
2.1 体积坐标（自然坐标）
对于四面体单元，使用体积坐标 $L_1, L_2, L_3, L_4$ ：
$L_i = \frac{V_i}{V}$
其中 V 是四面体体积， $V_i$ 是节点i对面所形成的小四面体体积。
关系： $L_1 + L_2 + L_3 + L_4 = 1$
2.2 线性形函数
对于4节点四面体单元：
$N_i = L_i \quad (i=1,2,3,4)$
2.3 坐标插值
物理坐标与体积坐标的关系：
$$x = \sum_{i=1}^4 N_i x_i, \quad y = \sum_{i=1}^4 N_i y_i, \quad z = \sum_{i=1}^4 N_i z_i$$
3. 位移插值
3.1 位移场表达式
$$u = \sum_{i=1}^4 N_i u_i, \quad v = \sum_{i=1}^4 N_i v_i, \quad w = \sum_{i=1}^4 N_i w_i$$
写成矩阵形式：
$\mathbf{u} = \begin{bmatrix} u \ v \ w \end{bmatrix} = \mathbf{N} \mathbf{a}^e$
其中：
$$\mathbf{N} = \begin{bmatrix} N_1 & 0 & 0 & N_2 & 0 & 0 & N_3 & 0 & 0 & N_4 & 0 & 0 \\ 0 & N_1 & 0 & 0 & N_2 & 0 & 0 & N_3 & 0 & 0 & N_4 & 0 \\ 0 & 0 & N_1 & 0 & 0 & N_2 & 0 & 0 & N_3 & 0 & 0 & N_4 \end{bmatrix}$$
$$\mathbf{a}^e = [u_1, v_1, w_1, u_2, v_2, w_2, u_3, v_3, w_3, u_4, v_4, w_4]^T$$
4. 应变-位移关系
4.1 几何方程
对于三维问题，应变向量：
$$\boldsymbol{\varepsilon} = [\varepsilon_x, \varepsilon_y, \varepsilon_z, \gamma_{xy}, \gamma_{yz}, \gamma_{zx}]^T$$
4.2 应变矩阵B
$\boldsymbol{\varepsilon} = \mathbf{B} \mathbf{a}^e$
其中 $\mathbf{B} = [\mathbf{B}_1, \mathbf{B}_2, \mathbf{B}_3, \mathbf{B}_4]$
每个子矩阵：
$$\mathbf{B}_i = \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}$$
4.3 形函数导数的计算
通过雅可比变换：
$$\begin{bmatrix} \frac{\partial N_i}{\partial x} \\ \frac{\partial N_i}{\partial y} \\ \frac{\partial N_i}{\partial z} \end{bmatrix} = \mathbf{J}^{-1} \begin{bmatrix} \frac{\partial N_i}{\partial L_1} \\ \frac{\partial N_i}{\partial L_2} \\ \frac{\partial N_i}{\partial L_3} \end{bmatrix}$$
雅可比矩阵：
$$\mathbf{J} = \begin{bmatrix} \frac{\partial x}{\partial L_1} & \frac{\partial y}{\partial L_1} & \frac{\partial z}{\partial L_1} \\ \frac{\partial x}{\partial L_2} & \frac{\partial y}{\partial L_2} & \frac{\partial z}{\partial L_2} \\ \frac{\partial x}{\partial L_3} & \frac{\partial y}{\partial L_3} & \frac{\partial z}{\partial L_3} \end{bmatrix} = \begin{bmatrix} x_1 & y_1 & z_1 \\ x_2 & y_2 & z_2 \\ x_3 & y_3 & z_3 \end{bmatrix} \begin{bmatrix} \frac{\partial N_1}{\partial L_1} & \cdots & \frac{\partial N_4}{\partial L_1} \\ \frac{\partial N_1}{\partial L_2} & \cdots & \frac{\partial N_4}{\partial L_2} \\ \frac{\partial N_1}{\partial L_3} & \cdots & \frac{\partial N_4}{\partial L_3} \end{bmatrix}$$
5. 单元刚度矩阵
5.1 刚度矩阵计算
$\mathbf{k}^e = \int_{V_e} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV$
其中 $\mathbf{D}$ 是弹性矩阵：
对于各向同性材料：
$$\mathbf{D} = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1 & \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & 1 & \frac{\nu}{1-\nu} & 0 & 0 & 0 \\ \frac{\nu}{1-\nu} & \frac{\nu}{1-\nu} & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2(1-\nu)} \end{bmatrix}$$
5.2 数值积分
对于常应变四面体，B矩阵是常数，因此：
$\mathbf{k}^e = \mathbf{B}^T \mathbf{D} \mathbf{B} V_e$
其中 $V_e$ 是单元体积。
四面体体积计算：
$$V_e = \frac{1}{6} \det \begin{vmatrix} 1 & x_1 & y_1 & z_1 \\ 1 & x_2 & y_2 & z_2 \\ 1 & x_3 & y_3 & z_3 \\ 1 & x_4 & y_4 & z_4 \end{vmatrix}$$
6. 等效节点力
6.1 体积力
$$\mathbf{f}_V^e = \int_{V_e} \mathbf{N}^T \mathbf{b} \, dV$$
其中 $\mathbf{b} = [b_x, b_y, b_z]^T$ 是体积力密度。
对于均匀分布：
$$\mathbf{f}_V^e = \frac{V_e}{4} [b_x, b_y, b_z, b_x, b_y, b_z, b_x, b_y, b_z, b_x, b_y, b_z]^T$$
6.2 表面力
对于作用在面上的分布力：
$$\mathbf{f}_S^e = \int_{A_f} \mathbf{N}^T \mathbf{t} \, dA$$
7. 求解步骤
步骤1：网格生成
- 将三维域离散为四面体网格
- 确保网格质量（避免过于扁平的四面体）
步骤2：单元分析
对每个单元：
- 计算形函数及其导数
- 计算B矩阵
- 计算单元刚度矩阵 $\mathbf{k}^e$
- 计算等效节点力 $\mathbf{f}^e$
步骤3：整体组装
$\mathbf{K} = \sum_e \mathbf{k}^e, \quad \mathbf{F} = \sum_e \mathbf{f}^e$
形成整体方程：
$\mathbf{K} \mathbf{a} = \mathbf{F}$
步骤4：边界条件处理
- 处理位移边界条件
- 处理力边界条件
步骤5：求解线性方程组
使用直接法（如LDLT分解）或迭代法（如PCG）求解。
8. 高阶四面体单元
10节点二次四面体
增加中间节点，位移模式为二次：
$u = \sum_{i=1}^{10} N_i u_i$
形函数：
- 角节点： $N_i = L_i(2L_i - 1)$
- 边中点： $N_5 = 4L_1 L_2 , N_6 = 4L_1 L_3$ , 等
9. 优缺点
优点：
- 几何适应性强，可离散复杂三维区域
- 网格生成相对容易
- 自动满足收敛条件
缺点：
- 计算精度较低（常应变）
- 单元数量通常较多
- 可能产生剪切锁死
10. 应用示例（MATLAB伪代码）
function [K, F] = TetrahedralFEM(nodes, elements, E, nu, force)
% nodes: N×3节点坐标
% elements: M×4单元连接
% E: 弹性模量
% nu: 泊松比
nNodes = size(nodes, 1);
nDOF = 3 * nNodes;
K = sparse(nDOF, nDOF);
F = zeros(nDOF, 1);
% D矩阵
D = E/(1+nu)/(1-2*nu) * [
1-nu, nu, nu, 0, 0, 0;
nu, 1-nu, nu, 0, 0, 0;
nu, nu, 1-nu, 0, 0, 0;
0, 0, 0, (1-2*nu)/2, 0, 0;
0, 0, 0, 0, (1-2*nu)/2, 0;
0, 0, 0, 0, 0, (1-2*nu)/2];
for e = 1:size(elements, 1)
% 提取单元节点
elemNodes = elements(e, :);
coords = nodes(elemNodes, :);
% 计算体积和B矩阵
[B, Ve] = computeBmatrix(coords);
% 单元刚度矩阵
ke = B' * D * B * Ve;
% 组装
dofs = zeros(12, 1);
for i = 1:4
dofs(3*i-2:3*i) = 3*elemNodes(i)-2:3*elemNodes(i);
end
K(dofs, dofs) = K(dofs, dofs) + ke;
end
% 施加边界条件和载荷
% ...
end
function [B, V] = computeBmatrix(coords)
% 计算四面体B矩阵和体积
x = coords(:,1); y = coords(:,2); z = coords(:,3);
% 计算体积
V = det([1, x(1), y(1), z(1);
1, x(2), y(2), z(2);
1, x(3), y(3), z(3);
1, x(4), y(4), z(4)]) / 6;
% 计算形函数导数
% ... 具体实现省略
B = zeros(6, 12);
% 填充B矩阵
end
参考代码 四结点四面体单元 www.youwenfan.com/contentalg/97774.html
11. 实际注意事项
- 网格质量：确保四面体形状良好（避免小角度）
- 应力平滑：常应力单元需要后处理平滑
- 收敛性：满足位移协调性和常应变条件
- 并行计算：大规模问题需要并行求解
