fft.mws

---

Solving a 1-DOF System Driven by Transient and Non-Periodic Forces

by Harald Kammerer, Germany, maple@jademountain.de

This worksheet compares the solutions of a single-degree-of-freedom-system derived by different techniques

KEYWORDS: FFT, frequency domain, time domain, DOF

The author expects that this worksheet will only be used for teaching and educational purposes and not for commercial profit without contacting the author for a licensed agreement.

>

restart;

>	with(plots):with(plottools):

Warning, the name changecoords has been redefined

Warning, the name arrow has been redefined

Introduction

We consider a simple mechanical system with one degree of freedom as shown in the following figure.

The system is loaded by a ground motion, described by the acceleration . In the following we consider different acceleration functions, described by the variables  and . The first one is a simple impulse motion, the second one is a series of impulses.

The motion of the mass m is described by the displacement  and the acceleration . The mass is supported on a combination of a spring with stiffness k  and a damper with the damping resistance d .

To solve the problem we use two different usual methods:

1. solve the equation of motion in the time domain

2. convert the ground acceleration into the frequency domain, multiply the result with the transfer function of the system and convert the solution back to the time domain

The goal is to show the difference between the results of both methods and describe the reason for this difference. The background for the idea for this example is , that in many cases people use the solution in the frequency domain to solve problems with transient excitation with non-periodic excitations. This can result in failures.

To convert the function from the time domain to the frequency domain and vice versa we use the Fast-Fourier-Transformation (see FFT ). To use this function we need a discrete number of points of the functions . This number must be  with p  a non-negative integer. Additionally we define here one more variable n2:

>	p:=8; n:=2**p; n2:=n/2+1;

Description of The Ground Motion

For the first case, the ground motion is given by a single impulse with the duration

>

T1:=0.5;

The ground acceleration is given by

>	spp0(t):=piecewise(t<T1/2 or t>3*T1/2,0,1);

Second we consider a series of impulses with the same duration. The time distance between the impulses is

>	Tmax:=1.0;

The ground motion is now described by the acceleration

>	sppa(t):=(1-cos(2*Pi*t)/abs(cos(2*Pi*t)))/2;

In the following figures these functions are shown

>	plot(spp0(t),t=0..5*Tmax,color=red,title="Excitation by a single impulse");

>	plot(sppa(t),t=0..5*Tmax,color=black,title="Excitation by a series of impulses");

Fourier Transform of the Ground Motion

To use the Fast-Fourier-Transformation (FFT) we need numerical values of the acceleration in discrete time steps. Further the FFT assumes that the function which will be transformed is periodic. To use the FFT for solving non-periodic systems we assume that the function is periodic.

We will see now what happens with the single-impulse ground motion after transforming into the frequency domain by use of the FFT. Therefore we calculate values of the function  for n  discrete time steps inside the time interval [0,1).

>	for i from 1 by 1 to n do

>	ti[i]:=evalf((i-1)/n*Tmax):

>	x[i]:=evalf(subs(t=ti[i],spp0(t))):

>

od:

This values are stored in the array X. Additionally we define the array Y, which contains n  times the value 0, because the function  is a real function.

>	X:=array([seq(x[i],i=1..n)]):

>	Y:=array([seq(0,i=1..n)]):

The next figure shows the n  discrete function points.

>	display(curve([seq([ti[i],x[i]],i=1..n)]),style=point,color=red,title="Discrete values of the excitation function");

The Fourier transform of the function  is now calculated by use of the Maple-function FFT

>	FFT(p,X,Y);

The result is shown in the following figure. The black dots show the real part and red dots show the imaginary part.

>	display(curve([seq([i,X[i]],i=1..n)],style=point),title="FFT-result of the excitation, real part");

>	display(curve([seq([i,Y[i]],i=1..n)],color=red,style=point),title="FFT-result of the excitation, imaginary part");

