
In this lecture we're going to explore one method of solving the coupled second order system of equations identified in the previous lecture. As the equations are coupled, the most practical approach to solving them will be a numerical integration. With some effort, it may be possible to develop analytical solutions in the case of simple external forcing functions. However, when the external force is described by a vector of numbers (e.g. $P_{\text{eff}}$, derived from a ground acceleration record), numerical integration is the only approach. This is why we explore it here.

There are various numerical integration techniques and indeed many lectures could be written on the topic. We're going to focus on implementing one of the most commonly deployed, the **central difference method**. We'll explain the technique in this lecture and implement it in the next.

### The central difference method

The central difference method is a time-stepping numerical integration technique in which we evaluate the response at time $t$ by computing the velocity and acceleration in terms of the response as some previous time. In this way we step through the time dimension, calculating values of response and advancing our solution by using the response values calculated during previous time steps. The means that we must have knowledge of the system initial conditions to kick the process off.

At first glance it can be a confusing concept so it's best explained by example. For the purposes of explanation, we'll revert back to the equation for a SDoF system,

$$
m\ddot{u} + c\dot{u} + k u = p(t)
\tag{1}
$$

The first task is to determine an expression for $\ddot{u}$ and $\dot{u}$ in terms of $u$. We'll work on the assumption that time is divided into equally spaced intervals or time-steps and that we can determine values for $\ddot{u}$, $\dot{u}$ and $u$ at these discrete time points.

To determine an expression for $\dot{u}$, consider the plot of $u$ versus $t$ and the three data points identified, $u_{n-1}$, $u_n$ and $u_{n+1}$, Fig 1.

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/direct-integration-of-the-coupled-eom/img1.jpg"
	alt="Plot of u verus t identifying \dot{u}. | EngineeringSkills.com"
	width="300"
	height="282"
/>
<Caption type="Fig 1. Plot of $u$ verus $t$ identifying">
	Fig 1. Plot of $u$ verus $t$ identifying $\dot{u}$.
</Caption>

The _central_ velocity, $\dot{u}_n$ is the rate of change of displacement. If we have a sufficiently small time-step such that the displacement can be assumed to vary linearly we can use displacements $u_{n-1}$ and $u_{n+1}$ to calculate the rate of change (slope) of the displacement. This gives the following central difference formula for velocity in terms of displacement,

$$
\boxed{\dot{u}_n = \frac{1}{2h}(u_{n+1}-u_{n-1})}
\tag{2}
$$

where, $h$ is the time interval between data points. We can determine an expression for $\ddot{u}_n$ using the same strategy. If we define the velocity $\dot{u}_a$ and $\dot{u}_b$ as shown in Fig 2. as,

$$
\dot{u}_a = \frac{u_n - u_{n-1}}{h}
\tag{3}
$$

$$
\dot{u}_b = \frac{u_{n+1} - u_n}{h}
\tag{4}
$$

then the acceleration $\ddot{u}_n$ is given by,

$$
\ddot{u}_n = \frac{\dot{u}_b - \dot{u}_a}{h}
\tag{5}
$$

$$
\boxed{\ddot{u}_n = \frac{1}{h^2}(u_{n+1}-2u_n + u_{n-1})}
\tag{6}
$$

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/direct-integration-of-the-coupled-eom/img2.jpg"
	alt=" Plot of u verus t identifying \ddot{u}. | EngineeringSkills.com"
	width="300"
	height="281"
/>
<Caption type="Plot of $u$ verus $t$ identifying $\ddot{u}$">
	Fig 2. Plot of $u$ verus $t$ identifying $\ddot{u}$.
</Caption>

We can now substitute our expressions for $\ddot{u}_n$ and $\dot{u}_n$ into the equation of motion at time point $n$,

$$
m\left[ \frac{1}{h^2}(u_{n+1}-2u_n + u_{n-1})\right] + c\left[\frac{1}{2h}(u_{n+1}-u_{n-1}) \right] + k u_n = p_n
\tag{7}
$$

