O metod ě konečných prvků Lect_ 1 5

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

M. Okrouhlík

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

FFT and FEM

Obsah

• Terminologie

• Fourierova řada

• CFT, DFT a FFT

• Pár příkladů

• Mudrování o platnosti MKP výsledků

TerminologieFourierova řada, CFT, DFT a FFT

• Fourier series

• Continuous Fourier Transform

• Discrete Fourier Transform

• Fast Fourier Transform

Fourierova řada, 20 členů T = Tmax = n * Timp, n = 2,4,8,16

Obdélníkový puls – CFTtři obdélníkové pulsy různé délky – různá normalizace

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0

1

2

3

4

5x 10

-5 C() = 2*abs(sin(T/2))/

Am

plit

ud

e s

pe

ctru

m

timp = 1e-5

timp = 2e-5timp = 5e-5

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0

0.2

0.4

0.6

0.8

1C()/T = 2*abs(sin(T/2))/(T)

frequency in Hz

Am

plit

ud

e s

pe

ctru

m

f3 3 fft c3.m, f ig1

timp = 1e-5


FFT je způsob výpočtu DFT

Základní vztahy_1

Sampling rate,vzorovací frekvence

0 5 10 15 20 25 30 35 40-20

-15

-10

-5

0

5

10

15

20

Continuous and sampled forms of a signal

0 5 10 15 20 25 30 35 40-20

-15

-10

-5

0

5

10

15

20Coarse sampling spoils the signal

timecontinuous and sampled fig1

Základní vztahy_2

Pro diskusi o významu termínu power spectrum jsou následující důležité vztahymean(y) = sum(y)/NArea = tmax*sum(y)/(N–1) = tmax*mean(y)*N/(N-1) ~ tmax*mean(y)

první členy řad Y a pYY(1) = mean(y) * NpY(1) = mean(y) * mean(y) * N

final

initial

dd0i0

t

t

t ttPtetPc

Je to vztah mezi číselnými hodnotami

Základní vztahy_3

Obdélníkový pulsFFT (Fast Fourier Transform) a MPS (Matlab Power Spectrum)

0 1 2 3 4 5 6 70

1

2

3

4sampled signal

time [s]

0 10 20 30 40 50 60-10

0

10

20FFT spectrum - real part

counter

Nyquist

0 10 20 30 40 50 60-10

0

10

20FFT spectrum - imaginary part

counter

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 50

1

2

3

4

5FFT power spectrum

frequency [Hz]

Nyquist

kkk c1 f ig1

Význam členů y(1) a Y(1)

Normalization of FFT data

Základní vztahy_4

Základní vztahy_5

Nyquistova frekvence závisí na tom, zda počet vzorků je sudý či lichý.

Je-li počet vzorků velký, rozdíl nestojí za řeč.

function [f, p_x, n_f] = power_spect_c2(t,x) % compute fft power spectrum for x signal in time domain % input % t ... time ..... t(i) ... i = 1 : N % time range ..... t = 0 : dt : tmax % x ... signal ... x(i) ... i = 1 : N % output % f frequency f(i) ..... i = 1 : NF % p_x complete fft power spectrum p_x(i) ... i = 1 : N % n_f ... Nyquist frequency % number of sampled points N = length(t); % dt .. time increment, assumed uniform dt = t(2) - t(1); % sampling frequency = sampling rate s_r = 1/dt; % Nyquist frequency n_f = 0.5*s_r; % is N odd or even? use the reminder function rm = rem(N,2); % calculate increment of frequencies and frequencies if rm == 0,, df = s_r/N; f = (0:N/2)*df; % N is even else df = s_r/(N-1); f = (0:(N-1)/2)*df; % N is odd end NF = length(f);

% Fourier spectrum f_x = fft(x); % complete power spectrum p_x = f_x.*conj(f_x)/N; % the first half of the power spectrum p_x = p_x(1:NF);

% end of power_spect_c2

Obdélníkový puls – CFTtři obdélníkové pulsy různé délky – různá normalizace

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0

1

2

3

4

5x 10

-5 C() = 2*abs(sin(T/2))/

Am

plit

