开源：三维弹塑性有限元分析程序

代码说明

  三维弹塑性有限元分析程序。采用C3D8单元，Von Mises屈服准则，关联流动法则，各向同性运动硬化法则，位移加载模式，映射返回算法，牛顿迭代。
  代码可以自动给长方体划分六面体网格，用户可以自定义长方体的三个尺寸，以及三个方向的单元尺寸，代码会自动生成节点编号和节点坐标信息。
  代码可以计算C3D8单元，如形函数、应变-位移矩阵、刚度矩阵等等。采用高斯积分方法计算刚度矩阵。
  弹塑性部分，采用Von Mises屈服准则，关联流动法则，各向同性运动硬化法则，应力试探和矫正采用映射返回算法。
  代码最后将计算的结果输出为Paraview可读取的Vtk文件，用于后续的可视化处理。

全部程序

% 三维弹塑性有限元分析程序 (立方体模型)
clear; clc; close all;  % 清除工作区、命令窗口和图形窗口
CoreNum=16; % 设置并行计算使用的CPU核心数（利用多核加速计算）
if isempty(gcp('nocreate')) % 检查并行计算池是否存在，不存在则创建
    parpool('local',CoreNum); % 创建包含16个工作线程的本地并行计算池
end
% ========================== 前处理 ==========================
% 材料参数 (基于J2塑性理论)
E = 200e3;         % 弹性模量 (MPa) - 胡克定律基础参数，公式：σ = Eε
nu = 0.3;          % 泊松比 - 横向应变与纵向应变比值，公式：ν = -ε横向/ε纵向
sigma_y0 = 200;    % 初始屈服应力 (MPa) - von Mises屈服准则临界值，σ_eq ≥ σ_y0时进入塑性
H = 0.1 * E;       % 硬化模量 (MPa) - 线性硬化模型参数，公式：σ_y = σ_y0 + Hε_p^eq

% 几何参数 (定义10×10×10 mm立方体模型)
Lx = 1; Ly = 1; Lz = 1;  % 模型尺寸 (mm)
nx = 10; ny = 10; nz = 10;  % 各方向单元数量，总单元数=10×10×10=1000

% 加载参数 (准静态位移控制加载)
total_disp = 1;   % 总位移量 (mm) - 模型右端沿X方向的总位移
n_steps = 10;       % 加载步数 - 控制位移增量为total_disp/n_steps

% 高斯积分设置 (2×2×2积分点，用于数值积分)
gauss_points = [-1/sqrt(3), 1/sqrt(3)];  % 一维高斯点坐标（对应ξ, η, ζ方向）
gauss_weights = [1, 1];                  % 高斯积分权重，用于加权求积

% 生成结构化六面体网格
[nodes, elements] = generate_mesh(Lx, Ly, Lz, nx, ny, nz);
n_nodes = size(nodes, 1);      % 节点总数：(nx+1)(ny+1)(nz+1)=11×11×11=1331
n_elements = size(elements, 1); % 单元总数：nx×ny×nz=1000

% 有限元求解变量初始化
dof_per_node = 3;              % 每个节点3个平动自由度（X, Y, Z方向）
total_dof = n_nodes * dof_per_node;  % 系统总自由度：1331×3=3993
U = zeros(total_dof, 1);       % 全局位移向量（待求解未知量）
F_ext = zeros(total_dof, 1);   % 外力向量（本例无体力，仅位移加载）

% 积分点历史变量（记录塑性状态）
n_gauss = 8;  % 每个单元8个高斯积分点（2×2×2布局）
plastic_strain = zeros(6, n_gauss, n_elements); % 塑性应变张量（Voigt六分量：[ε_xx, ε_yy, ε_zz, γ_xy, γ_yz, γ_zx]）
eq_plastic_strain = zeros(n_gauss, n_elements); % 等效塑性应变（硬化参数，ε_p^eq）

% ======================== 边界条件 ========================
% 固定X=0的面（全自由度固定）- Dirichlet边界条件
fixed_nodes = find(nodes(:,1) == 0);  % 查找X坐标为0的所有节点
% 将节点编号转换为自由度编号（每个节点3个自由度：X(3i-2), Y(3i-1), Z(3i)）
fixed_dof = sort([3*fixed_nodes-2; 3*fixed_nodes-1; 3*fixed_nodes]);

