IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Finite Elements: Modelling methods, Implicit
and Explicit analysis
IFB
Modelling hints and methods
Brief introduction to the Implicit method: Element stiffness and
assembly for bars
The Explicit integration scheme
Comparing Implicit and Explicit methods
Exercises on Explicit integration for a simple bar element using
Fortran and/or VBA Excel:
1.
Loading via velocity and force control
2.
Timestep and stability
3.
Nodal damping
4.
Contact treatment
PAM-CRASH Bar example: Force, velocity and damped loading
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Finite Elements: Some modelling hints
Contents:
Meshing guidelines.
Creating meshes and models.
Grading of meshes.
Poor elements - try to avoid them.
Convergence - refined meshes for better results.
Application of symmetry to save CPU costs.
Application of loading.
Example.
2
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Some comments on meshing from NAFEMS
.
Extract from: NAFEMS (National Agency for Finite Element
Methods and Standards) A finite element primer
Clear meshing rules are difficult to define, much depends on experience and instinct.
Carefully study and examine the results check for meshing/analysis errors and look
for possible mesh improvements.
3
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
How a pre-processor creates a mesh and analysis model (1/2)
Generally most simple preprocessors build up a mesh as follows:
1. Points are specified which will be used as reference points for lines. They can also be
future locations for loading, boundary conditions, output of information, etc.
2. Lines (straight, arcs, spline, etc.) are generated from the points.
3. For 2D meshes surfaces are generated from/between the lines.
4. For 3D meshes volumes are generated from the surfaces.
5. The 2D surfaces, or 3D volumes are meshed using the
selected element types.
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
How a pre-processor creates a mesh and analysis model (2/2)
You may want to regrade the mesh until you feel it is appropriate for the geometry
and loading.
Then apply boundary conditions; make sure you have sufficient constraints to
suppress all rigid body modes of the structure.
Impose loadings (point loads, pressures, perhaps gravity).
Define geometry information (thickness, inertias...) and material data (Modulus,
Poissons ratio, plasticity law...).
t, E, .
OR CAD data that defines the
geometry is used with
automatic/manual meshing tools
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Most pre-processors (mesh generator
programs) are able to clean-up the
CAD data and generate an FE mesh
e.g. Visual Mesh can automatically
surface mesh this complex structure in
about 1 minute!
But the initial CAD data should be good.
Only limited control is possible for the
mesh sizes and element distributions
Final FE Mesh
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Good and bad elements
Try to avoid distorted elements the formulation of [k] and the results they
yield can be poor(er)
Shape
good
poor
Aspect ratio
Avoid poor aspect ratio
elements
Warping
Avoid warped shell elements
ETC.
7
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Grading of FE meshes
Incorrect way of joining meshes
Incompatible displacement fields
between elements due to free nodes
holes will appear in the structure!
A
Either remesh (e.g. as below) or use
specialised tied constraint options that are
available in most codes to overcome this
and tie elements row A to element row B.
Extra constraint conditions are imposed
that tie free nodes to other structural nodes
Mesh grading can be done by introducing new
elements and triangles. Be careful the new
(distorted) elements and triangles are not
particularly good, so try to avoid doing this
grading in critical areas !
8
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Grading meshes using interface elements
E.g. detailed model of a spotweld: The local model
uses a carefully graded mesh from the boundary to the
nugget. For the box beam a coarse mesh of shell
elements is used to transfer loading to the area of
interest; namely the spotweld (modelled as 3D solids).
Note: This study uses the explicit FE crash simulation
code PAM-CRASH.
Special constraint conditions tie the
coarse shells to the fine solids
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Examples of meshing at a stress concentration
Try to distribute elements according to stress gradients. In stress concentrations fine
elements are needed; in areas of low stress gradients a coarse mesh is sufficient to save
CPU/memory. In between grade the elements carefully.
Exactly the same logic applies for flow type problems (e.g. heat or fluids): Use fine
meshes in areas of high flow gradients and coarse meshes in areas of low gradients to
save CPU/memory costs.
High flow gradient
close to input point
(heat, liquid...)
High stresses and
stress gradients are
located at the root of
the notch (use a FINE
mesh)
Areas of low stress
gradients (use a
COARSE regular mesh)
Low flow gradients remote from input point
(heat, liquid...).
But beware of any potential gradients at
boundary conditions (= FINE meshes again).
10
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Examples of meshing in a stress concentration
ST RESS
CONT OURS OF SY
Examples of meshes for a notch analysis
using automatic and manual meshing
0
0.23828
0.47656
0.71484
0.953121
1.1914
1.42968
1.66796
Poor
Better
F 1.90624
2.14452
2.3828
2.62108
2.85936
3.09764
3.33592
3.5742
C
E
M ax 3.632 at Node 1
M in -0.1801 at Node 40
B
A
Good
Note:
E
Z
Modern CPUs are so powerful that many
elements may be used if necessary.
BEWARE results in stress concentration areas will
be dependent on the mesh use the same mesh
in parametric type studies.
F
G
H
I
JK
L
X
MN
Local mesh
around the
notch
Study the results and see if you are happy with
the meshing strategy.
11
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Other examples where fine meshes are needed
Explicit: Crash
E.g. buckling:
In buckling problems:
1. Explicit codes for crash and folding/buckling
2. Implicit codes for buckling (Eigenmode) analysis
The mesh must be fine enough to capture the required
physical buckling modes.
Implicit: Eigenmodes
Z
12
True
(physical)
folding
mode
Best folding
mode with
coarse mesh
Necessary element
sizes to capture
the physical folds
MODE
1
EIGENVALUE
53.5783
FREQUENCY
ERROR NORM
1.16497
0.606236E-12
0.458210E-12
57.7202
1.20916
-57.7215
0.000000E+00 0.477816E-12
-62.1734
0.000000E+00 0.404732E-12
62.1760
1.25496
0.926765E-12
63.1296
1.26455
0.190411E-10
A finer mesh will also predict
the same true folds (it will not
be mesh dependent)
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Convergence: The more elements the better, but the CPU/memory
costs increase fast
13
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Symmetry: Always try to make use of this to reduce model size
1. Planar symmetry
If a symmetric part of the model is
used for analysis then the correct
boundary conditions are essential.
2. Circular symmetry
Beware: In these cases the geometry
allows definition of symmetry planes,
but check the loading/restraints also are
valid for this symmetry
A quarter symmetry could also be
used here. But further reductions are
possible using boundary conditions
defined in a skew axis system.
A so-called skew axis to which local
boundary conditions act
14
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Applied loading: Point and distributed loads
Beware point loads at nodes give concentrated loads that may not be what you want. If
you want distributed loading you may need to superimpose elements to physically model
the effect you want (or use local element stiffening).
P
Additional beam/shell elements
may help to distribute unrealistic
local loads (for both forces &
moments).
Pressure loads are distributed to nodes; note that nodal loads depend on the area they
act over.
Distributed 100N/mm2
50N
100 100
50
50
100
100 75 50 50 50
25
Equivalent nodal loads
0.5mm
1mm
NB FE codes have options to
automatically distribute pressures to
equivalent nodal loads.
15
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Apply the loadings carefully to represent what you really intend !
LOAD CASE =
1
Loadcase 1
RESULTS FILE =
1
DISPLACEMENT
CONTOURS OF DY
-12.5311
-11.7479
-10.9647
-10.1815
-9.39829
-8.6151
-7.83191
-7.04872
-6.26553
-5.48234
-4.69915
-3.91595
-3.13276
-2.34957
-1.56638
-0.783191
Dmax =
12.53mm
Y
X
These two examples have the
same total applied loads, but
give very different maximum
deflections in the structure !
LOAD CASE
RESULTS FILE =
DISPLACEMENT
1
1
CONTOURS OF DY
-1,59696
-1,49715
-1,39734
-1,29753
-1,19772
-1,09791
-0,998098
Max 0.0000E+00 at Node 1
Min -12.53 at Node 310
Loadcase 1
-0,898288
-0,798478
-0,698669
Distributed loadings at
mid span of a beam
Dmax =
1.6mm
-0,598859
-0,499049
-0,399239
-0,299429
-0,19962
-0,0998098
Max 0.0000E+00 at Node 1
Min -1.597 at Node 20
16
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Simplified FE models
A correct model of the part and loadings/BC can be difficult to make and require alot of effort
to get right (+ CPU expensive); but a simplified model can provide quick and valuable
information to help guide a design.
The simplified model can help understand:
Failure mechanisms.
Means to increase performance (even if it does not
accurately predict maximum load).
Or other effects (e.g. points of max stress )
Use correct composite
models (layups and
materials) possibly
with failure
17
IFB
Use special delamination
elements, if this failure
mode is to be modelled
Critical areas must be
modelled properly
Try to model
supports accurately
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Simplified FE models: Results
It is possible to
use the
simplified
model for a
design of
experiments
approach
Fixture (failure loads)
Variation
18
Maximum
load [kN]
Change in failure
load
Reference model
1,68
2 extra 0 plies in web
2,07
+23%
2 extra 0 plies in flange
2,08
+23%
4 extra 0 plies in web
2,55
+51%
4 extra 0 plies in flange
2,503
+49%
2 extra web + 2 extra
flange (0)
2,38
+42%
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Some basics of the Implicit method:
1. Computing a stiffness matrix
2. Element assembly
3. Obtaining a solution
Comparison of the Implicit and Explicit methods
The Explicit method and Exercises:
1. Fortran and Excel code for the simple single element case
2. Velocity and imposed force loading
3. Timestep controls
4. Element with damping
5. Contact: Rigid wall and penalty force methods
The Explicit algorithim in PAM-CRASH
19
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Some basics of the Implicit method:
1. Computing a stiffness matrix
2. Element assembly
3. Obtaining a solution
Comparison of the Implicit and Explicit methods
The Explicit method and Exercises:
1. Fortran and Excel code for the simple single element case
2. Velocity and imposed force loading
3. Timestep controls
4. Element with damping
5. Contact: Rigid wall and penalty force methods
The Explicit algorithim in PAM-CRASH
20
10
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Consider the simple (1D) Bar element
Length L, Modulus E and sectional area A
Only axial loads Fi , Fj and axial displacements ui, uj are allowed
It is said to have 1 d.o.f. (degree of freedom) per node. Torsion bending and transverse
shear are not permissible for this element type.
Sign convention
+ positive forces and displacement (to the right)
Case 1: ui>0 and uj=0
Fj(=-Fi)
Fi = (EA/L)ui, Fj= -(EA/L)ui
Fi
u j= 0
ui
Case 2: ui=0 and uj>0
Fi(=-Fj)
Fi = -(EA/L)uj, Fj= (EA/L)uj
ui = 0
uj
Fi EA / L EA / L ui
=
F j EA / L EA / L u j
In matrix form
For the general case ui0 and uj0
Fi = (EA/L)*( ui-uj)
{F }= [K ]{d }
Fi k
=
F j k
Fj = (EA/L)*(-ui+uj)
21
Spring stiffness
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Fj
Stiffness matrix
k ui
k u j
Composites modelling:
Modelling and Finite Elements (continued)
Connectivity of two elements (1/2)
L1
1
L2
E1, A1
E2, A2
k1 = E1A1/L1
k2 = E2A2/L2
u2
u1
2 ,
F2
Element 1
k1
k
1
Element 2
k2
k
2
2
F1=- f11
Force equilibrium
f11
u3
u 1 , F1
F1
3,
F3
dof
k1 u1 f
=
k1 u 2 f
1
1
1
2
Internal
nodal forces
dof
k 2 u2 f 22
=
k 2 u3 f 32
22
11
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Connectivity of two elements (2/2)
The stiffness matrices for each element are:
u1
u2
u1
u3
u2
u3
d.o.f.
0
0 0
0
k
2
k 2
Element 2 only
0 k 2 k 2
k 1 k 1 0
Element 1 only k 1 k 1 0
0
0 0
ASSEMBLY
Global stiffness matrix of the assembly
[K ]{U } = {F }
1
k1
k1
k k + k
1
2
1
k2
0
23
IFB
Force equilibrium at nodes requires:
Nodes 1: forces=0 -f11
dof
= F1
Nodes 2: forces=0 -f21 f22 = F2
0 u1 F1
k 2 u2 = F2
k 2 u3 F3
Nodes 3: forces=0 -f32
= F3
External loads
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Solution for unknowns (application of Boundary conditions)
Matrix [K] is singular (it can have any number of solutions since the structure is unsupported
and can undergo arbitrary Rigid Body Motion (RBM). At least one support must be imposed.
To suppress RBM set, for example, u3=0: (u1=0 or u2=0 would also suppress RBM).
k1
k1
k k + k
1
2
1
0 0 k 20
Operations done
to enforce u3 =0
= -P
0 0 u1 F1
k 20 u2 = F2
=0
0
k 2 1 u3 F3
u2 = - P / k2
1
P
2
1
3 dof
2
u1 = -P / k1 - P / k2
u3 = 0
Stresses and strains could also be found using:
1 = E1 1 = E1 (u2-u1) /L1 etc.
Note:
24
Ditto for 2
Of course changing the position of the constraint to suppress RBM changes
the problem and the results.
Exercise: Impose the constraint at d.o.f 2; what are the expressions for displacement
and strain? (hopefully 2 = 0).
12
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Structure of a simple Implicit FE program
1. Read in the nodal coordinates to a node array
2. Read in the element connectivity information
3. Read in material data for each element (usually done in groups - Material cards)
4. Compute the stiffness matrix [k] for each element and assemble to the global stiffness
matrix [K]
5. Assemble all loads to the load vector {P}
6. Apply boundary conditions to the Global stiffness matrix (making it positive definite) and
also to the load vector {P}
7. Invert the Global stiffness matrix [K] giving [K]-1
8. Compute the nodal displacements {u} from {u} = [K]-1 * {P}
9. Compute element strains from nodal displacements
10. Compute element stresses from strains with chosen material law
11. Check the solution: e.g results and deformations, load/reaction balance!
[k2]
[k1]
[k3]
[K]
25
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
After assembly the stiffness matrix is inverted for nodal displacements
{P}
{P} = [K ]{u}
{u} = [K ]1{P}
{u}
For STATIC linear
problems [K] = constant
{P}
constant
[K ] =requires
a simple
linear solution
Linear versus Non-linear
Most problems are obvious.
But for some it is less clear
and a comparative analysis
(Linear versus Non-linear) is
needed to check the required
analysis type.
26
[K ]
constant
requires a non-linear
(iterative) solution;
e.g. interative
Newton-Raphson
{u}
13
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
P
Linear versus Non-linear FE analysis
Force vs Displacement for the 4 point bending can be
predicted using beam theory (NB this is based on initial
geometry and assumes small displacements).
Linear FE results would be similar.
Both only have a limited range of validity.
P=
24dEI
a (3L2 4a 2 )
Linear FE analysis
Analytical formulae
Point A (limit of
linear assumption)
varies depending
on the problem:
loading
geometry
materials.
Explicit FE analysis
Non-linear FE analysis
Linear analysis
region
Non-linear
analysis region
27
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Some basics of the Implicit method:
1. Computing a stiffness matrix
2. Element assembly
3. Obtaining a solution
Comparison of the Implicit and Explicit methods
The Explicit method and Exercises:
1. Fortran and Excel code for the simple single element case
2. Velocity and imposed force loading
3. Timestep controls
4. Element with damping
5. Contact: Rigid wall and penalty force methods
The Explicit algorithim in PAM-CRASH
28
14
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
IFB
Explicit method
Implicit method
f(t),x, x& , &x&
{P} = [K ]{ }
Load
Stiffness
Displacements
(1) &x&n = m1(fn kxn )
(2) x& n+1/2 = x& n1/2 + tn&x&n
m&x& + cx& + kx = f(t)
(3)
m&x&n + kxn = fn (t)
x n+1 = xn + tn+1/2x& n+1/2
Advantages:
Advantages:
CPU efficient and robust
Static problems
Very large model sizes possible
Best for mildly non-linear
problems
Highly non-linear materials
Large deformations
CPU efficient for contact problems
29
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Typical commercial codes: Names, importance and applications
Implicit FE
Main codes
1. NASTRAN (large linear
structures)
2. ABAQUS (Non linear)
3. PAM-Implicit
4. MARC (Non linear)
5. ANSYS (general purpose)
6. IDEAS (meshing + FE)....
Literature
Industry
Car companies
Applications
30
Explicit FE
1.
2.
3.
4.
5.
6.
90%
1.
2.
3.
4.
5.
6.
DYNA3D
PAMCRASH
PAMSTAMP
PAMFORM
RADIOSS
ABAQUS explicit
Others
1.
2.
3.
4.
5.
6.
CFD codes
Vehicle dynamics
Welding simulation
Casting simulation
Vibro-acoustics
....
10%
90%
10%
50%
50% (due to importance of crash/safety/etc..)
Stress Analysis
Stiffness
Eigenvalue (vibrations)
Flow problems; eg
Heat, Magnetism.
Fatigue
.
1.
2.
3.
4.
5.
Crash
Safety
Stamping
Biomechanics
..
15
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Explicit FE codes are particularly strong on:
Large deformations and extreme material nonlinearities
Handling buckling
Efficient CPU times
Large structures
Contact (internal and to other structures)
31
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Typical Explicit FE examples
32
16
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Time integration: The Explicit integration scheme
m&x& + kx = f (t )
We start with the dynamic equation of equilibrium for the node
k
known
starting
point
&x&(t )
f (t ), x , x&, &x&
m&x&n + kxn = f n (t )
&x&n
&x&n 1
(1)
&x&n = m 1 ( f n kxn )
( 2)
x& n +1/ 2 = x& n 1/ 2 + t n &x&n
(3)
xn +1 = xn + t n +1/ 2 x&n +1/ 2
x& (t )
x&n +1 / 2
x&n 1 / 2
x (t )
tn
tn +1 / 2
xn +1
xn
tn 1
tn 1 / 2
tn
tn +1 / 2
tn +1
33
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Stable time step for the explicit scheme (Courant criteria)
m = AL/2
t
k = EA/L
L, E, ,
c= E
bar
spring-mass
L
c
is the speed of sound in
the material
is the Traversal time
Small elements (small L) are bad for the timestep (small) >> High CPU cost
Similarly, stiff (large E) or light elements (low ) also reduce the timestep
The timestep criteria ensures that the shock wave (information) does not pass more than
one element in one timestep
Note: The worst (smallest timestep) computed for all elements in the structure is
used for analysis! Unless specialised sub-cycling techniques are activated.
34
17
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Time integration methods for dynamic problems using Implicit integration
is also possible
m&x& + kx = f (t )
Known at n
m&x&n+1 + kxn+1 = f n+1 (t )
To find
( x& x& )
&x&n+1 = n+1 n
tn
f (t ), x, x&, x&&
&x&(t )
x&n+1 =
( xn+1 xn )
(x x )
x&n = n n1
tn
tn
&x&n+1 =
( xn+1 2xn + xn1 )
tn2
&x&n +1
known
&x&n
x& (t )
x&n +1
x&n
m
(2xn xn1 )
tn2
m
2 + k
tn
f n+1 +
x (t )
tn
xn +1
xn
tn 1
tn
tn +1
xn+1 =
t
k is function a function of xn+1
(iterative solution required)
35
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Implicit versus Explicit (for transient solutions)
Explicit methods use direct integration of the dynamic equations of motion to advance the
solution,
x n+1 =
f(
x n, x& n , &x&n , x n-1.......)
i.e. new information at time (n+1)T can be explicitly found from
previous known information at time nT.
A node-by-node solution is possible provided the stable timestep is used.
The solution is said to be conditionally stable (i.e. its stable provided a
valid timestep is used).
Implicit methods use solution methods having the form,
x n+1 =
36
f(
x& n+1, &x&n+1, x n , x& n, &x&n,.......)
i.e. new information at (n+1)T requires information of velocities and
accelerations at (n+1)T, which are not directly not known (they are a
function of xn+1 and can only be found implicitly).
xn+1 is determined using iterative convergence techniques.
The solution is unconditionaly stable (theoetically any timestep can be used).
18
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Comparison: Explicit versus Implicit
Number of time steps:
Although unconditionally stable, Implicit methods needs many time
steps in order to correctly trace the physical phenomena.
Explicit analysis requires a small time step. This leads to many, but
CPU cheap solutions: Contact, material and geometrical nonlinearities
are easily handled.
Equations to be solved:
Implicit analysis requires CPU expensive matrix inversion. The
solution of non-linear problems requires iterative solution strategies.
Explicit analysis requires no iterations and no matrix inversion.
Static problems:
Implicit methods can solve static problems.
Explicit methods can only solve quasi-static problems by careful
choice of slow and ramped loading and maybe some added damping.
37
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Propagation of shock wave information in a 1D explicit integration system
Node 1
Mass M
2
2M
3
2M
L o , F int = 0
4
M
L o , Fint = 0
L o , Fint = 0
Pext
A n = (Pext- Fint)/M
V n+ =V n-+A n * T
X n+1=X n + V n+ * T
n=0
An = 0
V n+ =0
X n+1= 0
An = 0
V n+ =0
X n+1= 0
L n+1 = Ln- X n+1
Fn+1 = K *X n+1
Pext
F int = (Fn+1)
A n = (Pext- Fint)/M
V n+ =V n-+A n * T
X n+1= Xn+ V n+ * T
F int = 0
=0
A n =(Pext - Fint)/ 2M
V n+ =V n- +A n * T
X n+1= Xn +V n+ * T
An = 0
V n+ =0
X n+1= 0
n=1
L n+1 = Ln - X n+1
Fn+1 = K *X n+1
Pext
F int = (Fn+1)
L n+1 = Ln - X n+1
Fn+1 = K *X n+1
Fint = (Fn+1)
Fint = 0
ETC
38
19
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Propagation of information in a 1D system and consequences of the stable
timestep
Node1
Axial force
Bar 3
Bar 2
Bar 1
Node 2
acceleration
velocity
Timestep
displacement
39
IFB
T=0
T=1
T=2
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Some basics of the Implicit method:
1. Computing a stiffness matrix
2. Element assembly
3. Obtaining a solution
Comparison of the Implicit and Explicit methods
The Explicit method and Exercises:
1. Fortran and Excel code for the simple single element case
2. Velocity and imposed force loading
3. Timestep controls
4. Element with damping
5. Contact: Rigid wall and penalty force methods
The Explicit algorithim in PAM-CRASH
40
20
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
The basic problem
A simple elastic spring, fixed at the base with velocity (or force) loading at the free end:
Velocity (or force) loading
Free end
Length = 100mm
Added mass = 50kg
E = 70 kN/mm2
Area = 100 mm2
Fixed end
Consistent units used: kg, mm, msec, kN
41
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 1. Explicit integration for the
simple elastic spring (Fortran version)
Program ElasticSpring
c
c... opening output files
Open(unit=3,file='Force_TimeHistory.xls')
Open(unit=4,file='Length_TimeHistory.xls')
c
c... initialization ( bar and integration data)
Alength_O = 100
Initialisation
Area
= 100.
and
Emod
= 70.0
meshing
AddedMass = 50.
c
TimeTotal = 0.0
Timestep = 0.001
Tend
= 10.
Fint = 0.
PreVold = -2.5
processing
Uold = 0.0
Fext = 0.0
c
42
Composites modelling:
Modelling and Finite Elements (continued)
c... start of integration: update acceleration,
c... velocities and displacement (explicit)
c
100 Fres = -Fint
Acel = Fres / AddedMass
Vnew = Vold + Acel * Timestep
Vold = Vnew
Unew = Uold + Vnew * Timestep
Uold = Unew
Solution
Alength = Alength_O + Unew
Strain = Unew / Alength_O
c
F=(EA/L)*
Fint = Emod * Strain * Area
TimeTotal = TimeTotal + Timestep
c
c... outputs (bar force and length time histories)
write(3,21)TimeTotal,Fres
write(4,21)TimeTotal,Alength
21 Format(E20.8,',',E20.8)
c
If(TimeTotal.gt.Tend)go to 200
go to 100
Postc
Processing
200 stop
end
21
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Exercises
There are seven simple Excel VBA codes I have prepared:
Ex 1. Explicit_SimpleSpring.xlsm
Explicit spring with velocity control
Ex 2. Explicit_SimpleSpring_Force. xlsm
Explicit spring with force control
Ex 3. Explicit_SimpleSpring_ForceTimestep.xlsm
Explicit spring to study stable timestep
Ex 4. Explicit_SimpleSpring_ForceDamp.xlsm
Explicit spring with nodal damping
Ex 5. Explicit_SimpleSpring_Contact.xlsm
Explicit spring with penalty contact
Ex 6. Explicit_SimpleSpring_ContactEnergies.xlsm
Explicit spring with energy balance check
Ex 7. Explicit_SimpleSpring_ContactEnergies_Timestep.xlsm
Explicit spring with energy balance check and stable/unstable timestep
43
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Exercises: Getting started
Try to:
Copy to a working directory and open the first example (Ex.1)
Study and understand it (see next slide for hints on using VBA Excel)
Compile, run and study the results (see additional handout notes)
Try to use this (Ex. 1) as a basis to do the other examples on the following pages:
1. Apply Force loading (Ex. 2)
2. Study timestep stability criteria for the explicit method (Ex. 3)
3. Apply nodal damping to get a quasi-static solution (Ex. 4)
Modifying Example 1 to do Examples 2,3,4; if you get stuck I have also supplied the
Excel VBA codes for all of these.
Try example 5, 6 and 7 yourself which cover contact analysis and energy balance for
stable and unstable solutions. Hints on doing these are given in the following slides.
44
22
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
How operate Microsoft Excel VBA
1. Start Micosoft Excel and read in the required file
(e.g. Ex 1. = Explicit_SimpleSpring.xlsm)
2. Press Alt with F11 to enter the programming mode
3. The coding is under Module; open this. Press the
right mouse key and code anzeigen to see the code
4. If either sicherheitswarnung appears activate it so the
macros (coding) can be seen/changed !
5. Any lines of the code can now be modified
6. To run the (modified) code press
7. If there are programming errors these will be
described and error line highlighted
8. To see the results go back to the spredsheet by
pressing Alt with F11
9. Note Alt with F11 toggles between code and
spreadsheet
45
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions: Ex 1. The simple (initial velocity) case
46
23
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 1. Results for the simple elastic spring (VBA Excel)
1. After running the code as previously described you should get the following spreadsheet
Time
101,5
101
100,5
100
99,5
99
98,5
47
IFB
0,01
0,02
0,03
0,04
0,05
0,06
0,07
0,08
0,09
0,1
0,11
0,12
0,13
0,14
0,15
0,16
0,17
0,18
0,19
0,2
0,21
0,22
0,23
0,24
0,25
0,26
0,27
0,28
0,29
0
0,3
Bar length
99,975
99,9500088
99,925035
99,9000875
99,8751749
99,8503061
99,8254896
99,8007342
99,7760486
99,7514413
99,7269211
99,7024964
99,6781758
99,6539679
99,6298811
99,6059238
99,5821045
99,5584314
99,5349129
99,5115571
99,4883723
99,4653666
99,442548
99,4199245
99,397504
99,3752944
99,3533035
99,3315389
99,3100082
2
99,2887191
Time
0,01
0,02
0,03
0,04
0,05
0,06
0,07
0,08
0,09
0,1
0,11
0,12
0,13
0,14
0,15
0,16
Bar length
0,17
0,18
0,19
0,2
0,21
0,22
0,23
0,24
0,25
0,26
0,27
0,28
0,29
6
0,3
Res. Force
0
1,75
3,4993875
5,24755021
6,99387629
8,7377545
150
10,4785745
12,215727
100
13,948604
50
15,676599
17,3991072
19,11552560
20,8252537 0
-50
22,5276929
24,2222474
-100
25,9083242
27,585333
-150
29,2526869
30,9098025
32,5560995
34,191002
35,8139376
37,4243383
39,0216405
40,6052851
42,1747179
43,7293895
45,2687559
846,7922782
10
48,2994231
Bar (nodal
force) versus
time
Res. Force
10
12
Mass = 50kg
Bar length
versus time
12
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 2. Simple spring with added end load of 1000 (kN)
Modify Example 1
Force = 1000kN
The added force creates a permanent end displacement which
should oscillate around ca. 115mm length
Mass = 50kg
There is no damping in the system so the oscillations continue.
This is typical for dynamic explicit integration and many short
duration impact problems do not consider damping; e.g. car crash
single element length time history
1,40E+02
Mean
length ca.
115mm
1,20E+02
length
1,00E+02
8,00E+01
Extract of modified code for additional
force:
Date
6,00E+01
4,00E+01
2,00E+01
0,00E+00
0,00E+00 2,00E+00 4,00E+00 6,00E+00 8,00E+00 1,00E+01 1,20E+01
time
48
24
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 3. Simple spring and
timestep stability study
Try analyses using:
Modify example 2.
2. Tfac = 0.9 (default PAM- safety factor)
Program and use the equation
for stable timestep, instead of
the fixed value.
3. Tfac = 1.1-1.5 (this should be poor/wrong). You
will need to increase Tend = 100; beware it may
also abort the Solution!!
1. Tfac = 0.001
1,40E+02
T = Tfac * Tnodal
1,20E+02
1,00E+02
Tfac =
Dat
0.001 e
8,00E+01
Tfac = 1.0
6,00E+01
4,00E+01
If Tfac 1.0 the solution is
stable and results are similar
1,60E+02
1,40E+02
2,00E+01
1,20E+02
1,00E+02
0,00E+00
0,00E+00 2,00E+00 4,00E+00
8,00E+016,00E+00 8,00E+00 1,00E+01 1,20E+01
Tfac = 1.1
6,00E+01
4,00E+01
2,00E+01
0,00E+00
0,00E+00
49
IFB
5,00E+00
1,00E+01
1,50E+01
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 4. Apply nodal damping
Modify example 2.
Introduce an additional
(resistive) damping force
(Fdamp) which is added to the
nodal forces (= Fdamp - Fint).
You will need the period of
oscillation (T) for qcrit which is
obtained from a pre-analysis.
= 100 114.2 = 14.2mm
1,18E+02
Typically a percent of critical
damping (0.1-0.3), i.e. 10-30 %
is applied.
Note the solution quickly damps
to a stable solution.
1,16E+02
1,14E+02
1,12E+02
1,10E+02
IMPLICIT solution:
1,08E+02
In practice damping should not
be applied in the loading phase.
It should only be switched on
once the maximum deformation
are imposed.
50
Date
F=k
1,06E+02
1000 = (EA/L)
1,04E+02
= 14.2mm
1,02E+02
1,00E+02
9,80E+01
0,00E+00
2,00E+00
4,00E+00
6,00E+00
8,00E+00
1,00E+01
1,20E+01
25
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Ex 4. Nodal damping: Some things you can try
A
Critical damping
Undamped solution
Damped solution (e.g. 20%)
Overdamping
1. Try critical damping 1.0 * qcrit
2. Try over damping (say 2.0 * qcrit)
3. Try under damping (say 0.2 * qcrit)
4. The usual (safe) way is to apply sub-critical damping only when the steady state
deformations are reached (after A) for about 2 cycles and then switch it off to see if the
structure starts moving again (i.e. This checks if the steady state solution been reached)
51
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Contact: Explicit FE codes (mostly) use two types of contacts
1. Contact at a rigid surface/wall
In this case nodal velocities and accelerations are set to
zero if contact occurs and the velocity is in the direction of
the wall.
Note the usual explicit accelerations and velocities are
computed. The Rigid wall condition checked and if contact
occurs the values are overwritten.
2. Contact between deformable bodies (and internal self
contact)
In this case accelerations and velocities cannot simply be
adjusted. Instead temporary penalty forces are introduced
to resist the contact. Contact with separation is possible.
Penalty force = Intrusion dist. * material stiffness *
empirical factor (e.g. 0.1-10)
Ex 5. Contact: This considers the simple spring problem with a penalty contact force
(temporary spring) to stop displacements on one side on the
oscillations (see next slide).
52
26
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Ex 5 Contact: Description of problem
Ca. 1mm
0.5mm
Mass = 10 kg
Ca. 1mm
Add additional spring force when displacements Unew
> say 0.5mm.
Fcont = factor (e.g. 10.)*Emod*(Unew -0.5)
Ex 5. Contact: You could try this with the simple spring example (Ex.1 is easiest).
Impose penalty force treatment at Unew = say 0.5mm. This is simply
an extra (stiff) spring that is only active when the displacement of the
free node is greater than 0.5mm.
53
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Materials laws
1
4
.
2
5
1.
2.
3.
4.
5.
Simple elastic law (programmed)
Non-linear elastic load with unloading (e.g. some rubbers)
Elasto plastic with unloading (e.g. most simple metals)
Strain rate effects to plasticity laws (most metals and plastics)
Load with Hysteresis unloading (e.g. foams)
54
27
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Energy computations
Numerical stability in the explicit method must be checked: check the general results and
especially from the energy balance.
Numerical instability is often caused by violation
of the timestep criteria. It may lead to:
Rapid generation of nonsense results (e.g.
large numbers) and failure
Or local instability causing e.g. artificial
plasticity; the analysis can then regain
stability. This can be identified by
generation of energy.
Example energies in the .out and results files
for PAM-CRASH. Overall structural and energies
per material are available.
55
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Energy computations
For a simple impact problem; such as car crash the initial Kinetic energy is converted to
Internal (plastic or damage) energy.
Indication of numerical instability
KE =
*Mass*
vel^2
Total Energy (KE +
IE = constant)
Internal Energy (IE)
Kinetic Energy (KE)
For problems where energy is fed into the system e.g. Application of external loading the
constant energy check is more difficult; but still look for sudden energy increases.
Contact can be responsible for such numerical instability.
Ex 6. Energy balance:
Compute and check energy balance for the contact example (Ex 5.). You will have to
compute energies for the bar (internal), contact spring and kinetic energy. Then add these
together to check the total is constant.
Ex 7. Contact and stability:
Combine Ex. 6 with the stable timestep criteria (Ex. 3) and try Tfac = 0.01, 0.1 and 1.1; take
56 a look at the results and the energy balance.
28
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Comparisons with PAM-CRASH explicit
Essentially the same operations are done in PAM-CRASH
The masses are assembled (lumped) at the nodes giving a 6*6
mass matrix.
1. The accelerations (3 translational + 3 rotational) are computed
2. The velocities and displacements are found (BCs etc. applied)
3. The element strains are computed from nodal displacements
4. The elements stresses are computed from the material law
5. Nodal forces Fint (balancing elements stresses) are found
6. Next cycle
57
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
mx
my
mz
0
Ix
Iy
Iz
The translational mass
and rotational inertias are
in different frames:
Translational Global fr.
Interia Principle fr.
This keeps the matrix
diagonal (= no matrix
inversion)!
Composites modelling:
Modelling and Finite Elements (continued)
PAM-CRASH FE models: Velocity, force and damping examples
A simple elastic spring, fixed the base with velocity (or force) loading at free end:
Velocity (or force) loading
Free end
Length = 100mm
PAM-CRASH example datasets for force,
velocity and damped loading:
Fixed end
Added mass = 50kg
E = 70 kN/mm2
Area = 100 mm2
Consistent units used:
kg, mm, msec, kN
Example results for
nodal displacement
(z) for the case of
damped loading
58
29
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions to VBA Excel
Explicit examples
59
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions: Ex 2. The simple (initial velocity and force) case
..
New external load
..
60
30
IFB
Composites modelling:
Modelling and Finite Elements (continued)
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Solutions: Ex 3. Timestep stability study
..
Still set but not used
..
61
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions: Ex 4. Application of nodal damping
..
New nodal force (see below)
T = 3 msec with 20% of critical damping
..
62
31
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions: Ex 5. Penalty contact force
..
Reduced to lower mass effect
Position of contact (Unew>0.5mm), Try
scale factors 1., 10., 100.
..
63
IFB
A.K. Pickett, 2013-2014
Institut fr Flugzeugbau, University of Stuttgart
Composites modelling:
Modelling and Finite Elements (continued)
Solutions: Ex 6. Check on energy balance
..
New code at the top for
output of headers
..
..
Additional code for
computation and output of
energies
64
..
32