This result is not the final solution of the transformation. To get this there are additional steps necessary.

The FFT yields the Fourier transform for discrete frequency values which are given by the time interval for which the function is given. The distance of this discrete values for the frequency is given by

>	df:=1/Tmax;

The values which are yielded by the function FFT have to be divided by n . So the constant parts are:

>	xs[1]:=X[1]/n;

>	ys[1]:=Y[1]/n;

In the considered cases of Fourier transforms of real functions, the constant part of the imaginary values is always 0.

For i from 2 up to n2 the X[i]/n and Y[i]/n are the amplitudes of the harmonic functions to the frequency . The values of X[i]/n and Y[i]/n for i from n2+1 to n are the amplitudes for the negative values for the frequency, that is . But we know that the real part, the X[i]  belongs to the cosine function and the Y[i] belongs to the sine function. So we use the relation

 and

, so we get only  pairs of harmonic terms

>	for i from 2 by 1 to n2 do

>	xs[i]:= 2*X[i]/n:

>	ys[i]:= 2*Y[i]/n:

>

od:

The discrete values for the frequencies are

>	for i from 1 by 1 to n2 do

>	f[i]:=(i-1)*df:

>

od:

Now we have the excitation transformed to series of harmonic terms:

.

The following figures show the  amplitudes of the amplitudes of the harmonic terms of the Fourier transform of the excitation

>	display(curve([seq([f[i],xs[i]],i=1..n2)],color=blue,style=point),title="Amplitudes of the cosine terms");

>	display(curve([seq([f[i],ys[i]],i=1..n2)],color=red,style=point),title="Amplitudes of the sine terms");

So the excitation is described as Fourier series by:

Sum of cosine terms:

>	xfc(t):= sum(xs[ij]*cos(2*Pi*f[ij]*t),ij=1..n2):

and sum of sine terms

>	xfs(t):=-sum(ys[ij]*sin(2*Pi*f[ij]*t),ij=1..n2):

The result is the sum of both parts

>	xf(t):=xfc(t)+xfs(t):

The solution is shown in the following figures  

>	display(plot(xfc(t),t=0..Tmax,color=blue),title="Fourier series, blue: cosine");

>	display(plot(xfs(t),t=0..Tmax),title="Fourier series, red: sine");

>	display(plot(xf(t),t=0..Tmax,color=black,thickness=2),title="Fourier series, black: cosine+sine");

You see that this series doesn't describe the actual function exactly. It is only an approximation. The greater the number n the better the approximation. Let's now compare the result with the actual ground acceleration at the same discrete time steps as used above to get this Fourier series.

>	for i from 1 by 1 to n do

>	xfi[i]:=evalf(subs(t=ti[i],xf(t))):

>

od:

The next figure show this discrete values

>	display(plot(spp0(t),t=0..Tmax,color=black),curve([seq([ti[i],xfi[i]],i=1..n)],color=red,style=point), title="black: actual excitation, red: Fourier series (discrete values)");

We have seen that the Fourier transform in reality only approximates periodic functions. In the following figure we see that the calculated Fourier series is not really an approximation for the single impulse excitation  but for the series of impulses .

>	display(plot(spp0(t),t=0..5*Tmax,color=black, thickness=2),plot(xf(t),t=0..5*Tmax,color=red),title="black: actual excitation, red: Fourier series");

Description of the Mechanical System

Here we consider a concrete system with the following characteristic

>	m:=2.*10**3;

>	k:=3.*10**6;

>	d:=5.*10**3;

The natural frequency and the damping ratio are

>	evalf(sqrt(k/m)/2/Pi);

>	evalf(d/(2*m*sqrt(k/m)));

The equation of motion for the considered system is

>	bgl0:=m*(diff(s(t),t$2)+a0(t))+d*diff(s(t),t)+k*s(t)=0;

with the general ground acceleration . Now we consider some different types for the solution.

Solution with transient vibration, excitation described by the Fourier series

