Transcript When a good agreement is not ehough

Solution of impact tasks,
assessment of result reliability
M. Okrouhlík
Institute of Thermomechanics
Scope of the lecture
• What is a ‘good agreement’ in dynamical
transient analysis in solid continuum mechanics
• Vehicle for comparison
• Continuum model and its limits
• Experimental and FE analysis
• Assessment of agreement quality
• Synergy of FE and experimental analyses
• Conclusions
A good agreement
of experimental and FE results
2
x 10
LOC1 - exp vs. FE 3D, NM
-4
G1R90L, mesh1, medium striker for all figures
0
-2
experiment
FE analysis, filtered
urmas data march 2007
-4
0
0.5
1
1.5
2
2.5
3
3.5
x 10
axial strain
2
x 10
-4
LOC2 - exp vs. FE 3D, NM
-4
0
experiment
FE analysis, filtered
-2
-4
0
2
0.5
x 10
1
-4
1.5
2
time
LOC3 - exp vs. FE 3D, NM
2.5
3
3.5
x 10
-4
0
-2
-4
0
dt = 1e-7; sr = 1/dt; nf = sr/2 = 5000000; rel cut off = 1/50 = 0.02
abs cut off = nf*rel cut off = 100000
f order = 2; [b,a] = butter(f order,rel cut off)
0.5
1
1.5
2
time
experiment
FE analysis, filtered
2.5
3
3.5
x 10
Experimental and FE treatment are limited by cut-off frequencies, which
are generally different – in this case the filter applied on FE results has
the same frequency limit as raw experimental data.
-4
Impact loading of a tube
with four spiral slots
Dimensions and points of interest
3D elements and node numbering
Location
for mesh1
Nodes in unfold location
Strain gauges
Strain
gauges
Nodes in one layer
Geometry, material and elements
G1R90L_NM ...
four slots with 15 deg. indentation, axial length = 17.28 mm,
slots and mesh rotate 90 deg
Geometry and elements
first part of the tube = 800 mm
second part, ie. axial length of slots = 17.28 mm
third part of the tube = 800 mm
total lenght of tube = 1617.28 mm
outer to inner diameters D/d = 22/16 mm
number of spiral slots = 4
800 layers, 144 elem in one layer
18 layers, 120 elem in one layer
800 layers, 144 elem in one layer
1618 total number of layers
Element properties
trilinear bricks elements, 8 nodes, Gauss order quadrature NG = 3
consistent mass matrix for NM (Newmark)
diagonal mass matrix for CD (central differences)
number of elements
232 560
number of nodes
310 576
number of degrees of freedom
931 728
max. front width
606
Material properties
Young modulus
Poisson's ratio
density
= 2.05e11 Pa
= 0.24
= 7800 kg/m^3
Loading and FE technology
Loading
Vertical arrangement of the loading is assumed.
In experiment the upper face of the tube is loaded by a vertically falling striker, which has been
released from a certain height.
Lower face of the tube is fixed. In FE model only axial displacements of the lower face are constrained.
The striker is made of the same material as the tube, has the same outer and inner diameters.
For computational purposes the upper side of the tube is loaded by uniform pressure,
whose time dependence is given by a rectangular pulse.
The loading pressure corresponds to the height from which the striker was released.
The time of the pulse corresponds to the length of the striker.
In this case the loading pressure is 88.5198 Mpa, which corresponds to the height of h = 1m.
Velocity of the striker just before the impact is sqrt(2*g*h) = 4.42944 m/s.
Material particle velocity immediately after the impact is v = 0.5*sqrt(2*g*h) = 2.2147 m/s.
The time length of the pulse is 15.6 microsec,
which corresponds to a medium length striker, Ls = 40 mm.
Pressure is evaluated from p = E*v/c0, with c0 = sqrt(E/ro).
Input energy from the striker is mgh = 0.548086 J. m = 05587 kg.
FE technology
Newmark with no algorithmic damping was used with consistent mass matrix
Central differences with diagonal mass matrix
timestep = 0.1 microsec = 1e-7 s
total number of time steps = 3500 (Newmark); 5000 (central differences)
this corresponds to total time = 350 microsec (Newmark); 500 microsec (central differences)
Wavefront timetable
FE strain distribution in location 1
Raw comparison – no tricks
-4
2
-4
LOC0 - exp vs. FE axisym, CD
x 10
2
1
mesh1, medium striker for all figures
LOC1 - exp vs. FE 3D, NM
x 10
1
axial strain
urmas data march 2007
0
0
-1
-1
-2
-2
-3
-3
experiment
FE analysis
-4
0
1
2
3
4
experiment
FE analysis
5
6
-4
0
0.5
1
1.5
2
2.5
3
-5
-4
x 10
-4
axial strain
2
x 10
-4
LOC2 - exp vs. FE 3D, NM
x 10
2
1
1
0
0
-1
-1
-2
-2
-3
LOC3 - exp vs. FE 3D, NM
x 10
-3
experiment
FE analysis
-4
3.5
0
0.5
1
1.5
2
time
2.5
3
experiment
FE analysis
3.5
-4
x 10
-4
0
0.5
1
1.5
2
time
2.5
3
3.5
-4
x 10
Are FE or experimental results closer to reality? Where is the truth?
Where is the truth?
Are experimental or FE results closer to reality?
•
When trying to reveal the ‘true’ behavior of a mechanical system we are using the
experiment
•
When trying to predict the ‘true’ behavior of a mechanical system we are accepting a
certain model of it and then solve it analytically and/or numerically
•
Physical laws (models) as
– Newton’s
– Energy conservation
– Theory of relativity
cannot be proved (in mathematical sense)
•
Often we say that it is the experiment which ultimately confirms the model
•
But experiments, as well as the numerical treatment of models describing the nature,
have observational thresholds. Sometimes, the computational threshold of
computational analysis are narrower than those of experiment
•
Let’s ponder about limits of applicability of continuum model as well as about limits of
modern computational approaches when applied to approximate solution of
continuum mechanics
Continuum mechanics
• Deals with response of solid and/or fluid
medium to external influence
• By response we mean description of
motion, displacement, force, strain and
stress expressed as functions of time and
space
• By external influence we mean loading,
constraints, etc. Expressed as functions of
time and space
Solid continuum mechanics
• Macroscopic model disregarding the corpuscular
structure of matter
• Continuous distribution of matter is assumed –
continuity hypothesis
• All considered material properties within the
observed infinitesimal element are identical with
those of a specimen of finite size
• Quantities describing the continuum behaviour
are expressed as piecewice continuous
functions of time and space
Governing equations
• Cauchy equations of motion
• Kinematic relations
 ijeng
