MATLAB 有限元程序：三节点三角形常应变单元

以下是一个三节点三角形常应变单元的 MATLAB 有限元程序：

clear all;
close all;
clc;

% 定义网格
nodes = [0,0; 1,0; 0,1; 0.5,0.5];
elements = [1,2,4; 1,4,3];

% 定义材料属性
E = 200e9; % 杨氏模量
nu = 0.3; % 泊松比
k = E/(3*(1-2*nu)); % 体积模量
mu = E/(2*(1+nu)); % 剪切模量

% 定义边界条件
u_fixed = [1,0; 2,0];

% 构建全局刚度矩阵和力向量
K = zeros(2*size(nodes,1));
F = zeros(2*size(nodes,1),1);
for e = 1:size(elements,1)
    % 获取单元节点
    n1 = elements(e,1);
    n2 = elements(e,2);
    n3 = elements(e,3);
    x1 = nodes(n1,1);
    y1 = nodes(n1,2);
    x2 = nodes(n2,1);
    y2 = nodes(n2,2);
    x3 = nodes(n3,1);
    y3 = nodes(n3,2);
    % 计算单元刚度矩阵
    B = [y2-y3, y3-y1, y1-y2; x3-x2, x1-x3, x2-x1];
    A = 0.5*det([1,x1,y1; 1,x2,y2; 1,x3,y3]);
    D = [2*mu+lambda, lambda, 0; lambda, 2*mu+lambda, 0; 0, 0, mu];
    Ke = B'*D*B*A;
    % 将单元刚度矩阵组装到全局刚度矩阵中
    K(2*n1-1:2*n1,2*n1-1:2*n1) = K(2*n1-1:2*n1,2*n1-1:2*n1) + Ke(1:2,1:2);
    K(2*n1-1:2*n1,2*n2-1:2*n2) = K(2*n1-1:2*n1,2*n2-1:2*n2) + Ke(1:2,3:4);
    K(2*n1-1:2*n1,2*n3-1:2*n3) = K(2*n1-1:2*n1,2*n3-1:2*n3) + Ke(1:2,5:6);
    K(2*n2-1:2*n2,2*n1-1:2*n1) = K(2*n2-1:2*n2,2*n1-1:2*n1) + Ke(3:4,1:2);
    K(2*n2-1:2*n2,2*n2-1:2*n2) = K(2*n2-1:2*n2,2*n2-1:2*n2) + Ke(3:4,3:4);
    K(2*n2-1:2*n2,2*n3-1:2*n3) = K(2*n2-1:2*n2,2*n3-1:2*n3) + Ke(3:4,5:6);
    K(2*n3-1:2*n3,2*n1-1:2*n1) = K(2*n3-1:2*n3,2*n1-1:2*n1) + Ke(5:6,1:2);
    K(2*n3-1:2*n3,2*n2-1:2*n2) = K(2*n3-1:2*n3,2*n2-1:2*n2) + Ke(5:6,3:4);
    K(2*n3-1:2*n3,2*n3-1:2*n3) = K(2*n3-1:2*n3,2*n3-1:2*n3) + Ke(5:6,5:6);
end

% 应用边界条件
for i = 1:size(u_fixed,1)
    K(2*u_fixed(i,1)-1,:) = 0;
    K(:,2*u_fixed(i,1)-1) = 0;
    K(2*u_fixed(i,1)-1,2*u_fixed(i,1)-1) = 1;
    F(2*u_fixed(i,1)-1) = u_fixed(i,2);
end

% 求解位移
u = K\F;

% 计算单元积分点处的应变和应力
for e = 1:size(elements,1)
    % 获取单元节点
    n1 = elements(e,1);
    n2 = elements(e,2);
    n3 = elements(e,3);
    x1 = nodes(n1,1);
    y1 = nodes(n1,2);
    x2 = nodes(n2,1);
    y2 = nodes(n2,2);
    x3 = nodes(n3,1);
    y3 = nodes(n3,2);
    % 计算单元刚度矩阵
    B = [y2-y3, y3-y1, y1-y2; x3-x2, x1-x3, x2-x1];
    A = 0.5*det([1,x1,y1; 1,x2,y2; 1,x3,y3]);
    D = [2*mu+lambda, lambda, 0; lambda, 2*mu+lambda, 0; 0, 0, mu];
    Ke = B'*D*B*A;
    % 计算单元位移向量
    ue = [u(2*n1-1); u(2*n1); u(2*n2-1); u(2*n2); u(2*n3-1); u(2*n3)];
    % 计算单元积分点处的应变和应力
    xi = 1/3;
    eta = 1/3;
    x = xi*x1 + eta*x2 + (1-xi-eta)*x3;
    y = xi*y1 + eta*y2 + (1-xi-eta)*y3;
    J = [x2-x1, x3-x1; y2-y1, y3-y1];
    B = J\B;
    epsilon = B*ue;
    sigma = D*epsilon;
    % 绘制单元并显示应力
    patch([x1,x2,x3],[y1,y2,y3],'w');
    hold on;
    if sigma(1) > 0
        patch([x2,x3],[y2,y3],'r','FaceAlpha',abs(sigma(1))/1e8);
    else
        patch([x3,x1],[y3,y1],'r','FaceAlpha',abs(sigma(1))/1e8);
    end
    if sigma(2) > 0
        patch([x3,x1],[y3,y1],'g','FaceAlpha',abs(sigma(2))/1e8);
    else
        patch([x1,x2],[y1,y2],'g','FaceAlpha',abs(sigma(2))/1e8);
    end
    if sigma(3) > 0
        patch([x1,x2],[y1,y2],'b','FaceAlpha',abs(sigma(3))/1e8);
    else
        patch([x2,x3],[y2,y3],'b','FaceAlpha',abs(sigma(3))/1e8);
    end
end
axis equal;
colorbar;