(no subject)
So, I think my Finite Element project might be finished. I still need to run it through Ansys and compare, but so far so good. The problem involves a 4x6 inch plate with a 1.5 inch diameter hole in the middle. The plate has a pressure of 2000psi outward applied at the left and right edges. The code was 196 lines in Matlab (including comments). Here's the solution:
1 % Stuart McDonald
2 % ME 461
3 % Final Project
4
5 % Clear everything
6 clear all;
7 clc;
8 clf;
9
10 % Physical dimensions
11 DELTA = .1;
12 DIAMETER = 1.5;
13 RADIUS = DIAMETER/2;
14 x = 0:DELTA:6;
15 y = 0:DELTA:4;
16
17 % Material properties
18 E = 30E6;
19 v = 0.3;
20 c = (E/(1-v^2))*[1 v 0; v 1 0; 0 0 (1-v)/2];
21
22 % Forces on the sides
23 f = 2000*4/length(y);
24
25 % The center of the circle
26 cx = (length(x)-1)/2+1;
27 cy = (length(y)-1)/2+1;
28
29 % Weighting values
30 w1 = 1;
31 w2 = 1;
32
33 disp('Setting up the mesh...');
34
35 % Define the array of nodes
36 % 0 if the node is in the circle
37 nodes = ones(length(y), length(x));
38
39 % Find the elements in the circle
40 for i1 = 1:length(x)
41 for i2 = 1:length(y)
42 if sqrt((i2-cy)^2 + (i1-cx)^2)*DELTA <= RADIUS
43 nodes(i2, i1) = 0;
44 end
45 end
46 end
47 clear cx cy i1 i2;
48
49 disp('Mapping nodes to elements...');
50
51 % Nodes at the corners of each element
52 El1 = zeros((length(y)-1),(length(x)-1));
53 El2 = zeros((length(y)-1),(length(x)-1));
54 El3 = zeros((length(y)-1),(length(x)-1));
55 El4 = zeros((length(y)-1),(length(x)-1));
56
57 % Map nodes to elements
58 for i1 = 1:length(x)-1
59 for i2 = 1:length(y)-1
60 tl = (i1-1)*length(y)+i2; bl = tl+1;
61 tr = (i1)*length(y)+i2; br = tr+1;
62 el = (i1-1)*(length(y)-1)+i2;
63 if (nodes(tl)==1) && (nodes(bl)==1) && (nodes(tr)==1) && (nodes(br)==1)
64 El1(el)=bl; El2(el)=br; El3(el)=tr; El4(el)=tl;
65 end
66 end
67 end
68 clear i1 i2 tl bl tr br el;
69
70 disp('Building the global stiffness matrix...');
71
72 % Build the global stiffness matrix
73 K = zeros(numel(nodes)*2, numel(nodes)*2);
74 for loop = 1:numel(El1)
75 if El1(loop) ~= 0
76 X = [0 DELTA DELTA 0; 0 0 DELTA DELTA];
77 ke = zeros(8,8);
78
79 % Build the elemental stiffness matrix
80 for i = -1:2:1
81 for e = -1:2:1
82 r1 = i*(1/sqrt(3));
83 r2 = e*(1/sqrt(3));
84 dn = (1/4)*[-(1-r2) (1-r2) (1+r2) -(1+r2); -(1-r1) -(1+r1) (1+r1) (1-r1)];
85 dns = (1/4)*[-(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2) 0; -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1) 0;...
86 0 -(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2); 0 -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1);];
87 j = X*dn';
88 j1t = inv(j)';
89 a = [j1t(1,1) j1t(1,2) 0 0; 0 0 j1t(2,1) j1t(2,2); j1t(2,1) j1t(2,2) j1t(1,1) j1t(1,2)];
90 b = a*dns;
91 ke = ke+(b'*c*b)*det(j)*w1*w2;
92 end
93 end
94 clear i e X r1 r2 dn dns j j1t a b k;
95
96 % Mapping of local nodes to global nodes
97 map=[2*El1(loop)-1; 2*El1(loop); 2*El2(loop)-1; 2*El2(loop);...
98 2*El3(loop)-1; 2*El3(loop); 2*El4(loop)-1; 2*El4(loop);];
99 for i = 1:8
100 for e = 1:8
101 K(map(i), map(e)) = K(map(i), map(e)) + ke(i, e);
102 end
103 end
104 end
105 end
106 clear i e ke map w1 w2;
107
108 disp('Building the force vector...');
109
110 % Setup the force vector
111 F = zeros(numel(nodes)*2, 1);
112 for i = 1:length(y)
113 left = i;
114 right = (length(x)-1)*length(y)+i;
115 F(2*left-1) = -f;
116 F(2*right-1) = f;
117 end
118 clear left right i;
119
120 disp('Reducing the problem...');
121
122 % Reduce the stiffness matrix and force vector
123 removed = [];
124 for loop = numel(nodes):-1:1
125 if nodes(loop) == 0
126 removed = [2*loop-1; 2*loop; removed(:)];
127 K(2*loop, :) = [];
128 K(2*loop-1, :) = [];
129 K(:, 2*loop) = [];
130 K(:, 2*loop-1) = [];
131 F(2*loop) = [];
132 F(2*loop-1) = [];
133 end
134 end
135 clear loop;
136
137 disp('Solving for the displacement vector...');
138
139 % Solve for the reduced displacement vector
140 Q = K\F;
141
142 % Rebuild the full displacement vector
143 for i=1:length(removed)
144 Q = [Q(1:removed(i)-1); 0; Q(removed(i):length(Q))];
145 end
146 clear i removed;
147
148 disp('Plotting the solution...');
149
150 % Find the initial locations and displaced locations
151 X = zeros(numel(nodes),1);
152 Y = zeros(numel(nodes),1);
153 Xd = zeros(numel(nodes),1);
154 Yd = zeros(numel(nodes),1);
155 for i = 1:numel(nodes)
156 xval = floor((i-1)/length(y))+1;
157 yval = i-length(y)*(xval-1);
158 X(i) = (xval-1)*DELTA;
159 Y(i) = (yval-1)*DELTA;
160 Xd(i) = X(i) + Q(2*i-1);
161 Yd(i) = Y(i) + Q(2*i);
162 end
163 clear xval yval i;
164
165 % Zero out the nodes that don't exist
166 X = X.*nodes(:); Xd = Xd.*nodes(:);
167 Y = Y.*nodes(:); Yd = Yd.*nodes(:);
168
169 % Solve for stresses and plot the displaced shape along with the stresses
170 hold on;
171 for loop = 1:numel(El1)
172 if El1(loop) ~= 0
173 xe = [0 DELTA DELTA 0; 0 0 DELTA DELTA];
174 q = [Q(2*El1(loop)-1); Q(2*El1(loop)); Q(2*El2(loop)-1); Q(2*El2(loop));...
175 Q(2*El3(loop)-1); Q(2*El3(loop)); Q(2*El4(loop)-1); Q(2*El4(loop));];
176 r1 = 0;
177 r2 = 0;
178 dn = (1/4)*[-(1-r2) (1-r2) (1+r2) -(1+r2); -(1-r1) -(1+r1) (1+r1) (1-r1)];
179 dns = (1/4)*[-(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2) 0; -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1) 0;...
180 0 -(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2); 0 -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1);];
181 j = xe*dn';
182 j1t = inv(j)';
183 a = [j1t(1,1) j1t(1,2) 0 0; 0 0 j1t(2,1) j1t(2,2); j1t(2,1) j1t(2,2) j1t(1,1) j1t(1,2)];
184 b = a*dns;
185 sigma = c*b*q;
186 stress = zeros(1,2);
187 stress(1) = ((sigma(1)+sigma(2))/2) + sqrt(((sigma(1)-sigma(2))/2).^2+sigma(3).^2);
188 stress(2) = ((sigma(1)+sigma(2))/2) - sqrt(((sigma(1)-sigma(2))/2).^2+sigma(3).^2);
189 fill([Xd(El1(loop)) Xd(El2(loop)) Xd(El3(loop)) Xd(El4(loop))],...
190 [Yd(El1(loop)) Yd(El2(loop)) Yd(El3(loop)) Yd(El4(loop))], ones(1,4)*max(stress));
191 end
192 end
193 axis auto; axis equal; colorbar;
194 clear loop xe r1 r2 dn dns j j1t a b sigma stress;
195
196 disp('Done');
1 % Stuart McDonald
2 % ME 461
3 % Final Project
4
5 % Clear everything
6 clear all;
7 clc;
8 clf;
9
10 % Physical dimensions
11 DELTA = .1;
12 DIAMETER = 1.5;
13 RADIUS = DIAMETER/2;
14 x = 0:DELTA:6;
15 y = 0:DELTA:4;
16
17 % Material properties
18 E = 30E6;
19 v = 0.3;
20 c = (E/(1-v^2))*[1 v 0; v 1 0; 0 0 (1-v)/2];
21
22 % Forces on the sides
23 f = 2000*4/length(y);
24
25 % The center of the circle
26 cx = (length(x)-1)/2+1;
27 cy = (length(y)-1)/2+1;
28
29 % Weighting values
30 w1 = 1;
31 w2 = 1;
32
33 disp('Setting up the mesh...');
34
35 % Define the array of nodes
36 % 0 if the node is in the circle
37 nodes = ones(length(y), length(x));
38
39 % Find the elements in the circle
40 for i1 = 1:length(x)
41 for i2 = 1:length(y)
42 if sqrt((i2-cy)^2 + (i1-cx)^2)*DELTA <= RADIUS
43 nodes(i2, i1) = 0;
44 end
45 end
46 end
47 clear cx cy i1 i2;
48
49 disp('Mapping nodes to elements...');
50
51 % Nodes at the corners of each element
52 El1 = zeros((length(y)-1),(length(x)-1));
53 El2 = zeros((length(y)-1),(length(x)-1));
54 El3 = zeros((length(y)-1),(length(x)-1));
55 El4 = zeros((length(y)-1),(length(x)-1));
56
57 % Map nodes to elements
58 for i1 = 1:length(x)-1
59 for i2 = 1:length(y)-1
60 tl = (i1-1)*length(y)+i2; bl = tl+1;
61 tr = (i1)*length(y)+i2; br = tr+1;
62 el = (i1-1)*(length(y)-1)+i2;
63 if (nodes(tl)==1) && (nodes(bl)==1) && (nodes(tr)==1) && (nodes(br)==1)
64 El1(el)=bl; El2(el)=br; El3(el)=tr; El4(el)=tl;
65 end
66 end
67 end
68 clear i1 i2 tl bl tr br el;
69
70 disp('Building the global stiffness matrix...');
71
72 % Build the global stiffness matrix
73 K = zeros(numel(nodes)*2, numel(nodes)*2);
74 for loop = 1:numel(El1)
75 if El1(loop) ~= 0
76 X = [0 DELTA DELTA 0; 0 0 DELTA DELTA];
77 ke = zeros(8,8);
78
79 % Build the elemental stiffness matrix
80 for i = -1:2:1
81 for e = -1:2:1
82 r1 = i*(1/sqrt(3));
83 r2 = e*(1/sqrt(3));
84 dn = (1/4)*[-(1-r2) (1-r2) (1+r2) -(1+r2); -(1-r1) -(1+r1) (1+r1) (1-r1)];
85 dns = (1/4)*[-(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2) 0; -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1) 0;...
86 0 -(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2); 0 -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1);];
87 j = X*dn';
88 j1t = inv(j)';
89 a = [j1t(1,1) j1t(1,2) 0 0; 0 0 j1t(2,1) j1t(2,2); j1t(2,1) j1t(2,2) j1t(1,1) j1t(1,2)];
90 b = a*dns;
91 ke = ke+(b'*c*b)*det(j)*w1*w2;
92 end
93 end
94 clear i e X r1 r2 dn dns j j1t a b k;
95
96 % Mapping of local nodes to global nodes
97 map=[2*El1(loop)-1; 2*El1(loop); 2*El2(loop)-1; 2*El2(loop);...
98 2*El3(loop)-1; 2*El3(loop); 2*El4(loop)-1; 2*El4(loop);];
99 for i = 1:8
100 for e = 1:8
101 K(map(i), map(e)) = K(map(i), map(e)) + ke(i, e);
102 end
103 end
104 end
105 end
106 clear i e ke map w1 w2;
107
108 disp('Building the force vector...');
109
110 % Setup the force vector
111 F = zeros(numel(nodes)*2, 1);
112 for i = 1:length(y)
113 left = i;
114 right = (length(x)-1)*length(y)+i;
115 F(2*left-1) = -f;
116 F(2*right-1) = f;
117 end
118 clear left right i;
119
120 disp('Reducing the problem...');
121
122 % Reduce the stiffness matrix and force vector
123 removed = [];
124 for loop = numel(nodes):-1:1
125 if nodes(loop) == 0
126 removed = [2*loop-1; 2*loop; removed(:)];
127 K(2*loop, :) = [];
128 K(2*loop-1, :) = [];
129 K(:, 2*loop) = [];
130 K(:, 2*loop-1) = [];
131 F(2*loop) = [];
132 F(2*loop-1) = [];
133 end
134 end
135 clear loop;
136
137 disp('Solving for the displacement vector...');
138
139 % Solve for the reduced displacement vector
140 Q = K\F;
141
142 % Rebuild the full displacement vector
143 for i=1:length(removed)
144 Q = [Q(1:removed(i)-1); 0; Q(removed(i):length(Q))];
145 end
146 clear i removed;
147
148 disp('Plotting the solution...');
149
150 % Find the initial locations and displaced locations
151 X = zeros(numel(nodes),1);
152 Y = zeros(numel(nodes),1);
153 Xd = zeros(numel(nodes),1);
154 Yd = zeros(numel(nodes),1);
155 for i = 1:numel(nodes)
156 xval = floor((i-1)/length(y))+1;
157 yval = i-length(y)*(xval-1);
158 X(i) = (xval-1)*DELTA;
159 Y(i) = (yval-1)*DELTA;
160 Xd(i) = X(i) + Q(2*i-1);
161 Yd(i) = Y(i) + Q(2*i);
162 end
163 clear xval yval i;
164
165 % Zero out the nodes that don't exist
166 X = X.*nodes(:); Xd = Xd.*nodes(:);
167 Y = Y.*nodes(:); Yd = Yd.*nodes(:);
168
169 % Solve for stresses and plot the displaced shape along with the stresses
170 hold on;
171 for loop = 1:numel(El1)
172 if El1(loop) ~= 0
173 xe = [0 DELTA DELTA 0; 0 0 DELTA DELTA];
174 q = [Q(2*El1(loop)-1); Q(2*El1(loop)); Q(2*El2(loop)-1); Q(2*El2(loop));...
175 Q(2*El3(loop)-1); Q(2*El3(loop)); Q(2*El4(loop)-1); Q(2*El4(loop));];
176 r1 = 0;
177 r2 = 0;
178 dn = (1/4)*[-(1-r2) (1-r2) (1+r2) -(1+r2); -(1-r1) -(1+r1) (1+r1) (1-r1)];
179 dns = (1/4)*[-(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2) 0; -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1) 0;...
180 0 -(1-r2) 0 (1-r2) 0 (1+r2) 0 -(1+r2); 0 -(1-r1) 0 -(1+r1) 0 (1+r1) 0 (1-r1);];
181 j = xe*dn';
182 j1t = inv(j)';
183 a = [j1t(1,1) j1t(1,2) 0 0; 0 0 j1t(2,1) j1t(2,2); j1t(2,1) j1t(2,2) j1t(1,1) j1t(1,2)];
184 b = a*dns;
185 sigma = c*b*q;
186 stress = zeros(1,2);
187 stress(1) = ((sigma(1)+sigma(2))/2) + sqrt(((sigma(1)-sigma(2))/2).^2+sigma(3).^2);
188 stress(2) = ((sigma(1)+sigma(2))/2) - sqrt(((sigma(1)-sigma(2))/2).^2+sigma(3).^2);
189 fill([Xd(El1(loop)) Xd(El2(loop)) Xd(El3(loop)) Xd(El4(loop))],...
190 [Yd(El1(loop)) Yd(El2(loop)) Yd(El3(loop)) Yd(El4(loop))], ones(1,4)*max(stress));
191 end
192 end
193 axis auto; axis equal; colorbar;
194 clear loop xe r1 r2 dn dns j j1t a b sigma stress;
195
196 disp('Done');