代码说明 三维弹塑性有限元分析程序。采用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=1331n_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=3993U = 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; endend% ======================== 后处理 ========================% 绘制变形图（放大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 endend% 八节点六面体单元的三线性形函数及其导数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/∂ζ endend% 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'); endend% 计算应力张量的偏应力部分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 来源：有限元先生
