-
Notifications
You must be signed in to change notification settings - Fork 1
/
formStiffness_simplysupportedbeam.m
58 lines (49 loc) · 2.26 KB
/
formStiffness_simplysupportedbeam.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
function [stiffness, force, displacements, reactions] = formStiffness_simplysupportedbeam(m, P, L, E, I)
% Parameters
nodeCoordinates = linspace(0, L, m+1)'; % Beam length in meters
numberNodes = size(nodeCoordinates, 1);
EI = E * I; % Modulus of elasticity (Pa) * Moment of inertia (m^4)
GDof = 2 * numberNodes; % Total number of degrees of freedom
% Initialize matrices and vectors
force = zeros(GDof, 1);
stiffness = zeros(GDof);
displacements = zeros(GDof, 1);
% Element connectivity
elementNodes = zeros(m, 2);
for i = 1:m
elementNodes(i, 1) = i;
elementNodes(i, 2) = i + 1;
end
% Assembly of stiffness matrix and force vector
for e = 1:m
indice = elementNodes(e, :);
elementDof = [2*(indice(1)-1)+1, 2*indice(1), 2*(indice(2)-1)+1, 2*indice(2)];
LElem = nodeCoordinates(indice(2)) - nodeCoordinates(indice(1));
% Uniformly Distributed Load (UDL) assumption
f1 = P * LElem * [1/2, LElem/6, 1/2, -LElem/6]';
force(elementDof) = force(elementDof) + f1;
% Stiffness matrix for each element
k1 = EI/(LElem)^3 * [12, 6*LElem, -12, 6*LElem; 6*LElem, 4*LElem^2, -6*LElem, 2*LElem^2; -12, -6*LElem, 12, -6*LElem; 6*LElem, 2*LElem^2, -6*LElem, 4*LElem^2];
stiffness(elementDof, elementDof) = stiffness(elementDof, elementDof) + k1;
end
% Apply boundary conditions for simply supported beam
fixedNodeU = [1, 2*m+1]'; % Vertical displacement constraints at both ends
fixedNodeV = []; % No rotational fixity constraints
prescribedDof = [fixedNodeU; fixedNodeV];
activeDof = setdiff([1:GDof]', prescribedDof);
% Solve for displacements
U = zeros(length(activeDof), 1);
if ~isempty(stiffness(activeDof, activeDof))
U = stiffness(activeDof, activeDof) \ force(activeDof);
end
displacements(activeDof) = U;
% Calculate reactions
F = stiffness * displacements;
reactions = F(prescribedDof);
% Display results
disp('Displacements');
plot(nodeCoordinates, displacements(1:2:GDof), '-o');
xlabel('Position along Beam');
ylabel('Vertical Displacement');
title('Displacement of Simply Supported Beam under Load');
end