t
t


uj
1   ui
2 0  0
 x
 xi
j





 tt ji
 tx j
 tfi  t txi
For transient problems, we are interested in, the
inertia forces are to be taken into account
t
t
t
t


u

u

u

uk 
j
Green_Lagrange
i
k
1
 ij
2 0  0  0
0
 x


x

x

x
j
i
i
j


• Constitutive relations

eng
ij
 Cijkl 
eng
kl
Sij  Dijkl klGreen _ Lagrange
In linear continuum mechanics
the equations are simplified to
 ui
1  ui u j 
 f i   2 ,  ij 

,  ij  Cijkl  kl
x j
t
2  x j xi 
 ji
2
• Inertia forces are considered
• loading is
• localized in space
• of short duration with a short rise time
The high frequency components are of utmost importance
Wave equation 2D - plane stress,
Lamé equations – expressed in displacements only
 2u
E
 2 
t
1  2
  2u
  2u  2 v 
 2v 
 2 
  G 2 
,
xy 
yx 
 x
 y
 2v
E
 2 
t
1  2
  2v
  2 v  2u 
 2u 
 2 
  G 2 
.
xy 
 y
 x xy 
Nondispersive and uncoupled solutions in unbound region
u  f x  c3t 
v  h y  c3t 
c3 
E
 1  2


Longitudinal
dilatational
irrotational
extension
P … primary
u  F  y  c2t 
v  H x  c2t 
c2 
G

Transversal
shear
rotational
distortion
equivolumetrical
S … shear
Velocities for 2D plane strain and 3D
  2G
c1 


1 
 1   1  2 
E
E

