有限元分析悬臂梁matlab程序

以下是一个简单的有限元分析悬臂梁的Matlab程序：

% 定义参数
L = 1; % 梁的长度
h = 0.1; % 梁的高度
b = 0.05; % 梁的宽度
E = 2e11; % 弹性模量
nu = 0.3; % 泊松比
rho = 7800; % 密度
g = 9.81; % 重力加速度

% 定义节点和单元数
n_nodes = 11;
n_elems = n_nodes - 1;

% 定义节点坐标
x = linspace(0, L, n_nodes);
y = zeros(1, n_nodes);

% 定义单元节点
elems = [1 2; 2 3; 3 4; 4 5; 5 6; 6 7; 7 8; 8 9; 9 10; 10 11];

% 定义材料矩阵
D = E/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2];

% 定义初始应变矩阵
epsilon_0 = zeros(3, n_elems);

% 定义初始速度和加速度矩阵
v_0 = zeros(3, n_nodes);
a_0 = zeros(3, n_nodes);

% 定义质量矩阵和刚度矩阵
M = zeros(3*n_nodes, 3*n_nodes);
K = zeros(3*n_nodes, 3*n_nodes);

% 计算质量矩阵和刚度矩阵
for i = 1:n_elems
    % 获取单元节点
    n1 = elems(i, 1);
    n2 = elems(i, 2);
    
    % 计算单元长度和角度
    L_e = norm([x(n2)-x(n1) y(n2)-y(n1)]);
    theta_e = atan2(y(n2)-y(n1), x(n2)-x(n1));
    
    % 计算单元转换矩阵
    T_e = [cos(theta_e) sin(theta_e) 0 0 0 0;
           -sin(theta_e) cos(theta_e) 0 0 0 0;
           0 0 1 0 0 0;
           0 0 0 cos(theta_e) sin(theta_e) 0;
           0 0 0 -sin(theta_e) cos(theta_e) 0;
           0 0 0 0 0 1];
    
    % 计算单元质量矩阵和刚度矩阵
    Me = rho*h*b*L_e/420*[156 22*L_e 54 -13*L_e;
                          22*L_e 4*L_e^2 13*L_e -3*L_e^2;
                          54 13*L_e 156 -22*L_e;
                          -13*L_e -3*L_e^2 -22*L_e 4*L_e^2];
                      
    Ke = h*b/(L_e*(1-nu^2))*[1/3*E 1/6*E*L_e -1/3*E 1/6*E*L_e;
                             1/6*E*L_e 1/3*E*L_e^2 -1/6*E*L_e -1/3*E*L_e^2;
                             -1/3*E -1/6*E*L_e 1/3*E -1/6*E*L_e;
                             1/6*E*L_e -1/3*E*L_e^2 -1/6*E*L_e 1/3*E*L_e^2];
    
    % 转换到全局坐标系
    Me = T_e' * Me * T_e;
    Ke = T_e' * Ke * T_e;
    
    % 将单元质量矩阵和刚度矩阵加入全局质量矩阵和刚度矩阵
    M(3*n1-2:3*n1, 3*n1-2:3*n1) = M(3*n1-2:3*n1, 3*n1-2:3*n1) + Me(1:3, 1:3);
    M(3*n1-2:3*n1, 3*n2-2:3*n2) = M(3*n1-2:3*n1, 3*n2-2:3*n2) + Me(1:3, 4:6);
    M(3*n2-2:3*n2, 3*n1-2:3*n1) = M(3*n2-2:3*n2, 3*n1-2:3*n1) + Me(4:6, 1:3);
    M(3*n2-2:3*n2, 3*n2-2:3*n2) = M(3*n2-2:3*n2, 3*n2-2:3*n2) + Me(4:6, 4:6);
    
    K(3*n1-2:3*n1, 3*n1-2:3*n1) = K(3*n1-2:3*n1, 3*n1-2:3*n1) + Ke(1:3, 1:3);
    K(3*n1-2:3*n1, 3*n2-2:3*n2) = K(3*n1-2:3*n1, 3*n2-2:3*n2) + Ke(1:3, 4:6);
    K(3*n2-2:3*n2, 3*n1-2:3*n1) = K(3*n2-2:3*n2, 3*n1-2:3*n1) + Ke(4:6, 1:3);
    K(3*n2-2:3*n2, 3*n2-2:3*n2) = K(3*n2-2:3*n2, 3*n2-2:3*n2) + Ke(4:6, 4:6);
end

% 定义边界条件
bc_nodes = [1];
bc_dofs = [1 2 3];

% 将边界条件应用到质量矩阵和刚度矩阵上
for i = 1:length(bc_nodes)
    node = bc_nodes(i);
    dofs = bc_dofs;
    
    M(3*node-dofs, :) = 0;
    M(:, 3*node-dofs) = 0;
    M(3*node-dofs, 3*node-dofs) = eye(length(dofs));
    
    K(3*node-dofs, :) = 0;
    K(:, 3*node-dofs) = 0;
    K(3*node-dofs, 3*node-dofs) = eye(length(dofs));
end

% 定义外力矩阵
F = zeros(3*n_nodes, 1);
F(3*n_nodes-1) = -1000;

% 求解加速度矩阵
a = inv(M) * (F - K * v_0(:));

% 求解位移矩阵
v = v_0(:) + a_0(:) + a(:);

% 输出结果
disp(v);

% 绘制梁的形状
figure;
hold on;
for i = 1:n_elems
    n1 = elems(i, 1);
    n2 = elems(i, 2);
    plot([x(n1) x(n2)], [y(n1) y(n2)], 'k');
end
axis equal;

该程序定义了一个悬臂梁的几何和材料参数，然后使用有限元方法计算了梁的位移响应。程序输出了梁的位移矩阵，并绘制了梁的形状。