ud

e s

pe

ctru

m

timp = 1e-5


0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

x 105

0

0.2

0.4

0.6

0.8

1C()/T = 2*abs(sin(T/2))/(T)

frequency in Hz

Am

plit

ud

e s

pe

ctru

m

f3 3 fft c3.m, f ig1

timp = 1e-5


CFT vs. FFTthree different rectangular pulses, the same sampling and the same window

0 0.2 0.4 0.6 0.8 1

x 10-4

0

0.2

0.4

0.6

0.8

1

time

do

ma

in

time [s]

timp = 0.1

0 0.2 0.4 0.6 0.8 1

x 10-4

0

0.2

0.4

0.6

0.8

1

timp = 0.2

time [s]0 0.2 0.4 0.6 0.8 1

x 10-4

0

0.2

0.4

0.6

0.8

1

timp = 0.5

time [s]

0 0.5 1 1.5 2

x 105

0

1

2

3

4

5

ma

tlab

po

we

r sp

ect

rum

0 0.5 1 1.5 2

x 105

0

5

10

15

20

0 0.5 1 1.5 2

x 105

0

50

100

0 0.5 1 1.5 2

x 105

0

10

20

30

mp

s/m

ea

n/ti

mp

frequency [Hz]F3 3 fft c3.m, f ig2 ... the same w indow and sampling rate, three different pulses

0 0.5 1 1.5 2

x 105

0

10

20

30

frequency [Hz]0 0.5 1 1.5 2

x 105

0

10

20

30

frequency [Hz]

1/dt1/timp

CFT vs. FFTthe same rectangular pulse, the same sampling and three different windows

0 0.2 0.4 0.6 0.8 1 1.2

x 10-4

0

0.2

0.4

0.6

0.8

1

tmax = 2*timp

time [s]

rect

angu

lar

puls

e

0 0.2 0.4 0.6 0.8 1 1.2

x 10-4

0

0.2

0.4

0.6

0.8

1

tmax = 4*timp

time [s]0 0.2 0.4 0.6 0.8 1 1.2

x 10-4

0

0.2

0.4

0.6

0.8

1

tmax = 8*timp

time [s]

0 1 2 3 4 5

x 105

0

20

40

60

80

pow

er s

pect

rum

frequency [Hz]

Nyquist = 5000000 [Hz]

0 1 2 3 4 5

x 105

0

10

20

30

40

frequency [Hz]


0 1 2 3 4 5

x 105

0

5

10

15

20

frequency [Hz]


0 1 2 3 4 5

x 105

0

0.1

0.2

0.3

0.4

0.5

norm

aliz

ed p

ower

spe

ctru

m

frequency [Hz]0 1 2 3 4 5

x 105

0

0.1

0.2

0.3

0.4

0.5

frequency [Hz]0 1 2 3 4 5

x 105

0

0.1

0.2

0.3

0.4

0.5

frequency [Hz]

mean1 = 0.5 mean3 = 0.125mean2 = 0.25

c02 = mean22 * N2c01 = mean12 * N1 c03 = mean12 * N3

df1

df2

Fyzikální význam Denoting the discretized time evolution of reference pulse T

21 nPPPP we can express its

mean value

n

in 1mean

1iPP and the approximate the area under the pulse

final

initial

dt

t

ttP by TmeanP .

The physical meaning of the last expression is the mean impulse. The area under the

pulse can also be expressed by means of the first term of the continuous Fourier’s transform,

CFT [Doyle, 1997]. This is the term corresponding to zero-th frequency and has the form

final

initial

dd0i0

t

t

t ttPtetPc .

Pokud je tedy vstupní veličinou síla,

pak power spectrum (výkonové)

by se mělo nazývat

momentum (hybnostní) spectrum

Momentum - hybnost

Four different input pulses and their average

The correlation between 0c and meanP values is shown in the first column of subfigures in

Fig. 12. They are plotted as functions of the pulse counter. Both of them are readily available at

the time of experiment and bear information related to the total momentum since the impulse of

the loading force is equal to the rate of momentum, i.e. vmtP dd . Observing these two