Because we have constructed the equation of motion at time point $n$, rather than $n+1$, the central difference method is referred to as an explicit integration method. If we wrote the equation of motion at $n+1$, we would have an implicit integration scheme. We note this here as you're likely to come across this classification of numerical integration methods in further study.

After rearranging we get,

$$
\boxed{\left(\frac{m}{h^2} + \frac{c}{2h} \right) u_{n+1} = p_n + \left(-k+\frac{2m}{h^2} \right)u_n + \left(\frac{c}{2h} - \frac{m}{h^2} \right)u_{n-1}}
\tag{8}
$$

Using this equation we can obtain $u_{n+1}$ with knowledge of the preceding displacement values. Therefore we must know what the initial displacement value is. Slightly more troubling, we also need to know what $u_{n-1}$ is.

To obtain an expression for $u_{n-1}$, we can evaluate the expressions for $\ddot{u}_n$ and $\dot{u}_n$ at $t=0$ and solve for $u_{n-1}$,

$$
\boxed{u_{n-1} = u_n + \frac{h^2}{2}\ddot{u}_n - h\dot{u}_n}
\tag{9}
$$

Thus we need the initial position, velocity and acceleration in order to initiate our time-stepping solution. The preceding equations can easily be implemented within a `for` loop to determine the dynamic response of a SDoF structure.

### MDoF systems

Our discussion so far has explained how to apply the central difference method to SDoF systems. Fortunately, application to MDoF systems is identical with the exception that our scaler equations identified above now become matrix equations,

$$
\dot{\mathbf{u}}_n = \frac{1}{2h}(\mathbf{u}_{n+1}-\mathbf{u}_{n-1})
\tag{10}
$$

$$
\ddot{\mathbf{u}}_n = \frac{1}{h^2}(\mathbf{u}_{n+1}-2\mathbf{u}_n + \mathbf{u}_{n-1})
\tag{11}
$$

$$
\left(\frac{1}{h^2}\mathbf{M} + \frac{1}{2h}\mathbf{C} \right) \mathbf{u}_{n+1} = \mathbf{p}_n + \left(-\mathbf{K}+\frac{2}{h^2}\mathbf{M} \right)\mathbf{u}_n + \left(\frac{1}{2h}\mathbf{C} - \frac{1}{h^2}\mathbf{M} \right)\mathbf{u}_{n-1}
\tag{12}
$$

$$
\mathbf{u}_{n-1} = \mathbf{u}_n + \frac{h^2}{2}\ddot{\mathbf{u}}_n - h\dot{\mathbf{u}}_n
\tag{13}
$$

Note that for brevity we've temporarily switched to alternative bold font matrix notation where $\begin{bmatrix}M\end{bmatrix} = \mathbf{M}$.

### Stability

Like any numerical technique, we need to be mindful of the stability of the scheme. As a general rule we need to ensure that the time-step is sufficiently small. This raises the question of how small?

If we imagine a SDoF system we would like the time-step to be small when compared to the natural period of the system. So a reasonable approach would be to have a time-step no larger than $T_n/5$ although you will also find reference to $T_n/\pi$ in the literature.

As we'll see in the next section, MDoF systems actually possess ascending natural frequencies. There is much to unpack here, but for now, we'll say that the time-step should be less than $T_n/5$ for the highest of the natural frequencies.

As we're solving the equations of motion through direct integration, the system natural frequencies remain obscured so we can't identify a time-step that meets our criteria. However a time-step of $0.005$ seconds is likely to be more than satisfactory for most civil engineering structures. We can always investigate the sensitivity of the calculated response to variations in time-step if we are unsure of its suitability.

In the next lecture we'll implement the central difference scheme and finally solve for the dynamic behaviour of our three-storey building.


Please log in or enroll to continue