一、有限元分析技术体系概述

有限元分析(FEA)作为数值计算领域的核心方法，通过离散化连续体为有限个单元，建立代数方程组求解复杂工程问题。MATLAB凭借其强大的矩阵运算能力和灵活的编程接口，成为实现自定义有限元算法的理想平台。相较于商业软件，MATLAB方案具有算法透明、可定制性强、适合教学研究等优势。

1.1 核心实现路径

典型有限元分析流程包含以下关键步骤：

几何建模与网格划分
单元类型选择与形函数构造
刚度矩阵组装与边界条件施加
线性方程组求解
结果可视化与后处理

二、一维结构分析实现

2.1 平面纯弯梁分析

采用欧拉-伯努利梁理论，建立4节点Hermitian单元模型：

function [K] = euler_beam_stiffness(E,I,L)
    % E:弹性模量 I:截面惯性矩 L:单元长度
    k = E*I/L^3;
    K = k*[12 6*L -12 6*L;
           6*L 4*L^2 -6*L 2*L^2;
           -12 -6*L 12 -6*L;
           6*L 2*L^2 -6*L 4*L^2];
end

通过节点位移求解弯矩分布，验证圣维南原理在局部载荷作用下的适用性。

2.2 铁木辛柯梁单元

考虑剪切变形影响，引入剪切修正系数κ：

function [K] = timoshenko_beam(E,G,A,I,L,kappa)
    % G:剪切模量 A:截面积 kappa:剪切修正系数
    phi = 12*E*I/(G*A*kappa*L^2);
    K_bend = euler_beam_stiffness(E,I,L);
    K_shear = G*A*kappa/L*[1 -1; -1 1];
    % 组装总刚度矩阵（需考虑耦合项）
    ...
end

对比欧拉梁结果，分析高跨比对剪切变形敏感度的影响。

三、二维板壳结构分析

3.1 矩形薄板小挠度理论

基于Kirchhoff假设建立Mindlin板单元，采用双线性形函数：

function [N,dNdx,dNdy] = thin_plate_shape(xi,eta)
    % 自然坐标系下的形函数及其导数
    N = 0.25*[(1-xi)*(1-eta); (1+xi)*(1-eta); 
              (1+xi)*(1+eta); (1-xi)*(1+eta)];
    dNdx = [-0.25*(1-eta); 0.25*(1-eta); 
             0.25*(1+eta); -0.25*(1+eta)];
    dNdy = [-0.25*(1-xi); -0.25*(1+xi); 
             0.25*(1+xi); 0.25*(1-xi)];
end

通过高斯积分计算单元刚度矩阵，验证四边简支板在均布载荷下的解析解。

3.2 Mindlin厚板理论

引入横向剪切变形自由度，修正刚度矩阵计算：

function [K] = mindlin_plate_stiffness(D,k,h,A)
    % D:弯曲刚度 k:剪切修正因子 h:板厚 A:单元面积
    % 需计算5个自由度对应的子矩阵
    K_bending = ... % 弯曲项
    K_shear = ...  % 剪切项
    K_coupling = ... % 耦合项
    K = K_bending + K_shear + K_coupling;
end

分析厚跨比对结果收敛性的影响，确定合理积分点数量。

四、三维实体结构分析

4.1 四面体单元实现

采用线性四面体单元进行悬臂梁分析：

function [Ke] = tetra_stiffness(E,nu,nodes)
    % nodes:4个节点的坐标矩阵(3x4)
    % 计算雅可比矩阵及其行列式
    J = (nodes(:,2)-nodes(:,1))*(nodes(:,3)-nodes(:,1))' - ...
        (nodes(:,2)-nodes(:,1))*(nodes(:,4)-nodes(:,1))';
    detJ = det(J);
    % 构造B矩阵（需处理3D应变-位移关系）
    B = zeros(6,12);
    % ...（具体推导过程省略）
    % 计算弹性矩阵D
    D = E/((1+nu)*(1-2*nu))*[...];
    Ke = B'*D*B*abs(detJ)/6;
end

通过网格加密验证收敛性，对比解析解误差。

4.2 六面体单元优化

采用20节点高阶单元提升计算精度：