First the constant part of the excitation is considered.

>	sppi[1]:=xs[1]-ys[1]:

This part is substituted into the equation of motion

>	bgl[1]:=subs(a0(t)=sppi[1],bgl0):

The solution is reached by use of the Maple function dsolve. The motion should start from the rest, so the initial conditions are

>	IC:=s(0)=0,D(s)(0)=0:

>	lsg[1]:=dsolve({bgl[1],IC},s(t)):

And the  relative displacement of the system caused by this part of the excitation is

>	ski[1]:=value(rhs(lsg[1])):

In the same manner we get the displacement for every term of the Fourier series for the ground acceleration

>	for i from 2 by 1 to n2 do

>	sppi[i]:=xs[i]*cos(2*Pi*f[i]*t)-ys[i]*sin(2*Pi*f[i]*t):

>	bgl[i]:=subs(a0(t)=sppi[i],bgl0);

>	lsg[i]:=dsolve({bgl[i],s(0)=0,D(s)(0)=0},s(t)):

>	ski[i]:=value(rhs(lsg[i])):

>

od:

The total relative displacement between the body and the ground is finally the sum of all terms

>	sk(t):=evalf(sum(ski[ij],ij=1..n2)):

The relative velocity of the body is

>	vk(t):=  diff(sk(t),t):

The relative acceleration is

>	ak(t):=  diff(sk(t),t$2):

The absolute acceleration is

>	aak(t):= diff(sk(t),t$2)+xf(t):

Solution with transient vibration, excitation described by the actual ground acceleration

Now we consider the case that the excitation is the single impulse ground acceleration. The equation of motion is the described by

>	bgla:=convert(subs(a0(t)=spp0(t),bgl0),Heaviside):

And the solution is

>	lsga:=dsolve({bgla,s(0)=0,D(s)(0)=0},s(t)):

>	ska(t):=value(rhs(lsga)):

The relative velocity of the body is

>	vka(t):= diff(ska(t),t):

The relative acceleration is

>	aka(t):= diff(ska(t),t$2):

The absolute acceleration is

>	aaka(t):=diff(ska(t),t$2)+spp0(t):

Solution without transient vibration (steady state motion), excitation described by the Fourier series

the method to get the solution is nearly the same. But here we don't need any initial conditions, only the integration constants have to be set to zero. This yields the steady state motion. This is again first made for the constant part of the excitation.

>	sppi[1]:=xs[1]-ys[1]:

>	bgk[1]:=convert(subs(a0(t)=sppi[1],bgl0),Heaviside):

>	sli[1]:=rhs(subs({_C1=0,_C2=0},dsolve(bgk[1],s(t)))):

And then the same way is used for all the other terms of the Fourier series

>	for i from 2 by 1 to n2 do

>	sppi[i]:=xs[i]*cos(2*Pi*f[i]*t)-ys[i]*sin(2*Pi*f[i]*t):

>	bgk[i]:=convert(subs(a0(t)=sppi[i],bgl0),Heaviside);

>	sli[i]:=rhs(subs({_C1=0,_C2=0},dsolve(bgk[i],s(t)))):

>

od:

The total relative displacement between the ground and the body is the sum of all terms of the solution.

>	sl(t):=evalf(sum(sli[ji],ji=1..n2)):

The relative velocity of the body is

>	vl(t):=  diff(sl(t),t):

The relative acceleration is

>	al(t):=  diff(sl(t),t$2):

The absolute acceleration is

>	aal(t):=  diff(sl(t),t$2)+xf(t):

Solution without transient vibration (steady state motion), excitation described by the actual ground acceleration

The solution of the equation of motion is now

>	sla(t):=rhs(subs({_C1=0,_C2=0},dsolve(bgla,s(t)))):

The relative velocity of the body is

>	vla(t):= diff(sla(t),t):

The relative acceleration is

