FEM/MatlabScript at main · vitspa/FEM · GitHub

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
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
L=0.5;
r=5e-2;
I=pi*r^4/4;
A=pi*r^2;
E=100e6;
beta=70;
lambda=55;
delta=30;
gamma=0;

angles=[beta,lambda,delta,gamma];

%bulding the rod without deformation
x = zeros(5,1);
z = zeros(5,1);

%expressing each node from the previous one
for i = 1:4
    x(i+1) = x(i) + L*cosd(angles(i));
    z(i+1) = z(i) + L*sind(angles(i));
end

plot(x,z);
title("Rod without deformation")


%local stifness matrix
ke=[E*A/L 0 0 -E*A/L 0 0
    0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 6*E*I/L^2
    0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L
    -E*A/L 0 0 E*A/L 0 0
    0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -6*E*I/L^2
    0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 4*E*I/L];

%transformation matrix

for j=1:4

 cosa=(x(j+1)-x(j))/L;
 sina=(z(j+1)-z(j))/L;

 T = [cosa sina  0 0 0 0;
         -sina cosa 0 0 0 0;
         0 0 1 0 0 0;
         0 0 0 cosa sina 0;
         0 0 0 -sina cosa 0;
         0 0 0 0 0 1]

     %local stifness matrix expressed in the global system
     K_local(:,:,j)= T.'*ke*T;

end

%assembling global stifness matrix 5*df=15->15x15
K_global=zeros(15);

K_global(1:6,1:6) = K_local(:,:,1);
K_global(4:9,4:9) = K_global(4:9,4:9) + K_local(:,:,2);
K_global(7:12,7:12) = K_global(7:12,7:12) + K_local(:,:,3);
K_global(10:15,10:15) = K_global(10:15,10:15) + K_local(:,:,4);
K_global;

%creating the force vector
F=zeros(15,1);
F(14,1)=-20;

K_global=K_global(4:15,4:15); %Remove first 3 rows and 3 columns because of boundary condition
F=F(4:15,1); %F will be a 12x1
%solving the linear system to obtain the global deformation vector
D=K_global\F;

%expressing the deflection for each node
dx = zeros(5,1);
dz = zeros(5,1);

for i = 1:4
    dx(i+1) = D(3*(i-1) + 1) + dx(i)
    dz(i+1) = D(3*(i-1) + 2) + dz(i)
end
% Plotting the displacement in x-direction

plot(x,z)
hold on
plot(x+dx,z);
title('Deflection in x-direction')
ylabel(' Z ')
xlabel("X")
legend({'Undeformed', 'Deformed'}, 'Location', 'southeast')
figure
hold off

% Plotting the displacement in z-direction
plot(x,z)
hold on
plot(x,z+dz)
title('Deflection in z-direction')
ylabel(' Z ')
xlabel("X")
legend({'Undeformed', 'Deformed'}, 'Location', 'southeast')
figure
hold off

% Plotting total linear displacement
plot(x,z)
hold on
plot(x+dx, z+dz)
title('Deflection of the rod')
ylabel(' Z ')
xlabel("X")
legend({'Undeformed', 'Deformed'}, 'Location', 'southeast')

%table of the deflections at the nodes

x_def = [dx];
z_def = [dz];
Node = [1;2;3;4;5];
node_deflection = table(Node,x_def,z_def)