matlab code for beam finite element analysis in matlab

main.m
% Define the properties of the beam
E = 210e9; % Young's modulus in Pascal
L = 5; % Length of the beam in meters
b = 0.1; % Width of the beam in meters
h = 0.2; % Height of the beam in meters
A = b*h; % Cross-sectional area of the beam
I = (b*h^3)/12; % Moment of inertia of the beam

% Discretize the beam by dividing it into elements
num_elements = 5;
num_nodes = num_elements + 1;
nodes = linspace(0, L, num_nodes); % Node locations

% Assemble the stiffness matrix
K = zeros(2*num_nodes);
for e = 1:num_elements
    n1 = e;
    n2 = e+1;
    Le = nodes(n2) - nodes(n1);
    ke = (E*A/Le) * [1, -1; -1, 1] + (E*I/Le^3) * [12, 6*Le, -12, 6*Le; 6*Le, 4*Le^2, -6*Le, 2*Le^2; -12, -6*Le, 12, -6*Le; 6*Le, 2*Le^2, -6*Le, 4*Le^2];
    
    % Assemble the element stiffness matrix
    idx = [2*n1-1, 2*n1, 2*n2-1, 2*n2];
    K(idx, idx) = K(idx, idx) + ke;
end

% Apply boundary conditions
fixed_dofs = [1,2]; % Assuming both ends are fixed
free_dofs = setdiff(1:2*num_nodes, fixed_dofs);

% Solve for displacements
f_ext = zeros(2*num_nodes, 1); % External force vector (assuming no external force)
u = zeros(2*num_nodes, 1);
u(free_dofs) = K(free_dofs, free_dofs) \ (f_ext(free_dofs) - K(free_dofs, fixed_dofs) * u(fixed_dofs));

% Calculate the reaction forces
f_react = K(fixed_dofs, :) * u;

% Extract element forces and plot the deformed shape of the beam
figure;
hold on;
for e = 1:num_elements
    n1 = e;
    n2 = e+1;
    Le = nodes(n2) - nodes(n1);
    idx = [2*n1-1, 2*n1, 2*n2-1, 2*n2];
    ue = u(idx);
    fe = (E*A/Le) * ([1, -1] * ue(1:2)) + (E*I/Le^2) * ([-6, -3*Le, 6, -3*Le] * ue(3:4));
    
    % Plot the deformed shape of the beam element
    plot([nodes(n1), nodes(n2)], [ue(1), ue(3)], '-o');
end
hold off;
1728 chars
54 lines