>	ala(t):= diff(sla(t),t$2):

The absolute acceleration is

>	aala(t):= diff(sla(t),t$2)+spp0(t):

At last we control whether the result fulfill the equation of motion

>	Ek(t):=  (m*diff(ska(t),t$2)+d*diff(ska(t),t)+k*ska(t))/m:

>	El(t):=  (m*diff(sla(t),t$2)+d*diff(sla(t),t)+k*sla(t))/m:

Plot the results

The results are not so easy to show in simple manner so we show them graphical.

Displacement

>	P[1]:=plot((sk(t)),t=0..Tmax,color=black,title="Displacement with transient motion (Fourier)"):

>	P[2]:=plot((ska(t)),t=0..Tmax,title="Displacement with transient motion (actual excitation)"):

>	P[3]:=plot((sk(t)),t=50*Tmax..51*Tmax,color=black,title="Displacement with transient motion (Fourier)"):

>	P[4]:=plot((ska(t)),t=50*Tmax..51*Tmax,title="Displacement with transient motion (actual excitation)"):

>	P[5]:=plot((sl(t)),t=0..Tmax,color=black,title="Displacement, steady state motion (Fourier)"):

>	P[6]:=plot((sla(t)),t=0..Tmax,title="Displacement, steady state motion (actual excitation)"):

Velocity

>	P[7]:=plot((vk(t)),t=0..Tmax,color=black,title="Velocity with transient motion (Fourier)"):

>	P[8]:=plot((vka(t)),t=0..Tmax,title="Velocity with transient motion (actual excitation)"):

>	P[9]:=plot((vk(t)),t=50*Tmax..51*Tmax,color=black,title="Velocity with transient motion (Fourier)"):

>	P[10]:=plot((vka(t)),t=50*Tmax..51*Tmax,title="Velocity with transient motion (actual excitation)"):

>	P[11]:=plot((vl(t)),t=0..Tmax,color=black,title="Velocity, steady state motion (Fourier)"):

>	P[12]:=plot((vla(t)),t=0..Tmax,title="Velocity, steady state motion (actual excitation)"):

Relative acceleration

>	P[13]:=plot((ak(t)),t=0..Tmax,color=black,title="Acceleration (relative) with transient motion (Fourier)"):

>	P[14]:=plot((aka(t)),t=0..Tmax,title="Acceleration (relative) with transient motion (actual excitation)"):

>	P[15]:=plot((ak(t)),t=50*Tmax..51*Tmax,color=black,numpoints=500,title="Acceleration (relative) with transient motion (Fourier)"):

>	P[16]:=plot((aka(t)),t=50*Tmax..51*Tmax,title="Acceleration (relative) with transient motion (actual excitation)"):

>	P[17]:=plot((al(t)),t=0..Tmax,color=black,title="Acceleration (relative), steady state motion (Fourier)"):

>	P[18]:=plot((ala(t)),t=0..Tmax,title="Acceleration (relative), steady state motion (actual excitation)"):

Absolute acceleration

>	P[19]:=plot((aak(t)),t=0..Tmax,color=black,title="Acceleration (absolute) with transient motion (Fourier)"):

>	P[20]:=plot((aaka(t)),t=0..Tmax,title="Acceleration (absolute) with transient motion (actual excitation)"):

>	P[21]:=plot((aak(t)),t=50*Tmax..51*Tmax,color=black,numpoints=500,title="Acceleration (absolute) with transient motion (Fourier)"):

>	P[22]:=plot((aaka(t)),t=50*Tmax..51*Tmax,title="Acceleration (absolute) with transient motion (actual excitation)"):

>	P[23]:=plot((aal(t)),t=0..Tmax,color=black,title="Acceleration (absolute), steady state motion (Fourier)"):

>	P[24]:=plot((aala(t)),t=0..Tmax,title="Acceleration (absolute), steady state motion (actual excitation)"):

Control

