Hello, dear friend, you can consult us at any time if you have any questions, add WeChat: zz-x2580
Programming FEM for Stress
Analysis using MATLAB
CVEN9820
Computational Structural Mechanics
Finite Element Program
All the codes are stored inside the folder
CVEN4304.MTLB.CODE
• Readability takes precedence over
code efficiency and conciseness
• Read the file inside the folder
DOCUMENTATION to find out what
are the files that are stored in each
folder and subfolder
• Study the Matlab functions in the
folder SEEPAGE2D. You are
expected to port the code to stress
analysis
Functions for Stress Analysis
Tasks for assignment uses the folder ELAS2D:
• Save any new files you create in the subfolder CODES
• The subfolder TEMPLATES contain three files
– Elasticity_Template_01.m
– Elasticity_Template_02.m
– CVE4304FEASolverElast2D.m
Use these templates to aid you in completing the
assignment
Use of Templates
Matlab paths:
The path is set up at the beginning of the following files
Elasticity_Template_01.m
Elasticity_Template_02.m
% SET THE PATHS TO THE FOLDERS
addpath(genpath('D:\FEM\CVEN4304.MTLB.CODE'));
Use of Templates
Elasticity_Template_01.m
Manual input similar to HANDIN1.m
Input mesh details, material properties and boundary
conditions for Problem 1
Format for input is similar to that in seepage analysis (see
HANDIN1.m or Section 3 in the documentation file)
Boundary conditions are slightly different for stress analysis
Use of Templates
Elasticity_Template_01.m
Boundary conditions: Surface traction
Boundary conditions: Applied displacement
BC_trac = [NODE 1, NODE 2, X-TRACTION, Y-TRACTION;…];
BC_disp = [NODE NUMBER, DIRECTION, VALUE OF DISPLACEMENT;…];
Use of Templates
Elasticity_Template_02.m
Automatic mesh generation similar to HANDIN1_AUTO.m
Input mesh details, material properties and boundary
conditions for Problem 2
Format for input is similar to seepage analysis (see
HANDIN1_AUTO.m or Section 4 in the documentation file)
Boundary conditions are slightly different for stress analysis
Use of Templates
Both Elasticity_Template_01.m and
Elasticity_Template_02.m calls the function
CVE4304FEASolverElast2D
Hints:
You must write your finite element code for stress
analysis inside this function
Use the structure of CVEN4304SolverSeepage2D for
seepage analysis to guide you
Automatic Mesh Generation
Example in HANDIN1_AUTO.m
• Requires the following information
– Mesh size
– Material models
– Geometry
– Boundary conditions as applied on the geometry
• Automatically generates the variables coord,
EleConnTbl, BC_v and BC_h required to perform
a seepage analysis
Automatic Mesh Generation
Mesh size
Controls the size of the triangular elements that
are created by the mesh generator
hdata.hmax = …
Material models
mat(1).E = [];
mat(1).pos = [];
mat(1).plane = [];
Automatic Mesh Generation
Geometry: KEYPOINTS
Coordinates of vertices defining the structure
Input in Matlab
POINT=[X-COORD Y-COORD;…]
Automatic Mesh Generation
Geometry: LINES
Connectivity of KEYPOINTS that define the structure
Input in Matlab
LINE=[KYPT1 KYPT2;…]
Automatic Mesh Generation
Geometry: FACES
Connectivity of LINES that define the structure
Input in Matlab
FACE{1}=[LINE1 LINE2
LINE3… ; …]
MaterialInFace = […]
Note: See documentation for
creating geometries with more
than one face
Automatic Mesh Generation
Boundary conditions (BC) on geometry:
Surface traction
GivenTracOnLine = [LINE NO, X-VALUE, Y-VALUE;...];
Displacement
GivenDispOnLine = [LINE NO, DIRECTION, VALUE;...];
Automatic Mesh Generation
Transfer BC from geometry
to mesh:
Surface traction (use EdgesOnLine):
BC_trac=[];
if ~isempty(GivenTracOnLine)
for iTrac=1:size(GivenTracOnLine,1);
edges = EdgesOnLine{GivenTracOnLine(iTrac,1)};
BC_trac = [BC_trac; edges ones(length(edges),1)*...
GivenTracOnLine(iTrac,2:3)];
end
end
Your task: convert surface tractions to nodal forces
9
12
13
1 152
16
17
3
4
7
5
1
2
6
7
14
89
3 8
10
11 12
11
13
6
14
105
15
16
17
4
18
19
Automatic Mesh Generation
Transfer BC from geometry
to mesh:
Displacement (use NodesOnLine):
BC_disp=[];
if ~isempty(GivenDispOnLine)
for iDisp=1:size(GivenDispOnLine,1);
nodes = NodesOnLine{GivenDispOnLine(iDisp,1)};
BC_disp = [BC_disp; nodes
ones(size(nodes,1),1)*GivenDispOnLine(iDisp,2:3)];
end
end
9
12
13
1 152
16
17
3
4
7
5
1
2
6
7
14
89
3 8
10
11 12
11
13
6
14
105
15
16
17
4
18
19
Assembly of global stiffness matrix
• Numbering of DOFs in element equation = { }
• Numbering of DOFs in global equation = {}
Node 1 2 3 … i …
DOF 1 1 2 2 3 3 … i i …
Index 1 2 3 4 5 6 … (2 − 1) 2 …
Node 1 2 3
DOF 1
()
1
() 2() 2() 3() 3()
Index 1 2 3 4 5 6
Assembly of global stiffness matrix
• Element connectivity
Local Node 1 2 3
Local DOF 1
()
1
() 2() 2() 3() 3()
Local Index 1 2 3 4 5 6
… … … … … … …
El
e
7 Global Node
2 7 9
Global DOF 2 2 7 7 9 9
Global Index 3 4 13 14 17 18
9
12
13
1 152
16
17
3
4
7
5
1
2
6
7
14
89
3 8
10
11 12
11
13
6
14
105
15
16
17
4
18
19
Assembly of global stiffness matrix
Matlab code
stff = zeros(2*nNode,2*nNode);
for ie = 1:nEle %loop over elements
EleNode = EleConnTbl(ie,:);
…
EleStff = …; %element stiffness matrix
% DOF to DOF connectivity
EleDOF = reshape( [ 2*EleNode-1; 2*EleNode], [],1 );
% assemble global stiffness matrix
stff( EleDOF, EleDOF ) = stff( EleDOF, EleDOF ) + EleStff