function [N,dNdx] = hex20_shape(xi,eta,zeta)
    % 三维Serendipity族形函数
    N = zeros(20,1);
    % 角节点项
    N(1) = 0.125*(1-xi)*(1-eta)*(1-zeta)*(-xi-eta-zeta-2);
    % ...（其他节点项）
    % 计算形函数对自然坐标的导数
    dN = zeros(20,3);
    % ...（导数计算过程）
    % 转换为物理坐标导数（需计算雅可比逆矩阵）
    % ...
end

分析减缩积分与完全积分对体积锁死现象的影响。

五、结构动力学分析

5.1 模态分析实现

通过特征值问题求解结构固有频率：

function [freq,modes] = modal_analysis(K,M,n)
    % K:总刚度矩阵 M:总质量矩阵 n:需求模态数
    [V,D] = eig(K,M);
    [freq,idx] = sort(sqrt(diag(D))/(2*pi));
    modes = V(:,idx(1:n));
end

验证框架结构前几阶振型与理论解的一致性。

5.2 时程分析方法

采用Newmark-β法进行地震响应分析：

function [u,v,a] = newmark_time_history(M,C,K,p,dt,beta,gamma)
    % 初始化位移、速度、加速度
    u = zeros(size(p,1),size(p,2));
    v = zeros(size(p,1),size(p,2));
    a = zeros(size(p,1),size(p,2));
    % 有效刚度矩阵
    K_eff = K + (gamma/(beta*dt))*C + (1/(beta*dt^2))*M;
    for i = 2:size(p,2)
        % 预测步
        u_pred = u(:,i-1) + dt*v(:,i-1) + 0.5*dt^2*(1-2*beta)*a(:,i-1);
        v_pred = v(:,i-1) + (1-gamma)*dt*a(:,i-1);
        % 修正步
        delta_p = p(:,i) - C*v_pred - M*(u_pred/(beta*dt^2) - v_pred/(beta*dt) - (1-2*beta)*a/(2*beta));
        delta_u = K_eff\delta_p;
        u(:,i) = u_pred + delta_u;
        v(:,i) = v_pred + gamma*delta_u/dt;
        a(:,i) = (u(:,i)-u(:,i-1)-dt*v(:,i-1))/(beta*dt^2);
    end
end

对比不同阻尼比对结构响应的影响。

六、结果可视化与后处理

6.1 不规则云图绘制

采用散点插值法生成连续应力云图：

function plot_stress_contour(nodes,elements,stress_values)
    % nodes:节点坐标矩阵 elements:单元连接矩阵
    % stress_values:节点应力值向量
    % 创建三角剖分
    tri = delaunay(nodes(:,1),nodes(:,2));
    % 绘制填充等高线
    tricontourf(tri,nodes(:,1),nodes(:,2),stress_values,20);
    colorbar;
    axis equal;
end

6.2 动画生成技术

通过循环绘制实现振动模态动画：

function animate_mode(nodes,elements,mode_shape,fps)
    figure;
    max_amp = max(abs(mode_shape));
    for t = linspace(0,2*pi,100)
        scaled_mode = sin(t)*mode_shape/max_amp;
        deformed_nodes = nodes + reshape(scaled_mode,[],2);
        clf;
        hold on;
        for e = 1:size(elements,1)
            elem_nodes = elements(e,:);
            plot(deformed_nodes(elem_nodes,1),...
                 deformed_nodes(elem_nodes,2),'b-');
        end
        axis equal;
        title(sprintf('Mode Shape Animation (%.1f FPS)',fps));
        drawnow;
        pause(1/fps);
    end
end

七、工程应用建议

模型验证：始终通过简单算例验证程序正确性
单元选择：根据结构特征选择合适单元类型（梁/板/壳/实体）
网格密度：在应力集中区域采用局部加密网格
计算效率：对于大规模问题，考虑使用稀疏矩阵存储格式
结果解读：关注相对值而非绝对值，重点分析分布规律

通过系统掌握上述方法体系，工程师可构建自主可控的有限元分析平台，既能满足教学研究需求，也可为复杂工程问题提供定制化解决方案。建议结合具体项目需求，逐步扩展分析维度（如热-力耦合、非线性材料等），构建完整的数值仿真能力体系。

MATLAB有限元分析全流程实践指南