matlab实现三维四面体单元的有限元解法

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：单元分析

对每个单元：

1. 计算形函数及其导数
2. 计算B矩阵
3. 计算单元刚度矩阵 $\mathbf{k}^e$
4. 计算等效节点力 $\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. 优缺点

优点：

1. 几何适应性强，可离散复杂三维区域
2. 网格生成相对容易
3. 自动满足收敛条件

缺点：

1. 计算精度较低（常应变）
2. 单元数量通常较多
3. 可能产生剪切锁死

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. 实际注意事项

1. 网格质量：确保四面体形状良好（避免小角度）
2. 应力平滑：常应力单元需要后处理平滑
3. 收敛性：满足位移协调性和常应变条件
4. 并行计算：大规模问题需要并行求解