1   1  2 
P … primary
E
G
21   
Within the scope of linear theory of elasticity,
P and S waves are uncoupled.
Typical values for steel in m/s
For E = 2.1e11 Pa,  = 7800 kgm^-3,  = 0.3
c0  E /   5189  1D wave, slender bar
c1  (2G   ) /   6020  P wave for plane strain, 3D
c2  G /   3218  S wave, shear
c3  E /(  (1   2 ))  5439  P wave for plane stress
cR  0.9274c2  2984  R wave, Rayleigh for   0.3
Analytical solutions for bounded regions
are difficult and exist
for bodies with simple initial and
boundary conditions only
One way to solve the problem is to apply the
Fourier transform in space and the Laplace
transform in time on equations of motion.
The subsequent inverse procedure leads to
infinite series of improper integrals.
To indicate amount of effort to be exerted, the
formulas derived by F. Valeš are shown.
A typical expression
for a stress component
in a transiently loaded thin
elastic strip is shown
Relations between integrand
quantities and integral variable
Roots of these so called dispersive relations have to be
evaluated numerically before the integration itself.
Roots of dispersive relations
After all,
the analytical solution ends up
with a numerical evaluation
• Only finite number of terms can be summed up
• Integration itself has to be provided numerically
• This led to development of approximate
numerical solutions
Today, approximate methods of
solution prevail
 Discretization
 Finite difference method
 Transfer matrix method
 Matrix method
 Displacement formulation
 Force formulation
 Finite element method
 Displacement formulation
 Force formulation
 Hybrid finite element method
 Mixed finite element method
Numerical methods in FEA of
continuum
 Equilibrium problems
K (q) q  Q
solution of algebraic equations
 Steady-state vibration problems
K  Ω 2 M q  0
generalized eigenvalue problem
 Propagation problems
  F int  F ext , F int   B T σ dV
Mq
V
step by step integration in time
In linear cases we have
  Cq  Kq  P(t )
Mq
Unbound frequency response of
continuum
For fast transient problems as shock and impact the high
frequency components of solutions are of utmost
importance. In continuum, there is no upper limit of the
frequency range of the response. In this respect
continuum is able to deal with infinitely high frequencies.
This is a sort of singularity deeply embedded in the
continuum model.
As soon as we apply any of discrete methods for the
approximate treatment of transient tasks in continuum
mechanics, the value of upper cut-off frequency is to be
known in order to ‘safely’ describe the frequencies of
interest.
Dispersive properties
of 1D and 2D constant strain elements
When looking for the upper frequency limit
of a discrete approach to continuum problems,
we could proceed as follows
s

  5s
  cT
c  5000 m/s
f  1/ T
f  c /(5s )
• Characteristic element size
• Wavelength to be registered
• How many elements into the
wavelength
• Wavelength to period relation
• Wave velocity in steel
• Frequency to period relation
• The limit frequency
• For 1 mm element we get
5000
f 
 1106 Hz  1 MHz
5  0.001
Limits of continuum, FE analysis and experiment
characteristic sizes and corresponding frequencies
atom size
austenite steel grain size
1 mm finite element
1 MHz level
1 GHz level
FE analysis range from 0.1mm to 100mm
maximum exp. sampling limit 100 MHz - 14 bits
8
10
6
frequency in [Mhz]
10
Where is the continuum limit?
4
10
2
10
0
10
-2
10
-10
10
-8
10
-6
10
size in [m]
-4
10
-2
10
Validity limits of a model
• Model is a purposefully simplified concept of
a studied phenomenon invented with the
intention to predict – what would happen if ..
• Accepted assumptions (simplifications)
specify the validity limits of the model
• Model is neither true nor false
• Regardless of being simple or complicated, it
is good, if it is approved by an experiment
Using a model outside its limits
is a blunder
• Using a model outside of its validity limits
leads to erroneous results and conclusions
• This is not, however, the fault of the
model, but pure consequence of a poor
judgment of the model’s user
• Model gives no warning. Lot of checks
might be satisfied and still …
Blunders are easy to commit
• Point force is
– prohibited in continuum,
– frequently used in FE analysis
• Employing smaller and smaller elements, leads
to singularity, since we are coming closer and
closer to ‘continuum’ revealing thus
unacceptable behavior of the point force in
continuum
• An example follow
Example of a transient problem
Classical Lamb’s problem
1m
radial
L or Q
1m
B
Elements
L, Q, full int.
Consistent mass
axisymmetric
Mesh
Coarse 20x20
Medium 40x40
Fine
80x80
Newmark
Loading
a point force
equiv. pressure
axial
C
A
p0
t
Timp
Primary wave
R waves
S waves
P waves
Again, where is the
first nonzero P-wave
appearance?
Lamb, presure loading, rectangular pulse, velocity distribution, FEA
Pollution-free energy production
How to avoid blunders
By knowing and understanding
– the assumptions
– instrumental limitations
By providing
– validity checks
• within the model itself
• comparing with other methods
FE self-check considerations
Numerický experiment může mez rozlišení simulovat a napomoci tak experimentu.
Dá se ukázat, jak zjištěná rychlost šíření rozruchu závisí na hodnotě meze pozorování
V pozorovaném místě
známe
- vzdálenost od
aplikace budící síly,
- čas, kdy dorazí
„nenulový“ signál
můžeme tedy
vypočítat rychlost
šíření rozruchu.
Experimental
threshold
Numerical analysis should be robust
It should inform us about its limits
Should be independent of
–
–
–
–
mass matrix formulation
the method of integration
meshsize
element type
It is not always so, very often our results are
method dependent. See the next example.
Validity self-assessments, NM vs. CD
Computational infinite speed of wave propagation can be explained by analyzing
  Kq  Pt 
