Practical Aspects of the Finite Element Method Manuel Pastor Pablo Mira, José Antonio Fernández Merodo Centro de Estudios y Experimentación de Obras Públicas CEDEX, Ministerio de Fomento, SPAIN E.T.S. Ingenieros de Caminos Universidad Politécnica de Madrid, SPAIN
Slide 1
1.- Introduction (I) • FEM has become a powerful engineering analysis tool • Many engineers use commercial codes : ANSYS, ABAQUS,...... • Beginners face questions such as : i) Which element ? ii) How many elements and how big ? iii) Which material model ? iv) Is there any saving in reduced integration • This chapter adresses some of these questions It does it not in an exhaustive way Besides there are many other questions It is just a pocket guide for travellers Slide 2
2
1.- Introduction (II) • The questions we will address are: 1) The misteries of bending 2) Volumetric locking 3) The risks of reduced integration 4) Failure loads 5) Why cannot we choose all elements for geotechnical analysis ? 6) Our favorite solver : do we have any?
Slide 3
3
2.- The misteries of bending (I) • Not all elements perform satisfactorily under bending contitions • Linear triangles and bilinear quadrilaterals produce a much stiffer response than they should • This extra stiffness produces an increase in natural frequencies • Let us consider a square x∈[-1,+1], y∈[-1,+1]
M
M
Slide 4
4
2.- The misteries of bending (II) Continuum mechanics solution • Plane stress problem • Displacement field : M u − EI xy
M v − EI 1 − x 2
• Horizontal displacements at nodes :
u 0
M EI
−u 0 y M y y 0 xy 0 • Strain field x − EI
0 0
Slide 5
5
2.- The misteries of bending (III) Finite Element solution with linear quadrilateral • Nodal displacements will be : û1
−u 0 0
û2
u0
û3
0
• Element displacement vector : ûT
−u 0 0
û4
u0 0
ûT1 ûT2 ûT3 ûT4 y 4
3 x
• Shape functions: N1 N3
1 4 1 4
1 − x − y xy 1 x y xy
N2
1 4
N4
1 x − y − xy
1 4
1 − x y − xy
Slide 6
1
2
6
2.- The misteries of bending (IV) • B matrix : −1 y
0
1−y
0
1y
0
−1 − y
0
0
−1 x
0
−1 − x
0
1x
0
1−x
B
−1 x −1 y −1 − x
1 − y 1 x 1 y 1 − x −1 − y
• Strain field at element level is therefore: −u 0 y Bû
0 −u 0 x
• Spurious shear strain Slide 7
7
2.- The misteries of bending (V) • Computing external and internal work W ext 2M 2Mu 0 W int T d T Dd 2 u 0 yEu 0 yd 2 u 0 xGu 0xd 2u 0 EIu 0 2u 0 GIu 0 2u 0 E GIu 0
and equating them we get
u0
M EGI
• Shear locking
Slide 8
8
2.- The misteries of bending (VI) • With linear triangle T3 spurious shear strain also appears
−u 0
u0
1
0
2
0 −u 0
u0
2 1
• Constant strain in each element
• With quadratic elements T6 and Q8 spurious shear strain does not appear Slide 9
9
2.- The misteries of bending (VII) • Poor bending performance for T3 and Q4 • Since they are simple and easy to program attempts have been made to improve them • Wilson&Taylor, Simo&Rifai: Q4+Incompatible modes Quadratic enrichment. Distortion • Pian-Sumihara: Additional stress field • Battoz & Dhat: Beam test problem • Test performed with T3,T6,Q4,Q8,Q4WT and Q4PS • The most satisfactory Q4WT , Q4PS ,Q8,T6
Slide 10
10
Volumetric locking (I) • Let us assume a material with a Poisson coeficient close to 0.5 E • Bulk modulus: K 31−2 →∞ F
Y
1 X
3
4
1
2
6 5
• With υ=0.5 and linear triangles none of the nodes move • Two problems: critical parameter (υ=0.5 ) and element limitations (T3) • Critical parameter causes ill conditioning of the coeficient matrix Slide 11
11
Volumetric locking (II) • Nearly all classical displacement formulations present this type of limitations under incompressibility, in different degrees
• In 2D continuum mechanics problems, in incompressibility conditions we have at each point 1) 2 degrees of freedom ndof=2 2) 1 restriction (εv=0)
→ ndof/nrestr = 2
nrestr=1
• It would be desirable that in FE problems such ratio be maintained
Slide 12
12
Volumetric locking (III) • A simple test : incompressible patch test (Nagtegaal,Hughes) ndof = 2 nrestr ?
u a 1 a 2 x a 3 y a 4 xy v b 1 b 2 x b 3 y b 4 xy
y b3 b4 x
x a2 a4 y
v a 2 b 3 b 4 x a 4 y
• Element incompressibility restriction 1 1 x 1 y1
a2
0
1 1 x 2 y2
b3
0
1 1 x 2 y2
b4
1 1 x 2 y2
a4
0 0
Slide 13
Coefficient matrix is rank 3
r
ndof nrestr
2 3 13
Volumetric locking (IV) • Ratio values for different classical element types with full integration, in plane strain and axisymmetric conditions from Sloan & Randolph (1982)
T3 T6 T10 T15 Q4 Q8 Q12 Q17 PlaneStrain
1
4 3
3 2
8 5
2 3
1
1
8 7
Axisymmetric
1 3
2 3
9 10
16 15
2 5
2 3
10 13
16 19
Slide 14
14
Risks of reduced integration (I) • Reduced integration consists of using an integration rule of smaller degree than required to integrate exactly polinomials existing in the stiffness matrix of undistorted elements • It has two advantages: a) Reducing computational cost b) Improving incompressibility performance (Volumetric locking) • Two popular reduced integration rules are: a) One point for bilinear quadrilaterals b) Two by two for 8 noded quadrilaterals • It is interesting to study the problem from the point of view of deformation modes (eigenvectors of stiffness matrix) • Let us consider the case of the bilinear quadrilateral Slide 15
15
Risks of reduced integration (II) • Deformations modes of bilinear quadrilateral
Tv
Th
Rot
Dil
B1
B2
S1
• Description • Rigid body modes : Zero energy
Slide 16
S2
16
Risks of reduced integration (III) • Reduced integration: 1 point integration rule • Stiffness matrix which in general is : K B T . D. B d Is computed as : K B T0 . D. B 0 . W 0 • B is 8x3 and D is 3x3 → K is rank 3 K is 8x8 and rank 3 → 5 zero energy modes • 5 zero energy modes = 3 rigid body modes +... 2 bending modes Reduced integration of 8 noded quadrilaterals also causes Additional zero energy modes
Slide 17
17
Risks of reduced integration (IV) • Non rigid body zero energy modes might be a problem • Spurious bending modes • Ill conditioning → Hourglassing • Hourglass control • Selective integration
Slide 18
18
Failure loads (I) • Elastoplastic problems also present quasi incompressible conditions • Footing on vertical slope Material properties Tipo de material E(Pa) Suelo
Von Mises
Zapata Elásticolineal
y (Pa)
1.0E5 0.35 200.0 1.0E8 0.35
Perfect plasticity : H = 0.00 Softening : H= -0.01E
- Theoretical Solución Using límit analysis: - Failure mechanism direction 45º - P = 1154.7 N Slide 19
19
Failure loads (II) • How does mesh orientation affect the finite element solution of a failure mechanism problem ?
Slide 20
20
Failure loads (III) t3st
t6st, t6bb, t7st, t7bb
t4st, t4bb, t4sr
“Right” Orientation (r)
“wrong” Orientation (w)
Slide 21
21
Failure loads (IV) • Force-displacement diagrams for H = 0.00 Quadrilaterals
Triangles
Slide 22
22
Failure loads (V) • Failure mechanisms for H=0.00 with quadrilaterals and 3n triangles
Slide 23
23
Failure loads (VI) • Failure mechanisms for H=0.00 with 6n and 7n triangles
Slide 24
24
5.- Why cannot we use all elements ..........? (I) • Let us consider the saturated consolidation problem σy
σ 'y
τ xy
τ xy
σx
σx
p
τ ' xy
τ ' xy
= σ 'x
σ 'x
τ ' xy
τ xy τ xy
τ ' xy
σy
p
p
+
p
σ 'y
• Space discretization
′
S − mp b 0 m Su̇ − ∇ k∇p T
T
ṗ
∇ k w b 0 T
′
B d − Qp̄ f u T
.
.
Q u Hp̄ C p̄ f • Discretizing in time and solving the non linear problem with NR KT Δt
Q∗
T
.
−QΔt
dΔ ū n
.
−Q T Δt −Δt HΔt C Slide 25
dΔp̄ n
−
Gu −Δt G p
25
5.- Why cannot we use all elements ..........? (II) • Undrained incompressible limit kw→0; Q*→∝, jacobian will be : KT
−Q
−QT
0
.
dΔ ū n
.
d Δp̄ n
−
Gu −Gp
• Doesnt look good, determinant might be zero
• A complete mathematical treatment of this problem can be found in Babuska(1973) and Brezzi(1973) • A simpler a approach to the problem, although less complete can be found in Zienkiewicz & Taylor patch test (1986)
Slide 26
26
5.- Why cannot we use all elements ..........? (III) 1 m.
• Saturated soil column subject to surface load • Only uz allowed • All boundaries impermeable except top horizontal boundary kw
10 −7 m/s
n
0.333
E
7.492 10 8 (Pa)
0.2
s
2. 0 10 3 N/m 3
w
1. 0 10 3 N/m 3
q = 100 exp(-iwt) tz = 0 ux= 0
30 m
uz = 0
Slide 27
Z
∂ pw =0 ∂x
O
X
ux= 0
∂ pw =0 ∂z 27
5.- Why cannot we use all elements ..........? (IV) • Problem was discretized with a column of 20 elements
Pore-pressure distribution 1
1
0.8
0.8
0.6
0.6
z/L
z/L
q4p4
0.4
0.4
(a)
• Results a) Q* = 104MPa
(b)
0.2
0.2
0. 0.
•b) Q* =
Pore-pressure distribution
0.5
109MPa
1
0. 0.
1.5
0.5
2
Pore-pressure distribution
1
1
0.8
0.8 Q8P4
Q8P4
0.6
z/L
z/L
0.6
0.4
Q* = Water and grain compressibility
1.5
p/q
Pore-pressure distribution
q8p8
1
p/q
0.4
(c)
(d)
0.2
0.
0.2
0.
0.5
1
p/q
Slide 28
1.5
0.
0.
0.5
1
1.5
p/q
28
5.- Why cannot we use all elements ..........? (V) • Zienkiewicz & Taylor patch test (1986) KT
−Q
−QT
0
dΔ ū n
−
.
d Δp̄ n
.
.
}
Q K. QdΔp̄ n − G p Q T K −1 G u n p × n p = (n p × nu )(nu × nu )(nu × n p ) T
−1
−Gp
−1
dΔ ū n K Q. dΔp̄ n K G u −1
Gu
u
p
u
p
nu ≥ np
• This is a necessary but not sufficient condition Singularity would have to be tested in all cases Slide 29
29
5.- Why cannot we use all elements ..........? (VI) • Single element patch q4p4 • Constraints
u
p
• nu = 0 < np = 3 Test not satisfied Precribed displacement Prescribed Pressure
• Single element patch t6p3 • nu = 0 < np = 2 Test not satisfied
u
Slide 30
p
30
5.- Why cannot we use all elements ..........? (VII) • Six element patch T6P3 • Constraints • nu = 14 > np = 6 Test satisfied u
p
• Four element patch Q4P4
• nu = 2 < np = 8 Test not satisfied
• Four element patch Q8P4
• nu = 10 > np = 8 Test satisfied
• Conclusions Slide 31
31
5.- Why cannot we use all elements ..........? (VIII)
KT
−Q
−Q
0
T
dΔ ū n
.
d Δp̄ n
−
Gu −Gp
• If somehow the 0 is avoided there is no problem • Stabilization methods - Divergence, Fractional step; methods based in ideas used in FEM for fluid mechanics -Special implementation of Simo-Rifai
Slide 32
32
6.- Our favorite solver (Do we have any?) • One of the main computational tasks in the FEM is linear Equation system solving
• The issue is especially important with transient, non linear, ....multifield,.....3D analysis
• Which equation solver to use is therefore a crucial question to answer when thinking of building your own FE code
Slide 33
33
6.- Our favorite solver (Do we have any?) • Iterative solvers like preconditioned conjugate gradient are excellent choices : Easy to program Does not require renumbering • Iterative methods are specially good for really large problems: Tens or hundreds of thousands degrees of freedom • For not so large problems (thousands of dregrees of freedom) direct methods with special storage and renumbering schemes are also good choices • Non symmetric systems : skyline → some changes conjugate gradient → GMRES Slide 34
34
6.1 Conjugate gradient methods (I) • Based on the idea that the solution of Ax = b minimizes total potential 12 x T Ax − x T b • It is straightforward to prove the previous idea • The procedure would therefore be: Given xk and πk Find xk+1 such that πk+1 < πk • Not only we want πk+1 < πk but πk+niter << πk With niter not too large → we want fast convergence • In the kth iteration we use a set of linearly independent vectors p1, p2,....., pk and calculate the minimum potential in the space spanned by these vectors, thus obtaining xk+1 and also pk+1 in the process Slide 35
35
6.1 Conjugate gradient methods (II) • The algorithm can be summarized as follows: • Choose the starting iteration vector x¹ • Calculate the residual r¹=b-Ax¹ if r¹=0 quit else set p¹=r¹ do for k=1,2,..., until ‖rk ‖ ≤
T
rk rk
k
T
p k Ap k
x k1 x k k pk rk1 rk − k Apk
k
T
r k1 r k1 T
r k Ar k k1
pk1 p
enddo
k pk
Slide 36
36
6.1 Conjugate gradient methods (III) • Two important orthogonality properties pTi Apkj 0
PTj rk1 0
where
P j p1 , . . . , pj
• Convergence in at least n iterations. It should be faster • Rate of convergence depends on condition number of A Cond(A)=λn/ λ1 where λn is the largest eigenvalue of A and λ1 is the smallest eigenvalue of A • Large condition number → Determinant closer to 0 Small condition number → Closer to identity matrix Slide 37
37
6.1 Conjugate gradient methods (IV) • To increase rate of convergence → decrease condition number −1
−1
• Instead of solvingAx b, solve A Ax A b Where A is the preconditioner, easy to invert and as close as possible to A-1 A
−1
= diag(A) → Jacobi conjugate gradient Good for diagonal dominant A matrices −1
• Algorithm with preconditioner introduces: z A rk k
k
k
T
z k rk T
p k Ap k T
z k1 r k1 T
z k rk
pk1 zk1 k pk Slide 38
38
6.2 Skyline storage scheme (I) • Using a direct solution scheme and storing the hole nxn coeficient matrix is not an option considered in FE analysis • To do so would imply : a) A huge waste of storage space given the sparsity of The coefficient matrix (many coefficients are 0) b) A huge waste of computational effort since matrix operations on null coefficients are trivial • Because of this it is important to use an efficient storage scheme • The skyline storage scheme is a good alternative
Slide 39
39
6.2 Skyline storage scheme (II) • It consists of storing the stiffness matrix in a vector including the diagonal terms and the off diagonal terms between the diagonal And the farthest non zero off diagonal term, for each row (column) • Additionally it will be necessary to use an N component integer Vector to store the position in the stiffness vector of each diagonal component 2 −2 −2
3 −2
−1
0 −1 0
0
5 −3
0
0 −3 10
4
0 −2 0
0
0
0
i
1 2 3 4 5 6
7
8 9 10 11 12
ai 2 −2 3 −2 5 −3 10 −1 0 0 i
4 10
1 2 3 4 5
jdiagi 1 3 5 7 12
4 10
Slide 40
40
7.- Concluding remarks • We have addressed a few practical aspects of FE analysis • Bending behaviour: Standard low order displacement formulations Quadratic formulations Improved low order formulations : mixed..... • Volumetric locking : Standard low order formulations are bad Higher order formulations are better but ...... • Reduced integration : It is possible solution for volumetric locking But careful with it : Hourglassing Hourglass control or selective integration • Coupled formulations : Standard u-p formulations → nu > np u8p4 Stabilized : fractional step, Simo-Rifai, div • Solvers : Iterative: easy to program, very good for large problems 41 Direct with special storage and renumbering, also good Slide 41