quantities one can easily judge the magnitude of the total impulse of the loading force.

The second couple of quantities that could be correlated (see the subfigures in the right-hand side column of Fig. 12) are the Euclidean norm of the time history of the loading

pulse

n

i 1

2iP and the square root of the total input energy. This follows from the fact

that the energy is proportional to the quadratic form of displacements and velocities, i.e. qMqKqq TT and could also be expressed as proportional to the quadratic form of

applied forces TT QMQQKQ 1T1T , where qq , are global arrays of displacements and velocities ‘measured’ at time finalt , KM, are mass and stiffness

matrices, T is the time interval during which the mean global loading vector, defined by T

mean 000000 PQ , is applied. So observing the norm value one can estimate

amount of energy being delivered by the pulse to the mechanical system.

Tube with four spiral slots, coarse mesh, dimensions

Four different meshes

Comparison of FE and EXP data for axial strains in different locationsC:\spiral_slot_2006_backup_from_2007_10_24\1_u_2007_december_changed_data\u_2007_data_from_urmas_march_2007\ToMila\plot_urmas_data_march_2007_c1, f6

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10

-4 LOC1 - exp vs. FE 3D, NM

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


axia

l str

ain

time

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


time

experiment

FE analysis

experiment

FE analysis

experiment

FE analysis

Comparison of filtered FE and EXP data for axial strains in different locations

C:\spiral_slot_2006_backup_from_2007_10_24\1_u_2007_december_changed_data\u_2007_data_from_urmas_march_2007\ToMila\plot_urmas_data_march_2007_c1, f5

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10

-4 LOC1 - exp vs. FE 3D, NM - filtered to experimental frequency cutoff

G1R90L, mesh1, medium striker for all figuresurmas data march 2007

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


axia

l str

ain

time

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


time

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)

experiment

FE analysis, filtered

experiment


experiment


FE data were filtered by a filter having the same upper frequency as the experiment sampling frequency

And this leads to a question.

That is, up to which frequency limit is the FE approach trustworthy?

We know that FE method is a model of continuum. Continuum – also a model – being based on the continuity hypothesis, disregards the corpuscular structure of matter. It is assumed that matter within the observed specimen is distributed continuously and its properties do not depend on the specimen size.

Quantities describing the continuum behavior are expressed as piecewise continuous functions of time and space. It is known, see [2], that such a conceived continuum has no upper frequency limit.

To find a ‘meaningful’ frequency limit of FE model, which is of discrete – not continuous – nature, one might pursue the following heuristic reasoning.

Rule of thump1mm corresponds to 1 MHz