>	P[25]:=plot((Ek(t)),t=0..Tmax,color=black,title="Control (with transient motion, actual axcitation)"):

>	P[26]:=plot((El(t)),t=0..Tmax,title="Control (steady state motion, actual excitation)"):

Now the above defined plots are shown. First we consider the relative displacement. This is the relative motion between ground and body.

>	for i from 1 by 1 to 2 do display(P[i]) od;

Here we see nearly no difference between the solution for the actual excitation and the approximation by the Fourier series. But now we consider the motion after a longer time.

>	for i from 3 by 1 to 4 do display(P[i]) od;

Of cause now there is a great difference, because the actual excitation is a single impulse and yields a fading away vibration after the impulse. In the considered time range the motion nearly vanishes. The Fourier series describes a periodic excitation and the result is a periodic motion. Only while the first period both ground motions yield nearly the same result. When the number of Fourier terms increases, the difference between both motions during this first period decreases.

Next we consider the steady state motion. Of course in the case of the single impulse it makes no sense to show this solution, because during the impulse the motion is a transient motion and the steady state is reached after longer time, depending on the damping ratio.

>	for i from 5 by 1 to 6 do display(P[i]) od;

Now we see that in the case of the excitation by the Fourier series the solution is nearly the same as the solution for the case of consider the solution with transient motion after a longer time. But the "steady state motion" in the case of the solution with the actual ground motion is not the really steady state motion.

So we see: The approximation of the excitation by the Fourier series cannot describe the right solution for a longer time period but only for the time where the excitation is approximated by the Fourier series.

The same result yields the consideration of the relative velocity, the relative acceleration and the absolute acceleration. This plots are shown in the following figures

>	for i from 7 by 1 to 24 do display(P[i]) od;

At last the both control plots are shown

>	for i from 25 by 1 to 26 do display(P[i]) od;

Both plots show the negative excitation, so the control is all right.

Solve the System in the Frequency Domain

First we deduce the general solution for the problem. To do this we consider the general form of the equation of motion:

>	BGL:=M*(diff(Z(t),t$2)+A(t))+R*diff(Z(t),t)+K*Z(t)=0;

To solve the problem in the frequency domain the excitation must be converted to the frequency domain. So we have to use the Fourier series to solve the problem. To get the general solution we consider one harmonic part of the excitation. Further we use the method of complex numbers to go on. Then one harmonic term of the displacement of the ground motion is

>	A(t):=diff(X0*exp(I*Omega*t),t$2);

Here we used the fact that the acceleration  is the second deviation of the displacement.

We know that the displacement of the body has the same mathematical form like the excitation, so we set

>	Z(t):=Z0*exp(I*Omega*t);

Then we get for the equation of motion

>

BGL;

This equation is solved with respect to Z0, this yields

>	LSG:=solve(BGL,Z0);

The transfer function is defined as the relation between the amplitude of the motion of the body and the amplitude of the excitation so we get the

Complex Transfer Function

>	V:=expand(LSG)/X0;

Now we have to do some additional steps to separate the real and the imaginary part of this function. This is here done without further comment, because this is standard math.

(Because the solution has general manner it is here not easy to use the Maple functions "Re" and "Im".)

>	Zaehler:=numer(V);

>	Nenner:=denom(V);

>	Nenner2:=expand((tcoeff(Nenner,I)+coeff(Nenner,I)*I)*(tcoeff(Nenner,I)-coeff(Nenner,I)*I));

>	Zaehler2:=expand(Zaehler*(tcoeff(Nenner,I)-coeff(Nenner,I)*I));

Now we have the imaginary part

>	Imaginaer:=simplify(coeff(Zaehler2,I)/Nenner2);

and the real part of the transfer function

>	Real:=simplify(Zaehler2/Nenner2-Imaginaer*I);

We substitute the concrete values of our problem and get

>	Imaginaerteil:=simplify(subs({M=m,K=k,R=d},Imaginaer));

