O metodě konečných prvkůLect_5.ppt

M. Okrouhlík

Ústav termomechaniky, AV ČR, PrahaPlzeň, 2010

Pitfalls of FE Computing

N u m e r i c a l m e t h o d s i n F E A E q u i l i b r i u m p r o b l e m s

QqqK )( s o l u t i o n o f a l g e b r a i c e q u a t i o n s

S t e a d y - s t a t e v i b r a t i o n p r o b l e m s 0qMK 2Ω

g e n e r a l i z e d e i g e n v a l u e p r o b l e m

P r o p a g a t i o n p r o b l e m s

VV

d, Tintextint σBFFFqM

s t e p b y s t e p i n t e g r a t i o n i n t i m e I n l i n e a r c a s e s w e h a v e

)( tPKqqCqM

Now a few examples

from• Statics• Steady-state vibration• Transient dynamics

using• Rod elements• Beam elements• Bilinear (L) and biquadratic (Q) plane elements

S t a t i c l o a d i n g o f a c a n t i l e v e r b e a m b y a v e r t i c a l f o r c e a c t i n g a t t h e f r e e e n d .

B e a m 4 - d o f e l e m e n t s ( E u l e r - B e r n o u l l i ) a r e u s e d .

1

2

3

4

2

22

3

2sym.

36

32

3636

2

l

l

lll

ll

l

EIk

PKq

EI

PLv

3

3

tip

1 1 . 9 0 4 7 6 1 9 0 4 7 6 1 9 0 5 e - 0 0 3 1 . 9 0 4 7 6 1 9 0 4 7 6 1 9 0 4 e - 0 0 3 1 0 1 . 9 0 4 7 6 1 9 0 4 7 6 0 4 3 e - 0 0 3 1 . 9 0 4 7 6 1 9 0 4 7 6 1 9 0 4 e - 0 0 3

1 2 3 k m a x

P

‘Exact’ formula for thin beams

A thin cantilever beam, vertical point load at the free end, four-node plane stress elements

FE technology

V

VdTEBBk

full integration 1 … anal, 2 … Gauss q., different types of underintegration

0 1 2 3 4 5 6 7 8 9 100

1

2

3

4

5

6x 10

-3 maximum displacements for kmax = 10

case

"EXACT"

KHGS 0 1 2 3 4 5 6

LDefault in ANSYSReduced integration in Marc

Exact integration Exact integrationdoes not yield‘exact’ results

Data obtained with10 elements

Gauss quadrature and underintegrated elements

Summary_2Isoparametric approach

Four-node bilinear element, plane stress, full integration

Four-node bilinear element, plane stress, full integration

Comparison of beam and L elements used for modelling of a static loading

of a cantilever beam

• Beam … even one element gives a negligible error

• L … too stiff in bending, tricks have to be employed to get “correct” results

1 2 3

4 5 6

7 8

0 5 100

1

2

3x 10

4GEIVP

consistent

1 2 3

4 5 6

7 8

0 5 100

5000

10000

15000GEIVP

diagonal

  A single four-node plane stress element Generalized eigenvalueproblem

Full integration

0qMK )(

1 2 3 kmax

0qMK )(

Free transverse vibration of a thin elastic cantilever beamFE vs. continuum approach (Bernoulli_Euler theory)

Generalized eigenvalue problem

Two cases will be studied- L (bilinear, 4-node, plane stress elements)- Beam elements

0 1 2 3 4 5 6 7 8 9 100

2000

4000

6000

8000

10000axial, bending and FE frequencies [Hz]

diagonal

1, 100.506 2, 586.4532 3, 1299.7804 4, 1501.8821

5, 2651.6171 6, 3862.638 7, 3931.881 8, 5258.5331

9, 6313.8768 10, 6567.2502 11, 7802.0473 12, 8571.0469

13, 8907.7934 14, 9804.6267 15, 10394.9981 16, 10520.9827

Natural frequencies and modes of a cantilever beam … diagonal mass m. Four-node plane stress elements, full integration

Sixth FE frequency is the second axialSeventh FE frequency is the fifth bending

Higher frequencies are uselessdue to discretization errors

Analytical axial Analytical bending

FE frequencies

This we cannot say without looking at eigenmodes

Frequencies in [Hz] versus counter

L

0 1 2 3 4 5 6 7 8 9 100

2000

4000

6000

8000

10000axial, bending and FE frequencies [Hz]

consistent

1, 101.2675 2, 615.4501 3, 1302.6992 4, 1663.1767

5, 3135.596 6, 3941.6886 7, 4996.7233 8, 6682.2631

9, 7215.1641 10, 9593.4882 11, 9739.1512 12, 12430.9913

13, 12740.3253 14, 14990.9545 15, 16164.1614 16, 16950.4893

Natural frequencies and modes of a cantilever beam … consistent mass m. Four-node plane stress elements, full integration

0 1 2 3 4 5 6 7 8 9 100

2000

4000

6000

8000

10000frekvence analyticke axialni (x), ohybove (o) a MKP frekvence (*)

1 1.5 2 2.5 3 3.5 4 4.5 50

5

10

15

20

25

konsistentni formulace

relativni chyba v % pro axialni (x) a ohybove (o) frekvence0 1 2 3 4 5 6 7 8 9 10

0

2000

4000

6000

8000

10000frekvence analyticke axialni (x), ohybove (o) a MKP frekvence (*)

1 1.5 2 2.5 3 3.5 4 4.5 5-20

-10

0

10

20

diagonalni formulace

relativni chyba v % pro axialni (x) a ohybove (o) frekvence

Relative errors for axial (x) and bending (o) freq.

Relative errors for axial (x) and bending (o) freq.

Eigenfrequencies of a cantilever beamFour-node bilinear element, plane strain

Relative errors [%] of FE frequencies

x … axial – continuum, o … bending – continuum, * … FE frequencies

1 2 3 kmax

1

2

3

4

2

2

22

4sym.

22156

3134

135422156

420

l

l

lll

ll

lAm

2

22

3

2sym.

36

32

3636

2

l

l

lll

ll

l

EIk

0qMK )(

Free transverse vibration of a thin elastic cantilever beam

FE (beam element) vs. continuum approach (Bernoulli_Euler theory)

Free transverse vibration of a thin elastic beamANALYTICAL APPROACH The equation of motion of a long thin beam considered as continuum undergoing transverse vibration is derived under Bernoulli-Euler assumptions,

namely-there is an axis, say x, of the beam that undergoes no extension,-the x-axis is located along the neutral axis of the beam, -cross sections perpendicular to the neutral axis remain planar during the

deformation – transverse shear deformation is neglected,-material is linearly elastic and homogeneous,-the y-axis, perpendicular to the x-axis, together with x-axis form a principal plane of

the beam.

These assumptions are acceptable for thin beams – the model ignores shear deformations of a beam element and rotary inertia forces.

For more details see Craig, R.R.: Structural Dynamics. John Wiley, New York, 1981 or Clough, R.W. and Penzien, J.: Dynamics of Structures, McGraw-Hill, New York,

1993.

T h e e q u a t io n i s u s u a l l y p r e s e n te d in th e f o r m

),(2

2

2

2

2

2

txpt

vA

x

vEI

x

w h e r e x i s a l o n g i tu d in a l c o o r d in a t e , v i s a t r a n s v e r s a l d i s p l a c e m e n t o f th e b e a m in y

d i r e c t io n , w h ic h i s p e r p e n d ic u la r t o x , t i s t im e , E i s t h e Y o u n g ’ s m o d u lu s , I i s t h e p l a n a r

m o m e n t o f in e r t i a o f t h e c r o s s s e c t io n , A i s t h e c r o s s s e c t io n a l a r e a a n d i s t h e d e n s i t y .

O n th e r i g h t h a n d s id e o f th e e q u a t io n th e r e i s t h e lo a d in g ),( txp - g e n e r a l l y a f u n c t io n

o f s p a c e a n d t im e - a c t in g in th e x y p l a n e .F o r f r e e t r a n s v e r s e v ib r a t io n s w e h a v e z e r o o n

th e r i g h t - h a n d s id e o f E q . ( 1 ) . I f t h e b e n d in g s t i f f n e s s E I i s i n d e p e n d e n t o f t im e a n d

s p a c e c o o r d in a te s w e c a n w r i t e

02

2

4

4

t

v

EI

A

x

v .

( 4 a ) A s s u m i n g t h e s t e a d y s t a t e v i b r a t i o n i n a h a r m o n i c f o r m

)cos()(),( txVtxv

w e g e t

0)(d

)(d 44

4

xVx

xV

( 4 b ) w h e r e w e h a v e i n t r o d u c e d a n a u x i l i a r y v a r i a b l e b y

)/(24 EIA ( 5 )

T h e gen era l so lu tio n o f E q . (4 ) c an b e assu m ed (se e K re ys ig , E .: A d v an ced E n g in e e rin g M ath em atic s , Jo h n W ile y & S o n s , N ew Y o rk , 1 9 9 3 ) in th e fo rm

xCxCxCxCxV cossincoshsinh)( 4321

(6 ) w h ere co n s tan ts 41 to CC d ep en d o n b o u n d ary co n d itio n s . A p p ly in g b o u n d a ry co n d itio n s fo r a th in can tilev e r b eam (c lam p ed – free ) w e ge t a

freq u en cy d e te rm in an t [0, 1, 0, 1 ] [lam, 0, lam, 0 ] [sinh(lam*L)*lam^2, cosh(lam*L)*lam^2, -sin(lam*L)*lam^2, -cos(lam*L)*lam^2] [cosh(lam*L)*lam^3, sinh(lam*L)*lam^3, -cos(lam*L)*lam^3, sin(lam*L)*lam^3]

F r o m t h e c o n d i t i o n t h a t t h e f r e q u e n c y d e t e r m i n a n t i s e q u a l t o z e r o w e g e t t h e f r e q u e n c y e q u a t i o n i n h e f o r m

01coscosh LL

R o o t s o f t h i s e q u a t i o n c a n o n l y b e f o u n d n u m e r i c a l l y , D e n o t i n g Lx ii w e g e t t h e

n a t u r a l f r e q u e n c i e s i n t h e f o r m

...,3,2,12

iE

A

I

L

x ii

Comparison of analytical and FE results counter continuum frequencies FE frequencies 1 5.26650 4690912090e+002 5.26650 9194371887e+002 2 3.300 462151726965e+003 3.300 571391657554e+003 3 9.24 1389593048039e+003 9.24 3742518773286e+003 4 1.81 0943523875022e+004 1.81 2669270993247e+004 5 2.99 3619402962561e+004 3.00 1165614576545e+004 6 4.4 71949023233439e+004 4.4 96087393371327e+004 7 6. 245945376065551e+004 6. 308228786109306e+004 8 8. 315607746908118e+004 8. 451287572802173e+004 9 1.0 68093617279631e+005 1.0 92740977881639e+005

0 5 10 15 200

1

2

3

4

5

6

7

8

9x 10

5 o - FEM, x - analytical

counter

angu

lar

frequ

enci

es

0 5 100

0.5

1

1.5

2

2.5relative errors for FE frequencies [%]

counter

Are analytically computedfrequencies exact to be usedas an etalon for error analysis?

To answer this you have to recallthe assumptions used for thethin beam theory

FE computation with10 beam elementsConsistent mass matrixFull integration

Comparison of Beam and Bilinear Elements Used for Cantilever Beam Vibration

• 10 beam elements … the ninth bending frequency with 2.5% error

• 10 beam elements … this element does not yield axial frequencies

• 10 bilinear elements … the first bending frequency with 20% error, the errors goes down with increasing frequency counter

• 10 bilinear elements … the errors of axial frequencies are positive for consistent mass matrix, negative for diagonal mass matrix

• Where is the truth?

T ransien t p roblem s in linear dynam ics, no dam ping

tPKqqM

M odelling the 1D w ave equation

2

2202

2

x

uc

t

u

0 0.5 1-0.5

0

0.5

1

1.5eps t = 1.6

0 0.5 10.4

0.6

0.8

1

1.2

1.4dis

0 0.5 1-0.5

0

0.5

1

1.5vel

L1 cons 100 elem Houbolt (red)0 0.5 1

-100

-50

0

50

100acc

Newmark (green), h= 0.005, gamma=0.5

Rod elements used here, the results depend on the method of integration

Classical Lamb’s problem

ra d ia l

B

A C

1 m

1 m

a xia l

p 0

T im p

t

L or Q

ElementsL, Q, full int.Consistent massaxisymmetricMeshCoarse 20x20Medium 40x40Fine 80x80Newmark

Loading a point force equiv. pressure

Example of a transient problem

Axial displacements for point force loading

Coarse, L

Fine, Q

Fine L and medium Q

Medium L and coarse Q

Axial displacements for pressure loading

0 1 2

x 10-4

0

2000

4000

6000

8000

10000

12000

Point and pressure loading of the solid cylinder by a rectangular pulse in time

Time [s]

Tot

al e

ner

gy

[J]

point Q fine

point Q medium

point L fine

point Q coarse

point L medium

point L coarse

FEM-models with the concentrated point-load

For all FEM-models with distributed loads

P-Lcoarse

P-Lmedium

P-Lfine

P-Qcoarse

P-Qmedium

P-Qfine

Total energy

Pollution-free energy production by a proper misuse of FE analysis