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

M. Okrouhlík

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

Disperzní vlastnosti

Dispersion

due to• space discretization

• time discretization• geometry

• material

It is known thatthe spatial and temporal discretizations always

accompany the finite element modelling of transient problems,

the topic is of academic interest but has important practical consequences as well.

In the contributionthe spatial dispersive properties of typical

Lagrangian and Hermitian elements will be presented,the dispersion due to time discretization will be

briefly mentioned andexamples of spatial and temporal dispersions

will be shown.

If )(k ... a so called dispersion relation ... is a linear function of k ... the wavenumber, then the system is called NONDISPERSIVE, otherwise it is DISPERSIVE. For non-dispersive systems the velocity of propagation is constant – does not depend on frequency.

L1 … constant strain element

consistent

diagonal

Ec

l

c

kl

kl

kl

00

00

0

*

0*

0

0*

withand where

matrix mass diagonalfor ))cos(1(2

matrix mass consistentfor )cos(2

))cos(1(6

Dispersion for 1D constant strain elements with consistent and diagonal mass matrix formulations

Podrobné odvozenípro L1C a L1D

See: c/telemachos/vyukove texty/disperse_c3.doc

% disp_L1

% dispersive properties of L1 element

clear

i = 0;

gamma_range = 0.01:0.01:pi;

length_gr = length(gamma_range);

for gamma = gamma_range

i = i + 1;

n = (1 - cos(gamma));

n_c = 6*n;

n_d = 2*n;

den = 2 + cos(gamma);

om_c_star(i) = sqrt(n_c/den);

w_c_star(i) = om_c_star(i)/gamma;

om_d_star(i) = sqrt(n_d);

w_d_star(i) = om_d_star(i)/gamma;

lambda(i) = 2*pi/gamma;

end

for i = 1:length_gr

lambda_rev(i) = lambda(length_gr - i + 1);

w_c_star_rev(i) = w_c_star(length_gr - i + 1);

w_d_star_rev(i) = w_d_star(length_gr - i + 1);

om_c_star_rev(i) = om_c_star(length_gr - i + 1);

om_d_star_rev(i) = om_d_star(length_gr - i + 1);

end

% points

figure(1)

subplot(3,1,1); plot(gamma_range,om_c_star, gamma_range,om_d_star)

title('frequency vs. gamma'); axis([0 3.14 0 4])

subplot(3,1,2); plot(gamma_range,lambda)

title('lambda vs. gamma')

v = [0 20 0 4];

subplot(3,1,3); plot(lambda_rev,om_c_star_rev, lambda_rev,om_d_star_rev); axis(v)

title('frequency vs. lambda')

xx = [0 100]; yy1 = [1.01 1.01]; yy2 = [0.99 0.99];

ww = [1 1];

figure(2)

subplot(2,1,1); plot(gamma_range,w_c_star, gamma_range,w_d_star, xx,ww)

title('phase velocity vs. gamma'); axis([0 3.14 0 1.3])

legend('consistent', 'diagonal', 'continuum', 3)

% subplot(3,1,2); plot(gamma_range,lambda)

% title('lambda vs. gamma')

v = [2 15 0.6 1.3];

subplot(2,1,2); plot(lambda_rev,w_c_star_rev, lambda_rev,w_d_star_rev, xx,yy1, xx,yy2); axis(v)

title('phase velocity vs. lambda')

xlabel('lambda is non-dimensional wave length - how many element lengths into the wave length of a given harmonics')

legend('consistent', 'diagonal', '+1%', '-1%', 4)

0 0.5 1 1.5 2 2.5 30

0.2

0.4

0.6

0.8

1

1.2

phase velocity vs. gamma

consistentdiagonalcontinuum

2 4 6 8 10 12 140.6

0.7

0.8

0.9

1

1.1

1.2

1.3phase velocity vs. lambda

lambda is non-dimensional wave length - how many element lengths into the wave length of a given harmonics

consistentdiagonal+1%-1%

0 0.5 1 1.5 2 2.5 3 3.50.6

0.7

0.8

0.9

1

1.1

1.2

phase velocity vs. frequency

consistentdiagonalcontinuum

1D continuum is a non-dispersivemedium, it has infinite number of frequencies, phase velocity is constantregardless of frequency.

Discretized model is dispersive,it has finite number of frequencies,velocity depends on frequency, spectrum is bounded,there are cut-offs.

overestimated consistentFrequency (velocity) is with mass matrix.

underestimated diagonal

L1 element

Eigenvectors, unconstrained rod, 24 L1C elements

Notice the rigid body mode

Higher modes havewrong shapes

Not in accordance with assumed harmonic solution

Nodes in opposition

The system cannot be forced to vibrate more violently

1D – velocity vs. frequency for higher order elements: L2, L3, H3

Band-pass filter

A mixed blessing of a longer spectrum

1D Hermitian cubic element, H3

Eigenvectors of a free-free rod modeled by H3C elements

transversal,shear,S-wave

primary,longitudinal,P-wave

Von Schmidt wavefront

2D wavefronts, Huygen’s principleMaterial points having been hitby primary wave becomesources of both types of wavesi.e. P and S

dp007

2D isotropic medium, plane stress, FE simulation

Dispersive properties of a uniform mesh(plane stress) assembled of equilateral elements, full integrationconsistent and diagonal mass matrices

wavenumberHow many elements to a wavelenght

Direction of wave propagation

Dispersion effects depend alsoon the direction of wave propagation

Artificial (false) anisotropy

Dispersive properties of a uniform mesh (plane stress) assembled of square isoparametric elements, full integration, consistent mass matrix

Hodograph of velocities

meshsize

wavelength

Hodograph of velocities

Hermitian elements are disqualified due to their ‘long’ spectrum

Eigenmodes and eigenvalues of a four-node bilinear element

Eigenmodes and eigenvalues of a four-node bilinear element

Time and space discretization errors vs. frequency

Mother nature is kind to us

No algorithmic damping

Small algorithmic damping

The effect of space and time discretization cancels out

0 0.2 0.4 0.6 0.8 1-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.2

0.4

0.6

0.8

1dis

0 0.2 0.4 0.6 0.8 1-0.2

0

0.2

0.4

0.6

0.8

1

vel

0 0.2 0.4 0.6 0.8 1-50

0

50

100

150

200acc

L1C, 100 elements, h = 0.01, gamma = 0.6, beta = 1/6

‘Dispersionless’ Newmark

E:\edu_mkp_liberec_2\mtl_prog\mtlstep\linbad_c1_dispersionless_newmark.m

0 0.2 0.4 0.6 0.8 1

-1

-0.5

0

eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

dis

0 0.2 0.4 0.6 0.8 1

0

0.5

1

vel

0 0.2 0.4 0.6 0.8 1

-50

0

50

acc

L1C, 100 elements, h = 0.01, gamma = 0.5, beta = 0.25*(0.5+gamma)^2

Newmark, no algorithmic damping

0 0.2 0.4 0.6 0.8 1-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

dis

0 0.2 0.4 0.6 0.8 1-0.2

0

0.2

0.4

0.6

0.8

1

vel

0 0.2 0.4 0.6 0.8 1-30

-20

-10

0

10

20acc

Newmark, with small algorithmic damping

L1C, 100 elements, h = 0.01, gamma = 0.6, beta = 0.25*(0.5+gamma)^2

0 0.2 0.4 0.6 0.8 1-1.2

-1

-0.8

-0.6

-0.4

-0.2

0

eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5

0.6

dis

0 0.2 0.4 0.6 0.8 1-0.2

0

0.2

0.4

0.6

0.8

1

vel

0 0.2 0.4 0.6 0.8 1

-20

-10

0

10

acc

Newmark, with medium algorithmic damping

L1C, 100 elements, h = 0.01, gamma = 0.7, beta = 0.25*(0.5+gamma)^2

0 0.2 0.4 0.6 0.8 1-1

-0.8

-0.6

-0.4

-0.2

0eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

0.5dis

0 0.2 0.4 0.6 0.8 1-0.2

0

0.2

0.4

0.6

0.8

1vel

0 0.2 0.4 0.6 0.8 1-8

-6

-4

-2

0

2

4acc

Newmark, with too big algorithmic damping

L1C, 100 elements, h = 0.01, gamma = 2, beta = 0.25*(0.5+gamma)^2

0 0.2 0.4 0.6 0.8 1-0.8

-0.6

-0.4

-0.2

0eps t = 0.7

0 0.2 0.4 0.6 0.8 10

0.1

0.2

0.3

0.4

dis

0 0.2 0.4 0.6 0.8 1

0

0.2

0.4

0.6

0.8vel

0 0.2 0.4 0.6 0.8 1-3

-2

-1

0

1

2acc

Newmark – s vaničkou jsme vylili i dítě

L1C, 100 elements, h = 0.01, gamma = 10, beta = 0.25*(0.5+gamma)^2