Easy to rememberImagine a uniform finite element mesh with a characteristic element size, say s . Trying to safely ‘grasp’ a harmonic component (having the wavelength ) by this element size we require that at least five-element length fits the wavelength. This leads to s5 . What is the frequency of this harmonics? Taking a typical wave speed value in steel of about m/s5000c and realizing that cT and Tf /1 , we get the sought-after ‘frequency limit’ in the form )5/( scf . For a one-millimeter element we

get MHz 1 Hz101001.05

5000 6

f . Let’s call it the five-element frequency, denoting

it elem5f in the text.

Bathe recommends 10 here

When looking for the upper frequency limit of a discrete approach to continuum problems,

we could proceed as follows

• 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 have

s

s5

m/s5000ccT

Tf /1)5/( scf

MHz 1 Hz101001.05

5000 6

f

10-10

10-8

10-6

10-4

10-2

10-2

100

102

104

106

108

size in [m]

freq

uenc

y in

[Mhz

]characteristic sizes and corresponding frequencies

atom sizeaustenite steel grain size1 mm finite element1 MHz level1 GHz levelFE analysis range from 0.1mm to 100mmmaximum exp. sampling limit 100 MHz - 14 bits

Where is the continuum limit?

Limits of continuum, FE analysis and experiment

All considered material properties within the observed infinitesimal element are identical with those of a specimen of finite size

10 MHz

10 kHz

Back to our example

Observing the original (or raw) FE signal in Fig. 4 we may notice its three significant characteristics. First, the negative peak representing the input rectangular pulse, as it was changed on its way from the loading face of the tube to the measurement location; second, the slow frequency variation of the tail of the signal and finally the high frequency components superimposed on the signal everywhere. To estimate the low frequencies, appearing in the signal, let’s consider the lowest frequency of the unsupported infinitely long thin shell of the radius r . In [21] there is

derived the formula)21)(1(

1

2

1

E

rf , which when applied to our case gives

the value of 93 kHz. Due to the corresponding mode of vibration, let’s call this frequency the lowest breathing frequency.

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


axia

l str

ain

time

0 0.5 1 1.5 2 2.5 3 3.5

x 10-4

-4

-2

0

2x 10


time

experiment

FE analysis

experiment

FE analysis

experiment

FE analysis

The first eigenfrequency of an infinitely long thin-walled tubeBreathing frequency

.

The tube has the thickness h, which is substantially smaller than its radius r. The material constants are Young’s modulus E, density and the Poisson’s ratio . Assuming that the symmetry, radial and circumferential axes are denoted x, y and z respectively, we can write

zz

xy

yy

xx

zz

xy

yy

xx

E

10

02/)21(00

01

01

)21)(1(

Let’s assume that in this case only radial displacement yu is non-zero.

Assume the free body diagram of the mass element whose elementary mass is zrhm ddd . The radial displacement evokes stresses zz acting on both sides of the

element having the surface area hdz. The corresponding elementary forces are zhP zz dd . The equation of motion for the radial direction with the resultant of P-forces is

0)2/dsin(d2d Pum y .

Furthermore

.2/d)2/dsin(;;)21)(1(

)1(

r

uE yzzzzzz

Assuming the harmonic motion in the form t

y Uu ie and its second time derivative,

substituting this and the previous relations into the equation of motion, and rearranging we get

)21)(1(

11

E

r.

Zig-zag frequencyThe faster frequency appearing in time distributions of displacements and strains is called the zig-zag frequency in the text.

For the zig-zag frequency estimation let’s pursue the following reasoning. The case we are dealing with is three-dimensional even if its axial dimension is predominant and the thickness of the tube is rather small comparing to its axial length. Also the applied loading is rather mild – meaning that the time length of the pulse is relatively long with respect to time needed for a wave to pass through the overall length of the tube. Still, in reality there is a fully 3D wave motion pattern appearing within the tube cross section that is dutifully detected by the FE modeling we are employing. According to Huygens’ Principle each point on the surface being hit by a wave is a source of two kinds of waves – the longitudinal and transversal (shear) waves, respectively. [3].

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

At the beginning of the loading process the frequencies of evoked waves can be crudely estimated the following way. Each type of wave, being emanated from the outer surface, propagates through the tube thickness, is reflected from the inner surface, and hits the outer surface after the time interval

LL

TS

2,

2

c

st

c

st ,

where the tube thickness is denoted by s . The process is repeated. The corresponding estimates of frequencies of S- and L-waves hitting the outer surface are

LL

SS

1,

1

tf

tf .

Considering the given geometry and material properties the numerical values for these frequencies are

MHz93.0,MHz54.0 LS ff .

In the text we will call them zig-zag frequencies with attributes L (for longitudinal waves) and with S (for shear waves) respectively.

Transfer function for mesh1. NM vs.CD. Limit frequencies.C:\miok_2007\spiral_slot_2006_backup_from_2007_10_24\axisym\mesh_refinement_medium_striker_nm\ check_data_10.m, f5

0 1 2 3 4

x 10-5

0

0.5

1

1.5

2

2.5x 10

4 imput pulse

time0 1 2 3 4

x 10-5

0

1

2

3

x 10-5 mesh1, disp, corner node

time

radial NMaxial NMradial CDaxial CD

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5

x 106

0

2

4

6x 10

-10

frequency from 0 to Nyquist

transfer function, input vs. rad. disp

NMCDbreathing frequencyzig-zag L frequencyzig-zag S frequencyFE limit frequency consFE limit frequency diag5elem frequency

Transfer functions for different meshes from 0 to NyquistC:\miok_2007\spiral_slot_2006_backup_from_2007_10_24\axisym\mesh_refinement_medium_striker_nm\ check_data_10.m, f13

0 1 2 3 4 5

x 106

0

1

2

3

4

5

6x 10

-10

tra

nsf

er

fun

ctio

n

axisym, corner, rad. disp, mesh1

0 2 4 6 8 10

x 106

0

0.5

1

1.5

2

2.5

3x 10

-5 axisym, corner, rad. disp, mesh2

0 0.5 1 1.5 2

x 107

0

0.5

1

1.5

x 10-5 axisym, corner, rad. disp, mesh3

frequency from 0 to Nyquist [Hz]

tra

nsf

er

fun

ctio

n

0 1 2 3 4

x 107

0

0.2

0.4

0.6

0.8

1

1.2

x 10-5 axisym, corner, rad. disp, mesh4

frequency from 0 to Nyquist [Hz]

NMCDbreathing frequencyzig-zag L frequencyzig-zag S frequencyFE limit frequency consFE limit frequency diag5elem limit frequency




Transfer functions for different meshes from 0 to 2 MHzC:\miok_2007\spiral_slot_2006_backup_from_2007_10_24\axisym\mesh_refinement_medium_striker_nm\ check_data_10.m, f14

0 0.5 1 1.5 2

x 106

0

1

2

3

4x 10

-10

tra

nsf

er

fun

ctio

n

axisym, corner, rad. disp, mesh1

0 0.5 1 1.5 2

x 106

0

0.2

0.4

0.6

0.8

1x 10


0 0.5 1 1.5 2

x 106

0

0.2

0.4

0.6

0.8

1x 10


frequency range from 0 to 2 MHz

tra

nsf

er

fun

ctio

n

0 0.5 1 1.5 2

x 106

0

0.2

0.4

0.6

0.8

1x 10


frequency range from 0 to 2 MHz

NMCDbreathing frequencyzig-zag L frequencyzig-zag S frequency




Observing Fig. 7 one should notice the subsequent ‘convergence’ of CD and NM peaks within the zig-zag frequency interval. The natural explanation is that with the finer meshsize, and with the correspondingly smaller timestep, both methods operate in ‘good’ frequency intervals where their spatial and temporal discretization errors are minimized.

The presented transfer spectra for four studied meshes show a distinct indication of the breathing and zig-zag frequencies, the ‘convergence’ of CD and NM responses, subsequent disappearance of ‘false’ CD responses and that the ‘dubious’ CD frequency peaks do not have their counterparts in NM

responses.

What remains to be compared is the ‘raw’ FE signal with that the frequencies higher than the five-element frequency were filtered out. The results for the ‘raw’ and filtered FE signals, for the mesh1 and the NM operator with consistent mass matrix, are presented in Fig. 8.

FE raw signal compared to that in which the frequencies higher than five-element ones were filtered out.

C:\spiral_slot_2006_backup_from_2007_10_24\1_u_2007_december_changed_data\u_2007_data_from_urmas_march_2007\ToMila\plot_urmas_data_march_2007_c1.m, f10

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

x 10-4

-4

-3

-2

-1

0

1x 10

-4 LOC1 - FE 3D, mesh1, NM raw signal vs. filtered to FE 5elem cutoff

time [s]

axi

al s

tra

in

FE analysis, filtered to 5elem cutoffFE analysis, raw data

In future these results might be confirmed by a new and more sophisticated experiment having a lower observational threshold, a higher sampling rate and also a higher frequency amplifier cut-off .

Conclusions

The FE analysis is a robust tool giving reliable results with a satisfactory engineering precision in standard tasks of continuum mechanics.

Nevertheless, employing the FE method in cases on borders of their applicability is tricky and obtained results have to be treated with utmost care, since they might be profoundly influenced by intricacies of finite element technology.

It should be emphasized, however, that testing the methods in the vicinity of borders of their applicability we do not want to discredit them, on the contrary, the more precise knowledge of their imperfections makes us – users – more confident in them.

Date post:	14-Jan-2016
Category:	Documents
Upload:	lynnea
View:	36 times
Download:	0 times

O metod ě konečných prvků Lect_ 1 5

Documents