
In this lecture we will develop the equations that will be implemented in our time-stepping loop. We'll actually implement the solution in code in the next lecture. Let's start by considering the change in force magnitude over a single time-step, between time $t_n$ and $t_{n+1}$, Fig 1.

<Image
	src="/images/courses/fundamentals-of-structural-dynamics-with-python/the-piecewise-exact-method-development/img1.jpg"
	alt="Parameterised linear approximation of the actual loading over one time-step. | EngineeringSkills.com"
	width="400"
	height="320"
	caption="Fig 1. Parameterised linear approximation of the actual loading over one time-step."
/>

If the duration $t_{n+1} - t_n$ is sufficiently small, the error introduced by the linear approximation is negligibly small. We can represent the force as a function of time over the interval $\tau$,

$$
p(\tau) = p_n + \frac{p_{n+1}-p_n}{h}\tau
\tag{1}
$$

The equation of motion that describes the behaviour of the system is therefore,

$$
m\ddot{u} + c\dot{u} + ku = p_n + \frac{p_{n+1}-p_n}{h}\tau
\tag{2}
$$

An analytical solution to this equation of motion can be determined and consists of three parts which will be combined through superposition. Remember, linear analysis is assumed over the duration of a time-step so superposition is valid. The three components of the response are,

- the free vibration component with initial position and velocity $u_n$ and $\dot{u}_n$ respectively.

- the component due to a constant force of magnitude $p_n$.

- a component due to a linearly-varying force,

$$
\frac{p_{n+1}-p_n}{h}\tau
\tag{3}
$$

The solution to the differential equation can be obtained for each of these three cases in turn using the same procedures demonstrated previously (solving to find homogeneous and complementary solutions to the various equations of motion). For brevity we will simply state them here rather than deriving each. The free vibration response is given by,

$$
u_1 = e^{-\xi\omega_n\tau}\left[u_0\cos(\omega_d\tau) + \frac{\dot{u}_0 + \xi\omega_nu_0}{\omega_d}\sin(\omega_d\tau) \right]
\tag{4}
$$

The constant force vibration response is given by,

$$
u_2 = \frac{p_n}{k}\left[1-e^{-\xi\omega_n\tau} \left(\cos(\omega_d\tau) + \frac{\xi\omega_n}{\omega_d} \sin(\omega_d\tau) \right) \right]
\tag{5}
$$

The linearly-varying forced response is given by,

$$
u_3 = \frac{p_{n+1}-p_n}{kh}\left[\tau - \frac{2\xi}{\omega_n} + e^{-\xi\omega_n\tau} \left(\frac{2\xi^2-1}{\omega_d}\sin(\omega_d\tau) + \frac{2\xi}{\omega_n} \cos(\omega_d\tau) \right) \right]
\tag{6}
$$

The system displacement is obtained by adding each of these three components. We would obtain the corresponding velocity by simply differentiating the resulting expression.

In principle, with knowledge of our the system's initial position, velocity and applied force, we could _step_ our way through each time-step to determine the system position and velocity at any time in the future.

### Computational Scheme

In order to implement this concept in a simple time-stepping algorithm, we'll first restructure the equations to make them more _algorithm friendly_. I've used the formulation presented by Humar, (referenced in the resource panel above). You can optionally refer to that text for more on this but everything you need to implement this technique is presented here.

If we differentiate expressions for $u_1$, $u_2$ and $u_3$ we can state the following expressions for the position and velocity at the end of a time-step as,

$$
u_{n+1} = A\:u_n + B\:\dot{u}_{n} +C\:p_n + D\:p_{n+1}
\tag{7}
$$

$$
\dot{u}_{n+1} = A_1\:u_n + B_1\:\dot{u}_n + C_1\:p_n + D_1\:p_{n+1}
\tag{8}
$$

where the constants $A$ to $D_1$ are given by;

$$
A = e^{-\xi\omega_nh}\left[\frac{\xi}{\sqrt{1-\xi^2}}\sin(\omega_dh) + \cos(\omega_dh) \right]
\tag{9}
$$

$$
B=e^{-\xi\omega_nh}\left[\frac{1}{\omega_d}\sin(\omega_dh) \right]
\tag{10}
$$

$$
C=\frac{1}{k}\left[\frac{2\xi}{\omega_nh} + e^{-\xi\omega_nh}\left(\left(\frac{1-2\xi^2}{\omega_dh} - \frac{\xi}{\sqrt{1-\xi^2}} \right)\sin(\omega_dh)-\left(1+\frac{2\xi}{\omega_n h} \right)\cos(\omega_dh)\right) \right]
\tag{11}
$$

$$
D=\frac{1}{k}\left[1-\frac{2\xi}{\omega_nh}+e^{-\xi\omega_nh}\left(\frac{2\xi^2-1}{\omega_dh}\sin(\omega_dh) + \frac{2\xi}{\omega_nh}\cos(\omega_dh) \right) \right]
\tag{12}
$$

$$
A_1=-e^{-\xi\omega_nh}\left[\frac{\omega_n}{\sqrt{1-\xi^2}}\sin(\omega_dh) \right]
\tag{13}
$$

$$
B_1=e^{-\xi\omega_nh}\left[\cos(\omega_dh)-\frac{\xi}{\sqrt{1-\xi^2}}\sin(\omega_dh) \right]
\tag{14}
$$

$$
C_1=\frac{1}{k}\left[-\frac{1}{h} + e^{-\xi\omega_nh}\left(\left(\frac{\omega_n}{\sqrt{1-\xi^2}} + \frac{\xi}{h\sqrt{1-\xi^2}} \right)\sin(\omega_dh) +\frac{1}{h}\cos(\omega_d h) \right) \right]
\tag{15}
$$

$$
D_1=\frac{1}{k}\left[\frac{1}{h} - \frac{e^{-\xi\omega_nh}}{h} \left(\frac{\xi}{\sqrt{1-\xi^2}}\sin(\omega_dh)+\cos(\omega_dh) \right)\right]
\tag{16}
$$

Notice that constants $A$ to $D_1$ are just that, constants. Therefore they can be calculated once at the outset of our solution and then left alone!

Although these expressions appear exceptionally complex. It's important to remember that fundamentally all we have done is implement the same procedures you are familiar with from the analysis of harmonic loading.

It's also worth remembering that the purpose of our work here is to develop a set of equations that is suitable for implementation in a computer algorithm. As such it is not expected that you would manually process these equations to determine the step-by-step response of the structure.

In the next lecture, we'll finally build these expressions into an algorithm for calculating the structural response.


Please log in or enroll to continue

Fundamentals of Engineering Structural Dynamics with Python

Leverage fundamental structural dynamics to build your own flexible numerical solutions in Python.

Please log in or enroll to continue

Fundamentals of Engineering Structural Dynamics with Python

Leverage fundamental structural dynamics to build your own flexible numerical solutions in Python.