Universidade Federal do Rio Grande do Sul (UFRGS)¶

Programa de Pós-Graduação em Engenharia Civil (PPGEC)

PEC00025: Introduction to Vibration Theory¶

Class 06 - Numerical integration: finite differences and Duhamel techniques¶

1. Introduction
2. The finite differences technique
2.1. Formulation
2.2. Example
3. The Duhamel's integral technique
3.1. Formulation
3.2. Example
4. Example: excitation recorded as file
5. Finite differences for nonlinear equations
6. Assignment

---

Prof. Marcelo M. Rocha, Dr.techn. (ORCID)
Porto Alegre, RS, Brazil

1. Introduction ¶

Besides the basic conceptual solutions presented in the previous class, in practical cases where the external force $F(t)$ is arbitrarily defined, the dynamic equilibrium equation must be numerically solved.

In this class we present two of the most important integration techniques:

• the central finite differences scheme and
• the Duhamel integral method.

There is also another useful technique that will be presented after we introduce the Fourier transform and the frequency domain analysis of linear systems, in the next classes. For the moment we can state that time domain techniques are mainly used for transient loads, while frequency domain techniques are maily used for periodic or stationary random loads.

2. The finite differences technique ¶

2.1. Formulation ¶

Although the central finite differences technique may be used to solve the dynamic equilibrium equation of linear systems, its most important application is in the solution of nonlinear systems. For linear systems it allows the integration scheme to be explicit, what means that the system state (displacement) in future can be explicitly stated as a function of system state in the past. To understand this concept, let us discretize the time domain as $t_i = i \Delta t$, where $\Delta t$ is the so-called time step. This is depicted in the figure below:

Consequently, both the sistem response, $u(t)$, as well as the external load, $F(t)$, can be now expressed at the discrete time instants, $t_i$, as:

$$\begin{align*} u_i &= u(t_i) \\ F_i &= F(t_i) \end{align*}$$

From this definitions, the system velocity, $\dot{u}(t)$, and acceleration, $\ddot{u}(t)$, can be approximated through a central finite differences scheme as:

$$\begin{align*} \dot{ u}_i &= \frac{1}{2\Delta t} \left( u_{i+1} - u_{i-1} \right) \\ \ddot{u}_i &= \frac{1}{ \Delta t^2} \left( u_{i+1} - 2u_ {i} + u_{i-1} \right) \end{align*}$$

Considering that the dynamic equilibrium equation holds also for a generic time instant $t_i$:

$$\ddot{u}_i + 2 \zeta \omega_{\rm n} \dot{u}_i + \omega_{\rm n}^2 u_i = \frac{F_i}{m}$$

the discrete kinematic parameters $\dot{ u}_i$ and $\ddot{ u}_i$ can be replaced to give:

$$\frac{1}{ \Delta t^2} \left( u_{i+1} - 2u_ {i} + u_{i-1} \right) + \frac{2 \zeta \omega_{\rm n}}{2 \Delta t} \left( u_{i+1} - u_{i-1} \right) + \omega_{\rm n}^2 u_i = \frac{F_i}{m}$$

By isolating the future system displacement, $u_{i+1}$, it results:

$$u_{i+i} = \frac{1}{\beta_1} \left( \frac{F_i}{m} + \beta_2 u_{i-1} - \beta_3 u_{i} \right)$$

with:

$$\begin{align*} \beta_1 &= \frac{1}{\Delta t} \left( \zeta\omega_{\rm n} + \frac{1}{\Delta t}\right) \\ \beta_2 &= \frac{1}{\Delta t} \left( \zeta\omega_{\rm n} - \frac{1}{\Delta t}\right) \\ \beta_3 &= \frac{1}{\Delta t} \left( \Delta t \omega_{\rm n}^2 - \frac{2}{\Delta t}\right) \end{align*}$$

which is an explicit expression of the system response at time $t_{i+1}$ as a function of system responses at times $t_{i-1}$ and $t_i$. The initial position, $u_0$, and initial velocity, $v_0 = \dot{u}_0$, must be provided for the calculation of system responses at instants $t_0$ and $t_1$. The position at $t_0$ is the provided value $u_0$ itself, while the position at $t_1$ can be calculated with:

$$u_1 = u_0 + v_0 \Delta t + \frac{F_0 \Delta t^2}{2}$$

The finite differences technique presents a severe restriction for the (arbitrarily chosen) time step, $\Delta t$:

$$\Delta t \leq \frac{T_{\rm n}}{4}$$

ou

$$f_{\rm S} \geq \frac{4}{T_{\rm n}}$$

Whenever this restriction is not respected, the integration scheme will diverge and become unstable, leading to increasingly large response. Furthermore, as one tries to increase $\Delta t$ to reduce computational time, the technique loose accuracy quite fast. For these reasons, Duhamel's technique presented in the next section is usually preferable for solving linear systems.

2.2. Example of application with MRPy¶

The formulation above is implemented as a method for the MRPy class, as exemplified in the following example with a unit step loading.

Firstly, it is necessary to define the mechanical properties of the sdof system. Let us assume that:

Once the system properties are specified, the unit step loading can be created as a MRPy instance by calling the appropriate constructor:

Observe that the number of time steps N and the total duration Td will define the sampling frequency fs and the time step will be $\Delta t = 1/f_{\rm s}$. Smaller values of N will then worse the integration accuracy for a given Td.