% 在X=Lx的面施加X方向位移- 位移控制边界条件
loaded_nodes = find(nodes(:,1) == Lx);  % 查找X坐标为Lx的所有节点
loaded_dof = 3*loaded_nodes - 2;        % X方向自由度编号（每个节点的第1个自由度）

% ====================== 主分析循环（增量加载） ======================
for step = 1:n_steps
    fprintf('加载步 %d/%d\n', step, n_steps);
    
    % 施加位移增量（控制加载速率）
    disp_inc = total_disp / n_steps;  % 每步位移增量
    U(loaded_dof) = U(loaded_dof) + disp_inc;  % 更新加载端X方向位移
    
    % 牛顿-拉夫逊迭代（求解非线性平衡方程 F_int(U) = F_ext）
    tol = 1e-6;       % 收敛容差：残余力范数小于1e-6 MPa
    max_iter = 20;    % 最大迭代次数，防止不收敛
    residual_norm = 1; % 初始残余力范数（用于循环条件）
    iter = 0;         % 迭代计数
    
    while residual_norm > tol && iter < max_iter
        iter = iter + 1;
        
        % 初始化全局刚度矩阵（稀疏存储）和内力向量
        K_global = sparse(total_dof, total_dof);  % 切线刚度矩阵（稀疏矩阵节省内存）
        F_int = zeros(total_dof, 1);              % 内力向量（初始为零）
        
        % 遍历所有单元，组装系统矩阵
        for el = 1:n_elements
            el_nodes = elements(el, :);  % 获取当前单元8个节点编号
            
            % 映射单元节点到自由度（每个节点3个自由度，共24个自由度/单元）
            el_dof = zeros(1, 24);
            for i = 1:8
                el_dof(3*i-2:3*i) = [3*el_nodes(i)-2, 3*el_nodes(i)-1, 3*el_nodes(i)];
            end
            
            % 提取单元位移向量（24×1）
            U_el = U(el_dof);
            
            % 初始化单元刚度矩阵和内力向量
            K_el = zeros(24, 24);  % 单元切线刚度矩阵
            F_el = zeros(24, 1);   % 单元内力向量
            
            % 高斯积分循环（2×2×2积分点）
            for ix = 1:2
                xi = gauss_points(ix);  % 自然坐标ξ
                for iy = 1:2
                    eta = gauss_points(iy);  % 自然坐标η
                    for iz = 1:2
                        zeta = gauss_points(iz);  % 自然坐标ζ
                        gpt = (ix-1)*4 + (iy-1)*2 + iz;  % 积分点索引（1-8）
                        
                        % 获取形函数及其对自然坐标的导数（八节点六面体等参单元）
                        [N, dN_dxi, dN_deta, dN_dzeta] = hex8_shape(xi, eta, zeta);
                        
                        % 计算雅可比矩阵（自然坐标→物理坐标的映射）
                        J = zeros(3,3);
                        for i = 1:8
                            x = nodes(el_nodes(i), 1); % 节点X坐标
                            y = nodes(el_nodes(i), 2); % 节点Y坐标
                            z = nodes(el_nodes(i), 3); % 节点Z坐标
                            % 雅可比矩阵元素：J = ∂(x,y,z)/∂(ξ,η,ζ)
                            J(1,1) = J(1,1) + dN_dxi(i)*x;
                            J(1,2) = J(1,2) + dN_dxi(i)*y;
                            J(1,3) = J(1,3) + dN_dxi(i)*z;
                            J(2,1) = J(2,1) + dN_deta(i)*x;
                            J(2,2) = J(2,2) + dN_deta(i)*y;
                            J(2,3) = J(2,3) + dN_deta(i)*z;
                            J(3,1) = J(3,1) + dN_dzeta(i)*x;
                            J(3,2) = J(3,2) + dN_dzeta(i)*y;
                            J(3,3) = J(3,3) + dN_dzeta(i)*z;
                        end
                        detJ = det(J);  % 雅可比行列式（体积缩放因子，detJ>0保证映射有效性）
                        invJ = inv(J);  % 雅可比逆矩阵（用于导数转换）
                        
                        % 构造应变-位移矩阵B（6×24，Voigt记号）
                        B = zeros(6, 24);
                        for i = 1:8
                            % 形函数对物理坐标的导数：∇N = J⁻¹·∂N/∂(ξ,η,ζ)
                            dN_dx = invJ * [dN_dxi(i); dN_deta(i); dN_dzeta(i)]; %#ok<MINV>
                            idx = 3*i-2;  % 节点自由度起始索引（如节点1对应自由度1,2,3）
                            % 填充B矩阵（对应Voigt应变分量）
                            B(1, idx) = dN_dx(1);   % ε_xx = ∂u/∂x
                            B(2, idx+1) = dN_dx(2); % ε_yy = ∂v/∂y
                            B(3, idx+2) = dN_dx(3); % ε_zz = ∂w/∂z
                            B(4, idx) = dN_dx(2);    % γ_xy = ∂v/∂x + ∂u/∂y（Voigt记号中γ_xy=2ε_xy）
                            B(4, idx+1) = dN_dx(1);
                            B(5, idx+1) = dN_dx(3);  % γ_yz = ∂w/∂y + ∂v/∂z
                            B(5, idx+2) = dN_dx(2);
                            B(6, idx) = dN_dx(3);    % γ_zx = ∂u/∂z + ∂w/∂x
                            B(6, idx+2) = dN_dx(1);
                        end
                        
                        % 计算应变增量：Δε = B·ΔU（当前迭代的位移增量）
                        delta_epsilon = B * U_el;
                        
                        % 提取上一增量步的历史变量
                        alpha_prev = eq_plastic_strain(gpt, el);  % 等效塑性应变ε_p^eq（旧值）
                        ep_prev = plastic_strain(:, gpt, el);     % 塑性应变张量ε_p（旧值）
                        
                        % 本构积分：隐式径向返回算法（求解塑性应力状态）
                        [stress, D, alpha_new, ep_new] = ...
                            radial_return(E, nu, sigma_y0, H, delta_epsilon, alpha_prev, ep_prev);
                        
                        % 更新历史变量（用于下次迭代或增量步）
                        plastic_strain(:, gpt, el) = ep_new;  % 存储新的塑性应变张量
                        eq_plastic_strain(gpt, el) = alpha_new;  % 存储新的等效塑性应变
                        
                        % 计算积分点贡献（高斯积分权重×体积元）
                        weight = gauss_weights(ix)*gauss_weights(iy)*gauss_weights(iz)*detJ;
                        % 单元刚度矩阵积分：K_e = ∫B^T·D·B dV（通过高斯积分近似）
                        K_el = K_el + (B' * D * B) * weight;
                        % 单元内力积分：F_int_e = ∫B^T·σ dV
                        F_el = F_el + (B' * stress) * weight;
                    end
                end
            end
            
            % 组装到全局系统（有限元核心步骤：单元矩阵→全局矩阵）
            K_global(el_dof, el_dof) = K_global(el_dof, el_dof) + K_el; %#ok<*SPRIX>
            F_int(el_dof) = F_int(el_dof) + F_el;
        end
        
        % 计算残余力：R = F_ext - F_int（外力-内力，非线性方程残差）
        residual = F_ext - F_int;
        
        % 应用边界条件：固定自由度置零，保留自由自由度
        free_dof = setdiff(1:total_dof, fixed_dof);  % 未约束的自由度索引
        K_red = K_global(free_dof, free_dof);         % 缩减的刚度矩阵（仅自由自由度）
        residual_red = residual(free_dof);             % 缩减的残余力向量
        
        % 检查收敛性（残余力范数是否小于容差）
        residual_norm = norm(residual_red);
        fprintf('  迭代 %d: 残差 = %.4e\n', iter, residual_norm);
        
        % 求解位移增量：ΔU = K_red⁻¹·R_red（线性方程组求解，稀疏矩阵法）
        delta_U = zeros(total_dof, 1);
        delta_U(free_dof) = K_red \ residual_red;  % 使用反斜杠运算符求解线性系统
        
        % 更新位移向量（牛顿-拉夫逊迭代核心：U_{k+1} = U_k + ΔU）
        U = U + delta_U;
    end
end

% ======================== 后处理 ========================
% 绘制变形图（放大10倍便于观察小变形）
figure;
plot_deformed(nodes, elements, U);
title('变形图');
xlabel('X'); ylabel('Y'); zlabel('Z');

% 绘制等效应力云图（基于节点平均应力）
figure;
node_stress=plot_stress(nodes, elements, U, plastic_strain);
title('等效应力云图');
xlabel('X'); ylabel('Y'); zlabel('Z');

output_vtk('result.vtk', nodes, elements, U, node_stress);

fprintf('分析完成!\n');


% 生成结构化六面体网格
function [nodes, elements] = generate_mesh(Lx, Ly, Lz, nx, ny, nz)
    dx = Lx/nx; dy = Ly/ny; dz = Lz/nz;  % 单元尺寸（物理坐标间隔）
    n_nodes = (nx+1)*(ny+1)*(nz+1);      % 节点总数（每个方向节点数=单元数+1）
    nodes = zeros(n_nodes, 3);            % 节点坐标矩阵（X,Y,Z）
    elements = zeros(nx*ny*nz, 8);        % 单元拓扑矩阵（每个单元8个节点）
    
    % 生成节点坐标（规则网格，按Z-Y-X顺序排列）
    node_id = 1;
    for k = 0:nz  % Z方向循环
        for j = 0:ny  % Y方向循环
            for i = 0:nx  % X方向循环
                nodes(node_id, :) = [i*dx, j*dy, k*dz];  % 节点坐标=单元尺寸×索引
                node_id = node_id + 1;
            end
        end
    end
    
    % 生成单元连接（八节点六面体，右手法则排序）
    el_id = 1;
    for k = 1:nz  % Z方向单元循环
        for j = 1:ny  % Y方向单元循环
            for i = 1:nx  % X方向单元循环
                % 计算底层4节点和顶层4节点的全局编号
                n1 = (k-1)*(nx+1)*(ny+1) + (j-1)*(nx+1) + i;       % 底层前左
                n2 = n1 + 1;                                        % 底层前右
                n4 = (k-1)*(nx+1)*(ny+1) + j*(nx+1) + i;           % 底层后左
                n3 = n4 + 1;                                        % 底层后右
                n5 = n1 + (nx+1)*(ny+1);                           % 顶层前左
                n6 = n2 + (nx+1)*(ny+1);                           % 顶层前右
                n7 = n3 + (nx+1)*(ny+1);                           % 顶层后右
                n8 = n4 + (nx+1)*(ny+1);                           % 顶层后左
                % 存储单元节点（右手法则：底层逆时针→顶层逆时针，保证法向量一致性）
                elements(el_id, :) = [n1, n2, n3, n4, n5, n6, n7, n8];
                el_id = el_id + 1;
            end
        end
    end
end

% 八节点六面体单元的三线性形函数及其导数
function [N, dN_dxi, dN_deta, dN_dzeta] = hex8_shape(xi, eta, zeta)
    N = zeros(8,1);          % 形函数值（N1到N8）
    dN_dxi = zeros(8,1);     % 形函数对ξ的导数
    dN_deta = zeros(8,1);    % 形函数对η的导数
    dN_dzeta = zeros(8,1);   % 形函数对ζ的导数
    
    % 参考坐标系中节点坐标（标准立方体顶点，ξ,η,ζ∈{-1,1}）
    xi_vec = [-1, 1, 1, -1, -1, 1, 1, -1];  % 节点1-8的ξ坐标
    eta_vec = [-1, -1, 1, 1, -1, -1, 1, 1];  % 节点η坐标
    zeta_vec = [-1, -1, -1, -1, 1, 1, 1, 1];  % 节点ζ坐标
    
    % 计算三线性形函数及其导数（N_i = (1+ξξ_i)(1+ηη_i)(1+ζζ_i)/8）
    for i = 1:8
        N(i) = (1/8)*(1+xi*xi_vec(i))*(1+eta*eta_vec(i))*(1+zeta*zeta_vec(i));
        dN_dxi(i) = (1/8)*xi_vec(i)*(1+eta*eta_vec(i))*(1+zeta*zeta_vec(i));  % ∂N_i/∂ξ
        dN_deta(i) = (1/8)*(1+xi*xi_vec(i))*eta_vec(i)*(1+zeta*zeta_vec(i));  % ∂N_i/∂η
        dN_dzeta(i) = (1/8)*(1+xi*xi_vec(i))*(1+eta*eta_vec(i))*zeta_vec(i);  % ∂N_i/∂ζ
    end
end

% J2塑性理论隐式本构积分（径向返回算法）
function [stress, D, alpha_new, ep_new] = radial_return(E, nu, sigma_y0, H, delta_epsilon, alpha_prev, ep_prev)
    % 计算弹性常数
    mu = E/(2*(1+nu));       % 剪切模量，公式：μ = E/(2(1+ν))
    lambda = E*nu/((1+nu)*(1-2*nu));  % Lamé常数，公式：λ = Eν/[(1+ν)(1-2ν)]
    
    % 弹性刚度矩阵（Voigt表示，三维弹性本构：σ = D_el·ε）
    D_el = [lambda+2*mu, lambda, lambda, 0, 0, 0;
            lambda, lambda+2*mu, lambda, 0, 0, 0;
            lambda, lambda, lambda+2*mu, 0, 0, 0;
            0,0,0,mu,0,0;
            0,0,0,0,mu,0;
            0,0,0,0,0,mu];
    
    % 弹性预测步：计算试探应力（假设当前增量步为弹性）
    stress_tr = D_el * delta_epsilon;  % 试探应力，σ_tr = D_el·Δε
    s_tr = deviatoric(stress_tr);      % 试探偏应力，s_tr = σ_tr - (trσ_tr/3)I
    sigma_tr_eq = sqrt(3/2 * (s_tr'*s_tr));  % von Mises等效应力，σ_eq = sqrt(3/2 s_tr:s_tr)
    
    % 当前屈服应力（考虑线性硬化）
    sigma_y = sigma_y0 + H*alpha_prev;  % σ_y = σ_y0 + Hε_p^eq（α_prev=ε_p^eq旧值）
    
    if sigma_tr_eq <= sigma_y  % 未屈服，弹性状态
        stress = stress_tr;          % 应力不变
        D = D_el;                   % 切线刚度为弹性刚度
        alpha_new = alpha_prev;      % 等效塑性应变不变
        ep_new = ep_prev;            % 塑性应变不变
    else  % 屈服，塑性状态
        n_tr = s_tr / sigma_tr_eq;  % 流动方向单位向量，n = s_tr/σ_tr_eq
        delta_gamma = (sigma_tr_eq - sigma_y) / (3*mu + H);  % 塑性乘子增量，Δγ = (σ_tr_eq - σ_y)/(3μ + H)
        
        % 应力修正：径向返回至屈服面，σ = σ_tr - 2μΔγ n
        stress = stress_tr - 2*mu*delta_gamma*n_tr;
        
        % 塑性应变更新，Δε_p = sqrt(3/2)Δγ n
        delta_ep = delta_gamma * sqrt(3/2) * n_tr;
        ep_new = ep_prev + delta_ep;  % 累积塑性应变
        alpha_new = alpha_prev + delta_gamma;  % 累积等效塑性应变
        
        % 计算一致切线模量（保证牛顿法二次收敛）
        theta = 1 - 3*mu*delta_gamma/sigma_tr_eq;  % 修正参数
        theta_bar = 1 / (1 + H/(3*mu)) - (1 - theta);  % 偏量部分修正参数
        I_dev = diag([1,1,1,0.5,0.5,0.5]) - [1;1;1;0;0;0]*[1,1,1,0,0,0]/3;  % 偏量投影算子，I_dev = I - (1/3)1⊗1
        
        % 切线刚度矩阵，D = D_el - 2μ[θ_bar(n⊗n) + θ(I_dev - n⊗n)]
        D = D_el - 2*mu*theta_bar*(n_tr*n_tr') - 2*mu*theta*(I_dev - n_tr*n_tr');
    end
end

% 计算应力张量的偏应力部分
function s = deviatoric(sigma)
    p = (sigma(1) + sigma(2) + sigma(3)) / 3;  % 静水压力，公式：p = (σ_xx + σ_yy + σ_zz)/3
    s = sigma;  % 复 制应力张量
    s(1:3) = s(1:3) - p;  % 正应力分量减去静水压力，得到偏应力分量
    % 剪应力分量不变（剪应力天然属于偏量部分）
end

% 绘制变形后的三维模型
function plot_deformed(nodes, elements, U)
    scale = 10;  % 变形放大系数，便于观察小变形（实际位移×scale）
    
    % 提取各方向位移分量（按自由度顺序）
    Ux = U(1:3:end);  % X方向位移
    Uy = U(2:3:end);  % Y方向位移
    Uz = U(3:3:end);  % Z方向位移
    
    % 计算变形后节点坐标：原始坐标 + 放大后位移
    deformed_nodes = nodes + scale * [Ux, Uy, Uz];
    
    % 绘制六面体单元（每个单元6个面）
    for el = 1:size(elements,1)
        el_nodes = elements(el,:);  % 单元节点编号
        X = deformed_nodes(el_nodes,1);  % 变形后X坐标
        Y = deformed_nodes(el_nodes,2);  % 变形后Y坐标
        Z = deformed_nodes(el_nodes,3);  % 变形后Z坐标
        
        % 六面体的6个面（右手法则定义顶点顺序，保证法向量一致性）
        faces = [1,2,3,4; 5,6,7,8; 1,2,6,5; 2,3,7,6; 3,4,8,7; 4,1,5,8];
        for f = 1:6
            % 绘制半透明面片（FaceAlpha=0.1）
            patch(X(faces(f,:)), Y(faces(f,:)), Z(faces(f,:)), 'b', 'FaceAlpha', 0.1);
        end
    end
    
    axis equal;  % 等比例显示坐标轴，避免图形变形
    grid on;     % 显示网格线辅助观察
    view(3);     % 切换至三维视角
end

% 绘制节点平均等效应力云图，并返回节点等效应力
function node_stress = plot_stress(nodes, elements, U, plastic_strain)
    n_nodes = size(nodes, 1);  % 节点总数
    node_stress = zeros(n_nodes, 1);  % 节点等效应力（初始为零）
    node_count = zeros(n_nodes, 1);   % 每个节点的积分点贡献计数

    % 遍历所有单元和积分点，插值应力到节点（单元平均法）
    for el = 1:size(elements,1)
        for gpt = 1:8
            % 取该积分点的塑性应变张量
            ep = plastic_strain(:,gpt,el);
            % 计算von Mises等效应力
            seq = sqrt( (ep(1)-ep(2))^2 + (ep(2)-ep(3))^2 + (ep(3)-ep(1))^2 ...
                + 6*(ep(4)^2 + ep(5)^2 + ep(6)^2) ) / sqrt(2);
            % 将积分点应力累加到所属节点（8个节点/单元）
            for i = 1:8
                n = elements(el,i);  % 节点编号
                node_stress(n) = node_stress(n) + seq;
                node_count(n) = node_count(n) + 1;
            end
        end
    end

    % 计算节点平均应力（总应力/贡献积分点数）
    node_stress = node_stress ./ max(node_count,1);

    % 计算变形后坐标（放大显示位移场）
    Ux = U(1:3:end);
    Uy = U(2:3:end);
    Uz = U(3:3:end);
    deformed_nodes = nodes + 10 * [Ux, Uy, Uz];

    % 绘制应力云图（散点图，颜色映射应力大小，点大小50）
    scatter3(deformed_nodes(:,1), deformed_nodes(:,2), deformed_nodes(:,3), ...
        50, node_stress, 'filled');
    colorbar;  % 显示颜色标尺
    title('等效应力分布');
    axis equal;  % 等比例坐标轴
    grid on;     % 显示网格
    view(3);     % 三维视角
end


% VTK导出函数（支持八节点六面体单元和节点标量）
function output_vtk(filename, nodes, elements, U, node_stress)
    % 变形后节点坐标（放大10倍）
    scale = 10;
    Ux = U(1:3:end);
    Uy = U(2:3:end);
    Uz = U(3:3:end);
    deformed_nodes = nodes + scale * [Ux, Uy, Uz];

    n_nodes = size(nodes,1);
    n_elements = size(elements,1);

    fid = fopen(filename, 'w');
    fprintf(fid, '# vtk DataFile Version 3.0\n');
    fprintf(fid, 'ElasPlastic result\n');
    fprintf(fid, 'ASCII\n');
    fprintf(fid, 'DATASET UNSTRUCTURED_GRID\n');
    % 节点坐标
    fprintf(fid, 'POINTS %d float\n', n_nodes);
    for i = 1:n_nodes
        fprintf(fid, '%f %f %f\n', deformed_nodes(i,1), deformed_nodes(i,2), deformed_nodes(i,3));
    end
    % 单元连接
    fprintf(fid, 'CELLS %d %d\n', n_elements, n_elements*9);
    for i = 1:n_elements
        % VTK八节点六面体单元类型，节点编号从0开始
        fprintf(fid, '8 %d %d %d %d %d %d %d %d\n', elements(i,:)-1);
    end
    % 单元类型
    fprintf(fid, 'CELL_TYPES %d\n', n_elements);
    for i = 1:n_elements
        fprintf(fid, '%d\n', 12); % 12=VTK_HEXAHEDRON
    end
    % 节点标量（等效应力）
    fprintf(fid, 'POINT_DATA %d\n', n_nodes);
    fprintf(fid, 'SCALARS von_mises float 1\n');
    fprintf(fid, 'LOOKUP_TABLE default\n');
    for i = 1:n_nodes
        fprintf(fid, '%f\n', node_stress(i));
    end
    fclose(fid);
    fprintf('VTK文件已导出: %s\n', filename);
end