the time marching algorithms for Mq
Explicit (central differences)
Implicit (Newmark)
Equilibrium is considered
at time t
t  t
and leads to repeated solutions of
1
~
Mq

Pt
t  t
2
t
ˆq
ˆ
K
t  t  Pt  t
generally K , M banded; K 1 , M 1  full; M 1 is diagonal, only if M is diagonal
ˆ can never be diagonal computatio nal speed of propagatio n is not bounded
K
2
2
~
Pt  Pt   K  2 M q t  2 Mqt  t
t
t


ˆ K 1 M
K
 t 2
ˆ


P
t  t  Pt  t  M c1q t  c 2 q t  c 3 q t 
Instrumental limitation
Floating point representation of real numbers
threshold
Memory size limitation
Meshsize and time step limitations
machine epsilon
FFT frequency analysis
Now, let’s concentrate on the frequency analysis frequency
analysis of the loading pulse and of axial and radial
displacements obtained in the outer corner node of
location C by means of NM and CD operators for the
mesh1. The normalized power spectra are plotted in the
range from 0 to Nyquist frequency together with the
power spectrum of the loading pulse.
mesh1
mesh2
mesh3
mesh4
timestep [s]
1e-7
1e-7/2
1e-7/4
1e-7/8
sampling rate [MHz]
10
20
40
80
Nyquist frequency [MHz]
5
10
20
40
Assessment by
frequency analysis
2.5
x 10
imput pulse
4
x 10
mesh1, disp, corner node
-5
radial NM
axial NM
radial CD
axial CD
3
2
1.5
2
1
1
0.5
0
0
0
6
x 10
1
2
time
3
4
0
1
2
time
-5
x 10
transfer function, input vs. rad. disp
-10
3
4
x 10
-5
NM
CD
breathing frequency
zig-zag L frequency
zig-zag S frequency
FE limit frequency cons
FE limit frequency diag
5
4
3
2
1
0
0
0.5
1
1.5
2
2.5
3
frequency from 0 to Nyquist
3.5
4
4.5
5
x 10
6
transfer function, rad. disp, medium striker
Validity self-assessments
Mesh- and timestep refinement
x 10
-10
axisym, mesh1, corner node
1
NM
CD
breathing frequency
zig-zag L frequency
zig-zag S frequency
3.5
3
2.5
-5
axisym, mesh2, corner node
NM
CD
breathing frequency
zig-zag L frequency
zig-zag S frequency
0.8
0.6
2
0.4
1.5
1
0.2
0.5
0
0
0.5
1
1.5
2
x 10
transfer function, rad. disp, medium striker
x 10
1
x 10
0.8
0.6
-5
0
0
axisym, mesh3, corner node
1
NM
CD
breathing frequency
zig-zag L frequency
zig-zag S frequency
0.6
0.2
0.2
2
x 10
x 10
0.8
0.4
0.5
1
1.5
frequency range from 0 to 2 MHz
1
1.5
6
2
x 10
0.4
0
0
0.5
6
0
0
-5
6
axisym, mesh4, corner node
NM
CD
breathing frequency
zig-zag L frequency
zig-zag S frequency
0.5
1
1.5
frequency range from 0 to 2 MHz
2
x 10
6
Synergy of experiment and FE
analysis
FE analysis needs input data for
computational models from experiment
Experimental analysis could benefit from FE
in ‘proper’ settings of observational
thresholds
Four pulses, their mean values, power spectra
and input energies
Mean value, c0 value, norm and square root of the
input energy for different pulses
Terms and statistical tools for
quality assessment of close solutions
• Standard deviation (směrodatná odchylka) of
a sample vector (also called series of
measurements, signal, list of samples,
variables) having n elements is a scalar
quantity defined by
• Variance (rozptyl) is the square of the
standard deviation.
 1 n