>	Realteil:=simplify(subs({M=m,K=k,R=d},Real));

At last we substitute the angular frequency  by the frequency  and get

>	C:=subs(Omega=f*2*Pi,Realteil);

>	S:=subs(Omega=f*2*Pi,Imaginaerteil);

The real part and the imaginary part of the transfer function are shown in the following figures dependent on the frequency

>	plot(C,f=0..10,title="Real part of the transfer function");

>	plot(S,f=0..10,title="Imaginary part of the transfer function");

The excitation is with the Fourier transformation with use of the FFT function given for discrete frequency values. So we calculate the transfer function for the same discrete frequencies

>	for i from 1 by 1 to n2 do

>	ci[i]:=evalf(subs(f=f[i],C));

>	si[i]:=evalf(subs(f=f[i],S));

>

od:

>	display(curve([seq([f[i],ci[i]],i=1..n2)],style=point),title="Real part of the transfer function, discrete values");

>	display(curve([seq([f[i],si[i]],i=1..n2)],style=point),title="Imaginary part of the transfer function, discrete values");

We repeat the plot of the Fourier transform of the excitation

>	display(curve([seq([f[i],X[i]],i=1..n2)],style=point),title="Real part of Fourier transfom of the excitation");

>	display(curve([seq([f[i],Y[i]],i=1..n2)],style=point),title="Imaginary part of Fourier transfom of the excitation");

We get the Fourier transform of the system answer by multiply every term of the Fourier transform of the excitation with the corresponding term of the transfer function. See that this terms are complex.

>	for i from 1 by 1 to n2 do

>	USC[i]:=(ci[i]+I*si[i])*(X[i]+I*Y[i]):

>	UC[i]:= evalf(Re(USC[i])):

>	US[i]:= evalf(Im(USC[i])):

>

od:

We want to use the inverse Fast Fourier Transformation ( iFFT ). So we have to store the amplitudes for every harmonic part in arrays of that form like the function FFT yields the result above. The members of the  second part of the arrays (from n2+1 up to n) are

>	for i from n2+1 by 1 to n do

>	UC[i]:= UC[2*n2-i]:

>	US[i]:=-US[2*n2-i]:

>

od:

Now the two arrays which are needed by the function iFFT are defined.

>	UCa:=array([seq(UC[i],i=1..n)]):

>	USa:=array([seq(US[i],i=1..n)]):

>	display(curve([seq([i,UCa[i]],i=1..n)]),title="Real part of the Fourier transform of the motion of the system");

>	display(curve([seq([i,USa[i]],i=1..n)]),title="Imaginary part of the Fourier transform of the motion of the system");

Now iFFT yields the result

>	iFFT(p,UCa,USa):

The following figures show the calculated relative acceleration of the body

>	display(curve([seq([ti[i],USa[i]],i=1..n)]),title="Imaginary part of the relative acceleration");

>	PF:=curve([seq([ti[i],UCa[i]],i=1..n)],color=green,title="Real part of the relative acceleration"):

>	display(PF,title="Real part of the relative acceleration");

We see that the imaginary part of the result nearly vanishes. This is right, because the excitation was real, too. The small values result from numerical errors.

Comparison of the Solutions

We compare this result with the solutions above. The result of the solution in the frequency domain is in the plots with green color.

>	display(P[13],PF);

>	display(P[14],PF);

>	display(P[17],PF);

>	display(P[18],PF);

You see that the only plot with acceptable quality of the conformity is the solution for the steady state motion with approximate the excitation by the Fourier series. The consideration in the frequency domain yields results of bad quality in case of transient motion. (Remember: We have seen above that the "steady state motion" in the last figure is not a real steady state)

Disclaimer:  While every effort has been made to validate the solutions in this worksheet, Waterloo Maple Inc. and the contributors are not responsible for any errors contained and are not liable for any damages resulting from the use of this material.