Project Presentation
« Element-Free Galerkin Method »
Alternative numerical methods in continuum mechanics by L. Noels
Written by:
Majérus Julien
01/03/2010
Verheylewegen Guillaume 1
Outline
Introduction
Meshless
Weak form equations
Integration Scheme
Examples
Meshless in fracture mechanics
Examples
Summary
Bibliography
2nd Master in Aerospace Ingineer 2009-
01/03/2010 2
2010
Introduction
In the last decades, new computational methods were
developed in order to solve effectively problems with
discontinuities, moving boundaries, large deformations. One
of the approaches was:
Element-Free Method
Meshless Method
Mesh-Free Method
No mesh support
2nd Master in Aerospace Ingineer 2009-
01/03/2010 3
2010
Introduction - History
1977 - Smooth particle hydrodynamics (SPH) (Lucy)
1992 - Diffuse element method (DEM) (Nayroles et al.)
1994 - Element free Galerkin method (EFG) (Belytschko et al.; Organ et al.)
1995 - Reproducing kernel particle method (RKPM) (Liu et al.)
1996 - Hp-clouds method (Duarte and Oden, 1996)
1997 - Partition of unity method (PUM) (Babuska and Melenk)
1998 - Local boundary integral equation method (LBIE) (Zhu, Zhang, and Atluri)
1998 - Meshless local Petrov±Galerkin method (MLPG) (Atluri and Zhu).
2004 - Meshless Natural Element method (NEM) (Yvonnet).
…
2nd Master in Aerospace Ingineer 2009-
01/03/2010 4
2010
Introduction - Utility of the meshless method ?
The famous discretized method, Finite Element Method (FEM), which has been
widely used in many engineering problem simulations, exhibits some
limitations:
Interpolation fails when elements become too distorted
FE re-meshing is time-consuming, but necessary in some cases :
New surface creation : fracture mechanics
Large element distortion : finite strains
Shear bands orientation are mesh dependent
2nd Master in Aerospace Ingineer 2009-
01/03/2010 5
2010
Meshless method - Formulation
Meshless and FE methods are both based on a discretized
displacement field:
Meshless consists in subdividing the structure in nodes. At each node,
we associate a weight function and a shape function defined on a small
domain. FEM uses shape functions which are defined on what we call
elements.
FE method - Meshless method
(Duflot PhD 2005)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 6
2010
Approximation definition (MLS)
Moving Least Squares (MLS) Field approximation:
Where pT(x) is the approximation base and a(x) the
coefficients.
The MLS consists in minimizing the weighted
norm:
MLS approximation
This last equation can be written as: A(x)a(x)=b(x)
where
2nd Master in Aerospace Ingineer 2009-
01/03/2010 7
2010
Approximation definition (MLS) (2)
Based on the MLS approximation. Field u(x) which is defined on all domain, can
be written as:
Where Φi(x) is denoted as the shape function associated to the node i
2nd Master in Aerospace Ingineer 2009-
01/03/2010 8
2010
Approximation definition (MLS) (3)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 9
2010
Interpolation definition (FEM)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 10
2010
Shape functions definition
The shape functions Φi(x) are independent of any mesh by
definition. Shape function depends on:
The approximation base
A weight function p(x)
The domain of influence wi(x)
Note that these parameters can be chosen arbitrarily and
even differ between different nodes in the domain.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 11
2010
Shape functions properties Φi(x)
Zero-value outside the domain of the weight function.
Non-Zero value on the whole domain of influence, if A(x) is not singular for
any point. (There is a minimum size for the domain of the weight function)
and thus .The shape functions don't interpolate the
displacement at nodes. This is a major difference with FEM where one gets
interpolation at nodes. Note that it exists different kinds of approximations, or in
some cases interpolation as the C-NEM method which uses "natural element"
based on Voronoi diagram)
Partition of unity: Rigid body modes can be computed.
Consistency : If the exact solution belongs to the base of the shape functions,
then the solution is exactly fitted by the MLS approximation.
Continuity: where and
2nd Master in Aerospace Ingineer 2009-
01/03/2010 12
2010
Base functions definition p(x)
Since base function are , shape functions have the same degree of
continuity as the weight function.
The choice of the base function will define the size of the matrix A. Since this
one has to be inverted at each node, it is time-consuming. In addition to this, the
node domain has to be bigger with a higher order of the base if we want to avoid
A to be singular.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 13
2010
Weight function definition wi(x)
The weight function is a function which decrease with the distance from
node. Outside the node domain it must be equal to zero.
Circular domain (2D)
Rectangular domain (2D)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 14
2010
Choice of Weight function (2)
If we define , in the literature we find functions:
Troncated Gaussian
Not !
Modified Gaussian
Spline of order 4
2nd Master in Aerospace Ingineer 2009-
01/03/2010 15
2010
Weak form of elasticity equations (FEM)
If we consider an arbitrary test function , we
can write the weak form of the static equilibrium
equations:
By part integration and Gauss theorem:
The Galerkin procedure assume that same
function are used :
2nd Master in Aerospace Ingineer 2009-
01/03/2010 16
2010
Weak form of elasticity equations (FEM) (2)
With the elastic material behavior D:
As is arbitrary, the last equation can be written on the form:
2nd Master in Aerospace Ingineer 2009-
01/03/2010 17
2010
Weak form of elasticity equation (Meshless)
FEM equation are still valid for meshless.
However, due to the not interpolating shape functions, the
Essential Boundary Conditions (EBC) can not be enforced
the same way. One can not impose a nodal displacement
as in FEM.
As a consequence, modification have to be taken into
account to be able to satisfy EBC. Three main methods
(a,b,c) have been developed.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 18
2010
(a) Weak form with Lagrangian multiplicator
In this weak form, Lagrangian multiplicators are used to impose the
displacement boundary conditions on
With the Galerkin procedure and the arbitrary value of , we can write in matrix
form
is an interpolating shape
function with the curvilinear variable
s along
S is a diagonal matrix whose
element ii is equal to 1 if ui is
imposed and 0 otherwise. (useful
if, e.g., ux is imposed but uy is free)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 19
2010
(b) Weak form without Langrangian multiplicator
This method is similar to the previous one. Lagrange multiplier have a physical
meaning: they represent reaction forces at the boundary. These reaction forces t
can be directly identified in the weak form:
Where
Thus the system is simply: K q = g
The system of equations is smaller than in the previous method
2nd Master in Aerospace Ingineer 2009-
01/03/2010 20
2010
(c) Weak form with penalty factor
Here we use a penalty factor in the equation:
The difficulty of this method is too chose the right β. Too high the system is badly
conditioned, and too low the imposed displacements are not accurate.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 21
2010
Integration Scheme
Integrals over the domain and have to be evaluated. In contrast with the
FEM, we have not any "element" which are a natural solution to integrate over the
domain. Three main solutions can be used.
Fig: Integration scheme, (a) "meshnode" , (b) "background mesh", (c) "node domain"
(Duflot PhD 2005)
(a) Use nodes distribution to form "elements" where a Gauss integration will be performed on each
element.
(b) Use an independent structured mesh over the domain where Gauss integration will be performed.
(c) Integrate over nodes "integration domain”
2nd Master in Aerospace Ingineer 2009-
01/03/2010 22
2010
Process
1 2
STIFFNESS:
LOOP1 (At each integration point over )
LOOP2 (at each node i)
END LOOP2 1
LOOP3 (at each node i)
END LOOP3
ADD to K
END LOOP1
LOOP4 (At each integration point over )
... 2
ADD to K
END LOOP4
EXTERNAL FORCES
SOLVE
GET DISPLACEMENT FIELD:
2nd Master in Aerospace Ingineer 2009-
01/03/2010 23
2010
Algorithm optimisation
The time consumed in the algorithm come from:
The number of integration points in each cell
The inversion of the matrix A. For a linear choice of base functions
the size is only 3x3, but it has to be inverted at each integration point.
The loop over each node (LOOP2,LOOP3). Since the weight function
is non zero only in a small domain around the related node, at each
integration point only few node around this point will get a weight
function different from zero. By consequence it's insane to loop over all
nodes. One way to improve this is to take into account the relation
between nodes by a quadtree formulation which will give an tree fast to
scan.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 24
2010
Miscellaneous
Truly meshless method:
Non-element interpolation technique
Non-element approach for integrating the weak form
Example
a truly meshless method = Meshless local Petrov-Galerkin method
(MLPG), no need of mesh or "integration mesh »
a meshless method = Element free Galerkin method (EFG), need
of "integration mesh".
2nd Master in Aerospace Ingineer 2009-
01/03/2010 25
2010
1D example : the bar
Analytical solution :
2nd Master in Aerospace Ingineer 2009-
01/03/2010 26
2010
2D example : bi-clamped plane strain body
2nd Master in Aerospace Ingineer 2009-
01/03/2010 27
2010
Fracture mechanics - reminder
A crack involves a singular displacement field close to the crack tip. Under the hypothesis
of Linear Elastic Fracture Mechanics (LEFM), it can be shown that the displacement is
asymptotically almost proportional to where r is the distance from the crack tip. For
example, for a crack opening in mode I only, the displacement is :
Note that if one can accurately compute the stress intensity factor KI, the displacement
field is (asymptotically) completely defined. Similar relationships can be shown for the two
other modes.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 28
2010
How to insert a crack in the model
A crack introduces a discontinuity in the displacement field. The model
should be able to reproduce this discontinuity. Two possibilities:
The crack belongs to the boundary of the domain. (ex: in FEM the
mesh has to be conform to the crack. In propagation, a new mesh is
necessary at each step)
The shape functions are discontinuous. As a consequence
will be discontinuous too. (e.g.: X-FEM & meshless)
In fact, the discontinuity in the shape function is introduced through a
discontinuity in the weight function. Remember that
2nd Master in Aerospace Ingineer 2009-
01/03/2010 29
2010
How to insert a crack in the model: a) visibility criterion
Belytschko et al. /Cornput. Methods Appl. Mech. Engrg. 139 (1996) 3-47
The crack acts as an opaque surface. If a point can't "see" a node, then the
weight (related to the node) of this point is set to zero, i.e.:
Easy to compute
But introduces a spurious discontinuity (along AB)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 30
2010
How to insert a crack in the model: b) Diffraction criterion
Idea : avoid the spurious discontinuity introduced by the visibility
criterion
Solution : Use the "diffracted" distance d' in order to compute the weight
w(d'). If the line between the point x and the node xi crosses the crack,
the diffracted distance is :
d’ = max(S0=d , S1+S2)
Belytschko et al. /Cornput. Methods Appl. Mech. Engrg. 139 (1996) 3-47
2nd Master in Aerospace Ingineer 2009-
01/03/2010 31
2010
How to deal with the singularity
The displacement field is singular close to the crack tip. The node density should
be increased in this region. Several methods exist, but the main idea is to
introduce the asymptotic solution in the model. Either we introduce new degrees
of freedom (as X-FEM does), or we can simply enrich de base :
The matrix A that has to be inverted to compute the shape functions is now a 7
by 7 matrix. Moreover, as shape functions are now able to represent a singular
field, a higher number of integration point is necessary.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 32
2010
Numerical example: 2D mixed-mode loading - Discretisation
Because the formulation is quite similar to the X-FEM one, results will be
compared to a problem solved with XFEM (the latter has been computed with
the FEM/X-FEM software Morfeo).
Meshless (Matlab) – XFEM (Morfeo)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 33
2010
Numerical example: 2D mixed-mode loading - Displacements
Imposed displacement (traction + shear)
Meshless (Matlab) – XFEM (Morfeo)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 34
2010
Numerical example: 2D mixed-mode loading - Von Mises stresses
Meshless (Matlab) – XFEM (Morfeo)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 35
2010
Stress Intensity Factors Computation
Stress intensity factors are extremely important because they completely characterise the
displacement and stress field near the crack tip. Moreover, they are necessary to predict
the crack propagation direction and the crack propagation speed in fatigue. The most
accurate method to compute them is the J integral. For a straight crack in the direction x:
In practice, J is evaluated over a domain, not a contour (divergence theorem). The
relationship between J and the SIF's is :
Mode I only ( ): directly determined with J
Mixed mode: More equations are necessary. It is possible to show that by superposing
(linear elasticity framework!) auxiliary virtual solutions to the numerical solution (use of
interaction integrals), one can extract the stress intensity factors. The demonstration of this
is beyond the scope of this presentation. However, we have applied this method :
2nd Master in Aerospace Ingineer 2009-
01/03/2010 36
2010
Stress Intensity Factors Computation
Meshless : KI = 0.17262 KII = 0.096375
XFEM : KI = 0.171621 KII = 0.0969399
Our model is not so bad!
2nd Master in Aerospace Ingineer 2009-
01/03/2010 37
2010
Summary
Only a set of node is needed
But a mesh is generally needed in order to compute a
Gauss quadrature
One can easily modify the base functions or the weight
functions where we want in the domain
Arbitrary discontinuities can be modelled without
modifying the set of node
Stress and strain computed with the meshless methods
are smooth (in FEM they are discontinuous between
elements)
2nd Master in Aerospace Ingineer 2009-
01/03/2010 38
2010
Other applications
2nd Master in Aerospace Ingineer 2009-
01/03/2010 39
2010
Bibliography
M. Duflot, Application des méthodes sans maillage en mécanique de la rupture, Thèse
de doctorat, Ulg, 2005
T. Belytschko, Y.Y. Lu, [Link], Element-free Galerkin methods, International Journal for
Numerical Methods in Engineering vol 37, 229-256, 1994
T. Belytschko, L. Gu, Y.Y. Lu, Fracture and crack growth by element free Galerkin
methods, Modelling and Simulation in Materials Science and Engineering vol 2, 519-534,
1994
T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an
overview and recent developments, Computer methods in applied mechanics and
engineering, vol 139, 3-47, 1996
J. Yvonnet, Nouvelles approches sans maillage basées sur la méthode des éléments
naturels pour la simulation numérique des procédés de mise en forme, Thèse de doctorat,
Ecole Nationale Supérieur d’Arts et Métiers Centre Paris, 2004.
2nd Master in Aerospace Ingineer 2009-
01/03/2010 40
2010