Title: Default for 23.7 An Application of Fourier Series
Course: HELM Workbooks 2023-2024

Modelling vibration by differential equation

Vibration problems are often modelled by ordinary differential equations with constant coefficients. For example the motion of a spring with stiffness \(k\) and damping constant \(c\) is modelled by

\(\qquad m\dfrac{d^2y}{dt^2}+c\dfrac{dy}{dt}+k y=0 \qquad \qquad \qquad \qquad\) (1)

where \(y(t)\) is the displacement of a mass \(m\) connected to the spring. It is well-known that if \(c^2<4mk\), usually referred to as the lightly damped case, then

\(\qquad y(t) = e^{-\alpha t}(A\cos\omega t+B\sin\omega t) \qquad \qquad \qquad\) (2)

i.e. the motion is sinusoidal but damped by the negative exponential term. In (2) we have used the notation

\(\qquad \alpha = \dfrac{c}{2m}\qquad\omega = \dfrac{1}{2m}\sqrt{4km-c^2} \quad\) to simplify the equation.

The values of \(A\) and \(B\) depend upon initial conditions.

The system represented by (1), whose solution is (2), is referred to as an unforced damped harmonic oscillator.

A lightly damped oscillator driven by a time-dependent forcing function \(F(t)\) is modelled by the differential equation

\(\qquad m\dfrac{d^2y}{dt^2} + c\dfrac{dy}{dt}+k y = F(t) \qquad \qquad \qquad\) (3)

The solution or system response in (3) has two parts:

\( \quad\) (a) A transient solution of the form (2),

\(\quad\) (b) A forced or steady state solution whose form, of course, depends on \(F(t).\)

If \(F(t)\) is sinusoidal such that

\(\qquad F(t)= A\sin (\Omega t+\phi) \, \) where \(\Omega\) and \(\phi\) are constants,

then the steady state solution is fairly readily obtained by standard techniques for solving differential equations. If \(F(t)\) is periodic but non-sinusoidal then Fourier series may be used to obtain the steady state solution. The method is based on the principle of superposition which is actually applicable to any linear (homogeneous) differential equation. (Another engineering application is the series \(LCR\) circuit with an applied periodic voltage.)

The principle of superposition is easily demonstrated:

Let \(y_1(t)\) and \(y_2(t)\) be the steady state solutions of (3) when \(F(t)=F_1(t)\) and \(F(t)=F_2(t)\) respectively. Then

\(\qquad m\dfrac{d^2y_1}{dt^2} + c\dfrac{dy_1}{dt}+k y_1 = F_1(t)\)

\(\qquad m\dfrac{d^2y_2}{dt^2} + c\dfrac{dy_2}{dt}+k y_2 = F_2(t)\)

Simply adding these equations we obtain

\(\qquad m\dfrac{d^2}{dt^2}(y_1+y_2) + c\dfrac{d}{dt}(y_1+y_2) + k(y_1+y_2) = F_1(t) + F_2(t)\)

from which it follows that if \(F(t) = F_1(t) + F_2(t)\) then the system response is the sum \( \, y_1(t) + y_2(t).\) This, in its simplest form, is the principle of superposition. More generally if the forcing function is

\(\qquad F(t )= \sum_{n=1}^NF_n(t)\)

then the response is \(\, \, y(t) = \sum_{n=1}^Ny_n(t)\) where \(y_n(t)\) is the response to the forcing function \(F_n(t).\)

Returning to the specific case where \(F(t)\) is periodic, the solution procedure for the steady state response is as follows:

\( \quad\) Step 1: Obtain the Fourier series of \(F(t).\)

\( \quad\) Step 2: Solve the differential equation (3) for the response \(y_n(t)\) corresponding to the \(n^{\text{ th}}\) harmonic in the Fourier series. (The response \(y_o\) to the constant term, if any, in the Fourier series may have to be obtained separately.)

\( \quad\) Step 3: Superpose the solutions obtained to give the overall steady state motion:

\(\qquad \qquad y(t) = y_0(t)+\sum_{n=1}^Ny_n(t)\)

The procedure can be lengthy but the solution is of great engineering interest because if the frequency of one harmonic in the Fourier series is close to the natural frequency \(\, \sqrt{\dfrac{k}{m}}\) of the undamped system then the response to that harmonic will dominate the solution.

Applying Fourier series to solve a differential equation

The following Task which is quite long will provide useful practice in applying Fourier series to a practical problem. Essentially you should follow Steps 1 to 3 above carefully.

Example

The problem is to find the steady state response \(y(t)\) of a spring/mass/damper system modelled by

\(\qquad m\dfrac{d^2y}{dt^2} + c\dfrac{dy}{dt}+k y = F(t) \qquad \qquad \qquad\) (4)

where \(F(t)\) is the periodic square wave function shown in the diagram.

[Image: Image of a periodic square wave of the function F(t). The maximum amplitude is F_0.]

Solution

Step 1: Obtain the Fourier series of \(F(t)\) noting that it is an odd function.

The calculation is similar to those you have performed earlier in this Workbook.

Since \(F(t)\) is an odd function and has period \(2t_0\) so that \(\, \omega = \dfrac{2\pi}{2t_0} = \dfrac{\pi}{t_0}\), it has Fourier coefficients:

\(\begin{align} b_n &= \frac{2}{t_0}\int_0^{t_o}F_0\sin\left(\frac{n\pi t}{t_0}\right)\,dt \qquad n=1,2,3,\dots \\ &= \left( \frac{2F_0}{t_0} \right) \left( \frac{t_0}{n\pi} \right) \left[ -\cos \frac{n\pi t}{t_0} \right]_0^{t_0}\\ &= \frac{2F_0}{n\pi}(1-\cos n\pi) = \left\{\begin{matrix}\dfrac{4F_0}{n\pi}&n\ \text{ odd}\\ 0&n\ \text{ even} \end{matrix}\right. \end{align}\)

so \( \quad F(t) = \dfrac{4F_0}{\pi}\sum_{n=1}^{\infty}\dfrac{\sin n\omega t}{n} \qquad\)(where the sum is over odd \(n\) only).

Step 2(a): Since each term in the Fourier series is a sine term you must now solve (4) to find the steady state response \(y_n\) to the \(n^{\text{ th}}\) harmonic input: \(\quad F_n(t) = b_n\sin n\omega t \qquad n=1,3,5,\dots\)

From the basic theory of linear differential equations this response has the form

\( \qquad y_n = A_n\cos n\omega t + B_n\sin n\omega t \qquad \qquad\) (5)

where \(A_n\) and \(B_n\) are coefficients to be determined by substituting (5) into (4) with \(F(t)=F_n(t).\) Do this to obtain simultaneous equations for \(A_n\) and \(B_n :\)

We have, differentiating (5),

\( \qquad y_n' = n\omega (-A_n\sin n\omega t+B_n\cos n\omega t)\)

\( \qquad y_n''=(n\omega )^2(-A_n\cos n\omega t-B_n\sin n\omega t) \)

from which, substituting into (4) and collecting terms in \(\cos n\omega t\) and \(\sin n\omega t,\) \[(-m(n\omega)^2A_n + cn\omega B_n+kA_n) \cos n\omega t + (-m(n\omega)^2B_n-cn\omega A_n + kB_n)\sin n\omega t = b_n\sin n\omega t \]

Then, by comparing coefficients of \( \cos n\omega t\) and \(\sin n\omega t\), we obtain the simultaneous equations:

\(\qquad (k-m(n\omega)^2)A_n + c(n\omega)B_n =0 \qquad \qquad \qquad \) (6)

\(\qquad -c(n\omega)A_n + (k-m(n\omega)^2)B_n =b_n \qquad \qquad \qquad \) (7)

Step 2(b): Now solve (6) and (7) to obtain \(A_n\) and \(B_n:\)

\[A_n = -\frac{c\omega_nb_n}{(k-m\omega_n^2)^2 + \omega_n^2c^2} \qquad \qquad (8)\]

\[B_n = \frac{(k-m\omega_n^2)b_n}{(k-m\omega_n^2)^2 + \omega_n^2c^2}\qquad \qquad (9)\]

where we have written \(\omega_n\) for \(n\omega\) as the frequency of the \(n^{\text{ th}}\) harmonic.

It follows that the steady state response \(y_n\) to the \(n^{\text{ th}}\) harmonic of the Fourier series of the forcing function is given by (5). The amplitudes \(A_n\) and \(B_n\) are given by (8) and (9) respectively in terms of the systems parameters \(k,c,m\), the frequency \(\omega_n\) of the harmonic and its amplitude \(b_n\). In practice it is more convenient to represent \(y_n\) in the so-called amplitude form:

\(\qquad y_n = C_n\sin (\omega_nt +\phi_n) \qquad \qquad \qquad \) (10)

where, from (5) and (10),

\( \qquad A_n\cos\omega_nt + B_n\sin\omega_nt = C_n(\cos\phi_n\sin\omega_nt + \sin\phi_n\cos\omega_nt).\)

Hence

\(\qquad C_n\sin \phi_n = A_n \qquad C_n\cos \phi_n=B_n\)

so

\(\qquad \tan \phi_n = \dfrac{A_n}{B_n} = \dfrac{c\omega_n}{(m\omega_n^2-k)^2} \qquad \qquad \qquad \qquad\) (11)

\(\qquad C_n = \sqrt{A_n^2+B_n^2} = \dfrac{b_n}{\sqrt{(m\omega_n^2-k)^2 + \omega_n^2c^2}} \qquad \qquad \) (12)

Step 3: Finally, use the superposition principle, to state the complete steady state response of the system to the periodic square wave forcing function:

\( \qquad y(t) = \sum_{n=1}^{\infty}y_n(t) = \sum_{{n=1}\atop (n \text{ odd})} C_n(\sin\omega_nt + \phi_n) \quad\) where \(C_n\) and \(\phi_n\) are given by (11) and (12).

In practice, since \(\, b_n=\dfrac{4F_0}{n\pi}\) it follows that the amplitude \(C_n\) also decreases as \(\, \frac{1}{n}\). However, if one of the harmonic frequencies say \(\omega_n'\) is close to the natural frequency \(\, \sqrt{\frac{k}{m}}\) of the undamped oscillator then that particular frequency harmonic will dominate in the steady state response. The particular value \(\omega_n'\) will, of course, depend on the values of the system parameters \(k\) and \(m.\)