2
xi  x  
y

 n  1 i 1

1
2
1 n
x   xi
n i 1
Covariance
• is the measure of how much two variables
vary together (as distinct from variance,
which measures how much a single
variable varies).
• The algorithm for covariance is
[n,p] = size(X);
X = X – ones(n,1)*mean(X);
C = X’*X/(n-1);
Correlation
• The correlation, also called correlation
coefficient, indicates the strength and
direction of a linear relationship between
two (or more) variables. The correlation
refers to the departure of two (or more)
variables from independence.
Respective means
n
rxy 
 ( xi  x )( yi  y )
i 1
(n  1) s x s y
and
Corresponding standard deviations
Geometric interpretation of
correlation
• The correlation coefficient can also be
viewed as the cosine of the angle between
the two vectors of samples drawn from the
two random variables.
xy
cos uncentered 
x y
Validity assessments, Axisymmetric and 3D elements
2.5
x 10
-6
Axisym vs. 3D elements, medium striker, CD
6
-5
x 10 Axisym vs. 3D elements, medium striker, CD
axi
3D
axi
3D
2
axial displacements, location 1, layer 1, node 1 [m]
radial displacements, location 1, layer 1, node 1[m]
5
1.5
1
0.5
0
4
3
2
1
-0.5
-1
0
500
1000
steps
1500
0
0
500
1000
steps
1500
Quantitative measure of
agreement quality
The agreement of solutions obtained by a different element
types is excellent – for a given loading and the employed time
and space discretizations, there is almost no ‘measurable’
difference.
One way to measure the difference between two solutions,
having the form of a vector in n-dimensional space (n is the
number of time steps in this case), is to compute the angle
between them. In our case the quality of agreement could be
quantified by
xy
cos  
 0.997847497832020
x y
This value is called the uncentered correlation coefficient.
radial displacements, corner node, locC
Assessment of solutions obtained by means
of NM and CD operators
and for different time and space discretizations
x 10
four meshes, original data, axisym, NM
-6
2
1
0
-1
m1
m2
m3
m4
-2
0
0.5
1
1.5
2
2.5
3
3.5
4
radial displacements, corner node, locC
time
x 10
4.5
x 10
-5
four meshes, original data, axisym, CD
-6
2
1
0
-1
-2
0
m1
m2
m3
m4
0.5
1
1.5
2
2.5
time
3
3.5
4
4.5
x 10
-5
Measure of sameness by correlation
relative directional error of m<i> to m<j> data, NM
0.025
0.02
m1, <i> index
m2, <i> index
m3, <i> index
0.015
0.01
0.005
0
2
3
m2 to m4, <j> index
4
relative directional error of m<i> to m<j> data, CD
0.025
0.02
m1, <i> index
m2, <i> index
m3, <i> index
0.015
0.01
0.005
0
2
3
m2 to m4, <j> index
4
Measure of sameness
by variance and covariance
variance and covariance for m1 to m4
-12
1.6
x 10
variance for NM
variance for CD
covariance NM. vs. CD
1.58
1.56
1.54
1.52
1.5
1.48
1.46
1
2
3
meshes m1 to m4
4
To reveal the true nature of reality
good-agreement solutions have to analyzed.
How?
• Avoid crude mistakes and omissions
• Avoid blunders – poor judgment
• The ‘sameness’ of close solutions should be assessed by statistical
tools
• Employed models should be robust, which means they should be
able give warning of their misuse. This is a pretty demanding task.
Presently, the robustness can only be achieved indirectly by means
of a posteriori checks by
– Comparing results obtained by different variant of the employed model
(coarse and fine meshes, different types of elements, coarse and fine
timesteps, different integration methods, etc)
– Analyzing the frequency contents of signals by Fourier analysis
• Results obtained by models should be checked by other models
and/or by experiments
• Since the experiment is just another tool for revealing the true nature
of reality, its results have to viewed by the prism of its validity limits
Conclusions
When trying to ascertain the reliability of modelling
approaches and the extent of their validity one
has to realize that the models as a rule do not
have self-correction features. That’s why we
have to let the models to check themselves, be
checked by independent models and let the
systematic doubt be our everyday companion.
Leftovers