The solution by finite differences is then available through the sdof_fdiff method:

It can be clearly observed that, after the step onset, the system will oscillate with decreasing amplitude around the static response. This static displacement can be calculated directly from system stiffness:

Let us see now what happens whenever the time step is too long in comparison to the system natural period of vibration.

Since we choose $\Delta t < T_{\rm n}$, the solution has diverged!

To make full use of the Jupyter notebook concept, try now changing the fraction of $T_{\rm n}$ used to define the time step in the example above. What is the worst acceptable accuracy?

3. The Duhamel's integral technique ¶

3.1. Formulation ¶

The numerical solution by the Duhamel's integral technique is restricted to linear systems. It relies on the superposition of system responses to a sequence on impulses $F(\tau) \, d\tau$. It has been seen that the general solution of the equilibrium equation of a sdof system subjected to a general load $F(t)$ is given by the convolution of this load with the impulse response:

$$u(t) = u_0(t) + \frac{1}{m \omega_{\rm D}} \int_0^t \exp \left[ -\zeta\omega_{\rm n}(t - \tau) \right] \; \sin \omega_{\rm D} (t - \tau) \; F(\tau) \; d\tau$$

where $u_0(t)$ is the system response to the initial conditions. Considering now the trigonometric identity:

$$\sin \omega_{\rm D} (t - \tau) = \sin \omega_{\rm D}t \; \sin \omega_{\rm D}\tau - \cos \omega_{\rm D}t \; \cos \omega_{\rm D}\tau$$

and disregarding for a moment the initial conditions leads to:

$$u(t) = \frac{1}{m \omega_{\rm D}} \; \left[ \frac{A(t) \sin \omega_{\rm D}t - B(t) \cos \omega_{\rm D}t} {\exp ( -\zeta\omega_{\rm n}t )} \right]$$

with:

$$\begin{align*} A(t) &= \int_0^t \exp ( -\zeta\omega_{\rm n}\tau ) \; \cos \omega_{\rm D} \tau \; F(\tau) \; d\tau \\ B(t) &= \int_0^t \exp ( -\zeta\omega_{\rm n}\tau ) \; \sin \omega_{\rm D} \tau \; F(\tau) \; d\tau \end{align*}$$

The reformulations above suggest that a recursive scheme can be used to spare computational time spent in the calculations of trigonometric functions. By using the discretized time $\tau_i = i \Delta \tau$ gives the following discrete evaluations:

$$\begin{align*} e_i &= \exp ( -\zeta\omega_{\rm n} \tau_i ) \\ s_i &= \sin \omega_{\rm D} \tau_i \\ c_i &= \cos \omega_{\rm D} \tau_i \\ F_i &= F(\tau_i) \end{align*}$$

The functions $A_i = A(t_i)$ and $B_i = B(t_i)$ can be readily calculated as cumulative summations:

$$\begin{align*} A_i &= \Delta \tau \; \sum_{j=0}^i e_j c_j F_j \\ B_i &= \Delta \tau \; \sum_{j=0}^i e_j s_j F_j \end{align*}$$

and the solution is finally obtained as:

$$u_i = u_{0i} + \frac{1}{m \omega_{\rm D} } \left( \frac{A_i s_i - B_i c_i}{e_i} \right)$$

with the response to the initial conditions (displacement $u_0$ and velocity $v_0$) from:

$$u_{0i} = \frac{1}{e_i} \left[ u_0 c_i + \left( \frac{v_0 + u_0 \zeta \omega_{\rm n}}{\omega_{\rm D}} \right) s_i \right]$$

This is a quite fast computational scheme, for Python can build de time vector and evaluate de functions through built-in numpy commands, including the cumulative summation. Furthermore, the Duhamel technique does not become unstable no matter how coarse the chosen time step is, although its accuracy may not be ensured in this case.

3.2. Example of application with MRPy¶

The formulation above is also implemented as a method for the MRPy class, as exemplified in the following. The same system properties and unit step loading from previous example are used. The class method is:

On the contrary of the diverged solution by finite differences, even by choosing $\Delta t < T_{\rm n}$ does not cause the solution to diverge.

4. Example: excitation recorded as file ¶

The MRPy class provides a constructor that reads a time series from file. There are some formatting options and others can be included as needed. In this example we will read a file recorded with the cell phone app iNVH by Bosch. The file has a .csv format with four columns, the first being the sampling time and the next three being the accelerations along the three measurement axes ($x$, $y$ and $z$).

Reading the file is a straighforward calling to the from_file constructor:

This data was obtained at the center of a metallic beam subjected to multiple impacts. We are interested in the acceleration along the $z$ axis, so it can be isolated from the complete set with:

Now this acceleration will be applied as a basis excitation to a sdof system with the same properties as defined in the previous sections. The solutions by finite differences and Duhamel are compared:

It seems that both techniques yield almos the same results, so let us take a look at the error:

This means that, for the provided sampling rate, the integration error is very small.

5. Assignments ¶

1. Utilizar as propriedades (frequência e amortecimento) do trabalho anterior.
2. Utilizar como entrada o registro de aceleração sísmica fornecido abaixo.
3. Calcular a resposta com MRPy.sdof_FDiff() e MRPy.sdof_Duhamel().
4. Relatório com descrição do objeto, gráficos e resultados.