Civil engineering coding problem ( I dont know where the problem in my code is and how to solve it, thanks in advance )
21 views (last 30 days)
Show older comments
clear, clc, close all
% Given parameters
a = 4; % width (m)
b = 3; % height (m)
E = 210E9; % E-modulus (Pa)
I = 0.5E-3; % Area moment of inertia (m^4)
R = 1000E3; % Point load (N)
% Element connectivity
edof = [1 2 3 4; 3 4 5 6 ]; % Corrected [EL1; EL2]
% Node coordinates
coord = [0 0; a 0; a b]; % Corrected coordinates
% Boundary conditions
fdof = [1 2 3 4]; % Free dofs
fixed_dofs = [5 6]; % Fixed dofs
% Initialize global stiffness matrix
K = zeros(6,6);
% Global load vector
F = zeros(6,1);
F(4) = -R;
% Assemble global stiffness matrix
for i = 1:size(edof,1)
n1 = edof(i,1)
n2 = edof(i,2)
x1 = coord(n1,1)
y1 = coord(n1,2)
x2 = coord(n2,1)
y2 = coord(n2,2)
L = sqrt((x2-x1)^2 + (y2-y1)^2);
c = (x2-x1)/L;
s = (y2-y1)/L;
ke = beam2d_stiff(E,I,L,c,s);
dofs = [2*n1-1 2*n1 2*n2-1 2*n2];
K(dofs,dofs) = K(dofs,dofs) + ke;
end
% Apply boundary conditions
Kf = K(fdof,fdof);
Ff = F(fdof) - K(fdof,fixed_dofs)*zeros(length(fixed_dofs),1);
% Solve for displacements
Df = Kf\Ff;
% Compute reaction forces
Rf = K(fdof,fixed_dofs)*zeros(length(fixed_dofs),1) - F(fdof);
% Compute axial forces
N = zeros(2,1);
for i = 1:size(edof,1)
n1 = edof(i,1);
n2 = edof(i,2);
x1 = coord(n1,1);
y1 = coord(n1,2);
x2 = coord(n2,1);
y2 = coord(n2,2);
L = sqrt((x2-x1)^2 + (y2-y1)^2);
c = (x2-x1)/L;
s = (y2-y1)/L;
ue = [Df(2*n1-1); Df(2*n1); Df(2*n2-1); Df(2*n2)];
[N(i),~,~] = beam2d_forces(E,I,L,c,s,ue);
end
% Display results
disp(' Node dx (mm) dy (mm)')
disp(' ---------------------------------')
disp([(1:3)' Df(1:2:end)*1E3 Df(2:2:end)*1E3])
disp(' ');
disp(' Node Rx (kN) Ry (kN)')
disp(' ---------------------------------')
disp([(1:3)' Rf(1:2:end)/1E3 Rf(2:2:end)/1E3])
disp(' ');
disp(' EL N (kN)');
disp(' -------------------');
disp([(1:2)' N/1E3]);
% Plot geometry and displacements
fac = 500; % scale factor, deformed shape
figure(1), clf
x = [0 a a+b]; y = [0 b 0]; % coordinates, undeformed geometry
dx = Df(1:2:end); dy = Df(2:2:end); % nodal displacements
plot(x,y,'ko-','LineWidth',1.5), hold all, axis equal
plot(x+fac*dx', y+fac*dy','Color',[0.4 0.4 0.4])
title('\rm 2D bar model')
legend('undeformed','deformed')
% Beam stiffness matrix function
function ke = beam2d_stiff(E,I,L,c,s)
ke = E*I/L^3 * [ c^2*L^2 c*s*L^2 -c^2*L^2 -c*s*L^2;
c*s*L^2 s^2*L^2 -c*s*L^2 -s^2*L^2;
-c^2*L^2 -c*s*L^2 c^2*L^2 c*s*L^2;
-c*s*L^2 -s^2*L^2 c*s*L^2 s^2*L^2 ];
end
% Beam forces calculation function
function [N,V,M] = beam2d_forces(E,I,L,c,s,ue)
N = E*I/L^2 * [-c -s c s] * ue;
V = E*I/L^2 * [-s c -s c] * ue;
M = E*I/L * [ s -c -s c] * ue;
end
1 Comment
Ayush Modi
on 16 Apr 2024 at 11:00
Hi Mahmoud,
You haven't mentioned "What" is the problem. Please elaborate on your problem.
Also refer the below link to see how to ask a question (on Answers) and get a fast answer -
See Also
Categories
Find more on Configure Serial Link Projects in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!