辅导Matlab Assignmen、辅导A Matlab program程序 辅导A Matlab program for a 3D structure
- 首页 >> Matlab编程Assignment 1: A Matlab program for a 3D structure
Please refer the attached 2D frame program (Appendix A) with line-by-line explanations to develop a Matlab program
for a 3D Frame, whose nodal and elemental information are shown in Appendix B. The shape of three types of cross
section is given in Appendix C. The loading and boundary conditions are shown in Appendix D. Appendix E shows a
schematic of frame element in 3D, the local stiffness matrix and material properties for the 3D structure. Noted that the
gravity is not considered. In Appendix F, a sample Transmission tower is plotted and your program can be used
directly for this complex structure. In Appendix G, analysis of 3D frame in Matlab is given via two examples.
The assignment should:
1) Write in terms of https://emedia.rmit.edu.au/learninglab/content/reports-0. A sample report is in Appendix
B (10 point)
2) Introduce the theory of Finite Element Analysis for 3D frame. (10 point)
3) Plot the flow chart of the program. (10 point)
4) The structure is composed of 3 types of members. (15 point)
5) Plot the 3D Frame before deformation in a 3D view. The type of different members should be clearly seen in
the figure. (15 point)
6) List the deformation of all nodes in a table. (20 point)
7) Write your understanding in Finite Element Analysis through this program. (10 point)
8) Please submit the report, Matlab program, input files (node.txt, elem.txt and etc.). (10 point)
Please upload the assignment to the blackboard before 00:00 19 Aug 2018. Late submissions will attract a penalty of
5% per day.
APPENDIX A
Finite element analysis for a 2D Frame in Matlab
Line Code explanation
1 clc; Clear command window.
2 clear; removes all variables from the workspace.
3 close all; closes all the open figure windows.
4 % blank line
5 N = textread('node.txt'); Read x and y coordinate for the node and assign
their values to N
6 E = textread('elem.txt'); Read the near and far end for the beam and
assign
7 % blank line
8 NN = length(N); Get the number of node
9 NE = length(E); Get the number of element
10 % blank line
11 E0 = 2.1e6; Set Young’s modulus
12 A0 = 22.8; Set the area of cross section
13 I0 = 935; Set the second moment of inertia
14 % blank line
15 hold on;
holds the current plot and all axis properties so
that subsequent graphing commands add to the
existing graph.
16 for i=1:NE Begin an iteration in which i begins from 1 to NE
17 plot([N(E(i,1),1),N(E(i,2),1)],[N(E(i,1),2),N(
E(i,2),2)],'k','linewidth',2);
Plot a line starts from the near end to the far
end of an element in black with the linewidth of
2.
18
text((N(E(i,1),1)+N(E(i,2),1))/2,(N(E(i,1),2)+
N(E(i,2),2))/2,num2str(i),'color','b','fontsiz
e',12);
Plot the element number at its middle in blue
and 12 font size.
19 end End the iteration
20 % blank line
21 for i=1:NN Begin an iteration in which i begins from 1 to NN
22 if i<3 If i<3
23 plot(N(i,1),N(i,2),'rs','markersize',12,'mark
erfacecolor','r');
Plot a red square marker at the node i with the
size of 12 and red marker face color.
24 else Else
25 plot(N(i,1),N(i,2),'ko','markersize',8,'mark
erfacecolor','k');
Plot a black circle marker at the node i with the
size of 8 and black marker face color.
26 end End if
27 text(N(i,1),N(i,2),num2str(i),'color','g','font
size',12);
Plot a green node number at the node i with the
size of 12.
28 end End the iteration
29 % blank line
30 F=zeros(3*NN,1); Create the force vector and all members are
zero.
31 FN = [11 12 17 18 23 24]; Set the node number at which a load is applied
32 F(3*FN-2) = -20; The x component of the load is -20
33 F(3*FN-1) = -10; The y component of the load is -10
34 % blank line
35 scale = 20; Set the scaling factor
36 for i=1:length(FN) Begin an iteration in which i begins from 1 to the
number of loads
37 plot( N(FN(i),1),N(FN(i),2),'b*','markersize'
,10,'markerfacecolor','b');
Plot a blue star with the size of 10 and blue
marker face color at the node where a load is
applied.
38
plot([N(FN(i),1),N(FN(i),1)-
100*sin(20/180*pi)],[N(FN(i),2),N(FN(i),2)
-100*cos(20/180*pi)],'-
>r','linewidth',2,'markerfacecolor','r');
Plot red “>” symbols to show the applied load in
20 degree angle with line width of 2 and red
marker face color
39 % blank line
40 ang=-1.5*pi:0.5:-.5*pi;
41 xp=scale*cos(ang);
42 yp=scale*sin(ang);
43 plot(N(FN(i),1)+xp,N(FN(i),2)+yp,'--
r>','linewidth',2,'markerfacecolor','r');
Plot 7 red “>” symbols to show the applied
moment with line width of 2 and red marker
face color
44 end End the iteration
45 % blank line
46 title('A schematic of the power
tower','fontsize',12);
Plot the title of this figure with 12 font size
47 set(gca,'fontsize',12); Set the font size of current axis as 12
48 axis equal tight Set the axis property as “equal” and “tight”
49 print('-dtiff','-r400','frame.tif'); Print the figure in the name of “frame.tif” with
tiff format and 400 resolution
50 % blank line
51 K=zeros(3*NN,3*NN); Determine the global stiffness matrix, its size is
3NN by 3NN and all members are zero.
52 for n=1:NE Begin an iteration in which i begins from 1 to NN
53 L = sqrt((N(E(n,1),1)-N(E(n,2),1))^2 +
(N(E(n,1),2)-N(E(n,2),2))^2);
Calculate the length of element
54 lx = (N(E(n,2),1) - N(E(n,1),1))/L; Calculate cos(theta)
55 ly = (N(E(n,2),2) - N(E(n,1),2))/L; Calculate sin(theta)
56 a=A0*E0/L; Calculate some parameters for the local matrix
57 b=6*E0*I0/L^2; Calculate some parameters for the local matrix
58 c=12*E0*I0/L^3; Calculate some parameters for the local matrix
59 d=E0*I0/L; Calculate some parameters for the local matrix
60 % blank line
61 Ke = [ a 0 0 -a 0 0 Determine the local matrix at local coordinate
system
62 0 c b 0 -c b
63 0 b 4*d 0 -b 2*d
64 -a 0 0 a 0 0
65 0 -c -b 0 c -b
66 0 b 2*d 0 -b 4*d];
67 % blank line Determine the transformation matrix
68 T = [ lx ly 0 0 0 0
69 -ly lx 0 0 0 0
70 0 0 1 0 0 0
71 0 0 0 lx ly 0
72 0 0 0 -ly lx 0
73 0 0 0 0 0 1];
74 % blank line
75 Ke = T'*Ke*T; Determine the local matrix at global coordinate
system
76 i = E(n,1); Get the near end number of an element
77 j = E(n,2); Get the far end number of an element
78 dof = [3*i-2:3*i 3*j-2:3*j]; 6 degrees of the element
79 K(dof,dof) = K(dof,dof) + Ke; Assemble the global matrix from the local matrix
80 end End the iteration
81 % blank line
82 fixeddofs = [1 2 3 4 5 6]; Fixed degrees
83 alldofs = [1:3*NN]; All degrees
84 freedofs = setdiff(alldofs,fixeddofs); Get the active degrees
85 U(freedofs,:) =
K(freedofs,freedofs)\F(freedofs,:);
Calculate the equation Ax=b to get the
deformations for all active degrees.
86 U(fixeddofs,:)= 0; Make the fixed degrees equal zero
87 % blank line
88 figure; Plot a new figure
89 hold on;
holds the current plot and all axis properties so
that subsequent graphing commands add to the
existing graph.
90 ND = N + scale*[U(1:3:end) U(2:3:end)]; Plot the node after the deformation. Noted the
deformation is multiplied by the scaling factor.
91 % blank line
92 for i=1:NE Begin an iteration in which i begins from 1 to NN
93 plot([N(E(i,1),1),N(E(i,2),1)],[N(E(i,1),2),N(
E(i,2),2)],'linewidth',2,'color','k');
Plot a line starts from the near end (before
deformation) to the far end (before
deformation) of an element in black with the
line width of 2.
94 plot([ND(E(i,1),1),ND(E(i,2),1)],[ND(E(i,1),2
),ND(E(i,2),2)],'--b','linewidth',2);
Plot a line starts from the near end (after
deformation)to the far end (after deformation)
of an element in blue with the dashed line width
of 2.
95 end End the iteration
96 % blank line
97 title(['The deformation magnified '
int2str(scale) ' times'],'fontsize',12);
Plot the title of this figure with 12 font size
98 axis equal tight off Set the font size of current axis as 12
99 set(gca,'fontsize',12); Set the axis property as “equal” and “tight”
100 print('-dtiff','-r400','U_frame.tif'); Print the figure in the name of “U_frame.tif”
with tiff format and 400 resolution
APPENDIX B
Node x y z
1 -0.5 -0.3 0
2 0.5 -0.3 0
3 0.5 0.3 0
4 -0.5 0.3 0
5 -0.5 -0.3 2
6 0.5 -0.3 2
7 0.5 0.3 2
8 -0.5 0.3 2
Element
number
Element
Type
Final
Node
Initial
Node
1 1 1 5
2 1 2 6
3 1 3 7
4 1 4 8
5 2 5 6
6 2 6 7
7 2 7 8
8 2 8 5
9 3 1 6
10 3 2 7
11 3 3 8
12 3 4 5
APPENDIX C
The cross section of three elements in the shape of Equal Angles, Welded Columns, and Circular Rings
APPENDIX D
Nodal External Loads Vector Nodal External Loads Vector
Node Degrees of
Freedom External Load
5
25 Px 1000
26 Py 2000
27 Pz 3000
28 Mx 0
29 My 0
30 Mz 0
6
31 Px 1000
32 Py 2000
33 Pz 3000
34 Mx 0
35 My 0
36 Mz 0
7
37 Px 1000
38 Py 2000
39 Pz 3000
40 Mx 0
41 My 0
42 Mz 0
8
43 Px 1000
44 Py 2000
45 Pz 3000
46 Mx 0
47 My 0
48 Mz 0
Fixed Node Fixed Node xed Node
Node Degrees of
Freedom External Load
1
1 u1 0
2 v1 0
3 w1 0
4 θx1 0
5 θy1 0
6 θz1 0
2
7 u2 0
8 v2 0
9 w2 0
10 θx2 0
11 θy2 0
12 θz2 0
3
13 u3 0
14 v3 0
15 w3 0
16 θx3 0
17 θy3 0
18 θz3 0
4
19 u4 0
20 v4 0
21 w4 0
22 θx4 0
23 θy4 0
24 θz4 0
APPENDIX E:
Structural Steelwork (AS 4100-1998) pa
Young’s modulus: E=2e11 pa
Shear Modulus: 8e10 pa
APPENDIX F: Sample Transmission tower
APPENDIX G: Finite Element Analysis of 3D Frame
104 8 Analysis of 3D frames
After transformation to the global axes, the stiffness matrix in global coordinates
is obtained as
K = RT K
R
where the rotation matrix R is defined as
R =
⎡
⎢
⎢
⎣
r000
0r00
00r0
000r
⎤
⎥
⎥
⎦
(8.2)
being
r =
⎡
⎣
CXx CY x CZx
CXy CY y CZy
CXz CY z CZz
⎤
⎦ (8.3)
and
CXx = cosθXx
where angles θXx, θY x, and θZx, are measured from global axes X, Y , and Z,
with respect to the local axis x, respectively. The two-node 3D frame element has
six degrees of freedom per node. Given the nodal displacements, it is possible to
calculate the reactions at the supports by
F = KU (8.4)
where K and U are the structure stiffness matrix and the vector of nodal displacement,
respectively. The element forces are also evaluated by axes transformation as
fe = keRUe (8.5)
8.3 First 3D frame example
The first 3D frame example is illustrated in figure 8.1. We consider E = 210
GPa, G = 84 GPa, A = 2 × 10−2 m2, Iy = 10 × 10−5 m4, Iz = 20 × 10−5 m4,
J = 5 × 10−5 m4.
Code problem12.m solves this problem, and calls function formStiffness3D
frame.m, that computes the stiffness matrix of the 3D frame element.
%................................................................
% MATLAB codes for Finite Element Analysis
% problem12.m
% antonio ferreira 2008
8.3 First 3D frame example 105
1
2 3
4
10kN 15kN
4m
x 3m 3m
y
z
Fig. 8.1 A 3D frame example (problem12.m)
% clear memory
clear all
% E; modulus of elasticity
% I: second moment of area
% L: length of bar
E=210e6; A=0.02;
Iy=10e-5; Iz=20e-5; J=5e-5; G=84e6;
% generation of coordinates and connectivities
nodeCoordinates=[0 0 0; 3 0 0; 0 0 -3; 0 -4 0];
xx=nodeCoordinates(:,1);
yy=nodeCoordinates(:,2);
zz=nodeCoordinates(:,3);
elementNodes=[1 2;1 3;1 4];
numberNodes=size(nodeCoordinates,1);
numberElements=size(elementNodes,1);
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
% GDof: global number of degrees of freedom
GDof=6*numberNodes;
U=zeros(GDof,1);
force=zeros(GDof,1);
stiffness=zeros(GDof,GDof);
106 8 Analysis of 3D frames
%force vector
force(1)=-10;
force(3)=20;
% stiffness matrix
[stiffness]=...
formStiffness3Dframe(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,E,A,Iz,Iy,G,J);
% boundary conditions and solution
prescribedDof=[7:24];
% solution
displacements=solution(GDof,prescribedDof,stiffness,force)
% displacements
disp(’Displacements’)
jj=1:GDof; format long
f=[jj; displacements’];
fprintf(’node U\n’)
fprintf(’%3d %12.8f\n’,f)
Listing of formStiffness3Dframe.m:
function [stiffness]=...
formStiffness3Dframe(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,E,A,Iz,Iy,G,J);
stiffness=zeros(GDof);
% computation of the system stiffness matrix
for e=1:numberElements;
% elementDof: element degrees of freedom (Dof)
indice=elementNodes(e,:) ;
elementDof=[6*indice(1)-5 6*indice(1)-4 6*indice(1)-3 ...
6*indice(1)-2 6*indice(1)-1 6*indice(1)...
6*indice(2)-5 6*indice(2)-4 6*indice(2)-3 ...
6*indice(2)-2 6*indice(2)-1 6*indice(2)] ;
x1=nodeCoordinates(indice(1),1);
y1=nodeCoordinates(indice(1),2);
z1=nodeCoordinates(indice(1),3);
x2=nodeCoordinates(indice(2),1);
8.3 First 3D frame example 107
y2=nodeCoordinates(indice(2),2);
z2=nodeCoordinates(indice(2),3);
L = sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1) +...
(z2-z1)*(z2-z1));
k1 = E*A/L;
k2 = 12*E*Iz/(L*L*L);
k3 = 6*E*Iz/(L*L);
k4 = 4*E*Iz/L;
k5 = 2*E*Iz/L;
k6 = 12*E*Iy/(L*L*L);
k7 = 6*E*Iy/(L*L);
k8 = 4*E*Iy/L;
k9 = 2*E*Iy/L;
k10 = G*J/L;
a=[k1 0 0; 0 k2 0; 0 0 k6];
b=[ 0 0 0;0 0 k3; 0 -k7 0];
c=[k10 0 0;0 k8 0; 0 0 k4];
d=[-k10 0 0;0 k9 0;0 0 k5];
k = [a b -a b;b’ c b d; (-a)’ b’ a -b;b’ d’ (-b)’ c];
if x1 == x2 & y1 == y2
if z2 > z1
Lambda = [0 0 1 ; 0 1 0 ; -1 0 0];
else
Lambda = [0 0 -1 ; 0 1 0 ; 1 0 0];
end
else
CXx = (x2-x1)/L;
CYx = (y2-y1)/L;
CZx = (z2-z1)/L;
D = sqrt(CXx*CXx + CYx*CYx);
CXy = -CYx/D;
CYy = CXx/D;
CZy = 0;
CXz = -CXx*CZx/D;
CYz = -CYx*CZx/D;
CZz = D;
Lambda = [CXx CYx CZx ;CXy CYy CZy ;CXz CYz CZz];
end
108 8 Analysis of 3D frames
R = [Lambda zeros(3,9); zeros(3) Lambda zeros(3,6);
zeros(3,6) Lambda zeros(3);zeros(3,9) Lambda];
stiffness(elementDof,elementDof)=...
stiffness(elementDof,elementDof)+R’*k*R;
end
Results are listed as follows.
Displacements
node U
1 -0.00000705
2 -0.00000007
3 0.00001418
4 0.00000145
5 0.00000175
6 0.00000114
7 0.00000000
8 0.00000000
9 0.00000000
10 0.00000000
11 0.00000000
12 0.00000000
13 0.00000000
14 0.00000000
15 0.00000000
16 0.00000000
17 0.00000000
18 0.00000000
19 0.00000000
20 0.00000000
21 0.00000000
22 0.00000000
23 0.00000000
24 0.00000000
8.4 Second 3D frame example
The next 3D problem is illustrated in figure 8.2 and considers E = 210 GPa, G =
84 GPa, A = 2×10−2 m2, Iy = 10×10−5 m4, Iz = 20×10−5 m4, J = 5×10−5 m4.
The MATLAB code for this problem is problem13.m.
8.4 Second 3D frame example 109
1
3
4
2
15kN
5m
x 4m 4m
y
z
Fig. 8.2 A second 3D frame example (problem13.m)
%................................................................
% MATLAB codes for Finite Element Analysis
% problem13.m
% antonio ferreira 2008
% clear memory
clear all
% E; modulus of elasticity
% I: second moments of area
% J: polar moment of inertia
% G: shear modulus
% L: length of bar
E=210e6; A=0.02;
Iy=10e-5; Iz=20e-5; J=5e-5; G=84e6;
% generation of coordinates and connectivities
nodeCoordinates=[0 0 0;
0 0 4;
4 0 4;
4 0 0;
0 5 0;
0 5 4;
110 8 Analysis of 3D frames
4 5 4;
4 5 0;
];
xx=nodeCoordinates(:,1);
yy=nodeCoordinates(:,2);
zz=nodeCoordinates(:,3);
elementNodes=[1 5;2 6;3 7; 4 8; 5 6; 6 7; 7 8; 8 5];
numberNodes=size(nodeCoordinates,1);
numberElements=size(elementNodes,1);
% for structure:
% displacements: displacement vector
% force : force vector
% stiffness: stiffness matrix
% GDof: global number of degrees of freedom
GDof=6*numberNodes;
U=zeros(GDof,1);
force=zeros(GDof,1);
stiffness=zeros(GDof);
%force vector
force(37)=-15;
% calculation of the system stiffness matrix
% and force vector
% stiffness matrix
[stiffness]=...
formStiffness3Dframe(GDof,numberElements,...
elementNodes,numberNodes,nodeCoordinates,E,A,Iz,Iy,G,J);
% boundary conditions and solution
prescribedDof=[1:24];
% solution
displacements=solution(GDof,prescribedDof,stiffness,force);
% displacements
disp(’Displacements’)
jj=1:GDof; format long
f=[jj; displacements’];
fprintf(’node U\n’)
fprintf(’%3d %12.8f\n’,f)
8.4 Second 3D frame example 111
%drawing mesh and deformed shape
U=displacements;
clf
drawingMesh(nodeCoordinates+500*[U(1:6:6*numberNodes)...
U(2:6:6*numberNodes) U(3:6:6*numberNodes)],...
elementNodes,’L2’,’k.-’);
drawingMesh(nodeCoordinates,elementNodes,’L2’,’k--’);
Results are obtained as: