
Over the last two lectures we've discussed SDoF motion with respect to a stationary global reference frame. In this lecture, we'll develop the equation of motion that, when solved, yields SDoF motion _relative to the moving support_. This simple shift of reference frame allows us to model the influence of ground motion as an **effective (ground motion) force**, applied to the mass.

Let's start by defining $u_r(t)$ as motion of the SDoF system relative to the moving support structure. $u_s(t)$ will again represent the motion of the base/support structure. The absolute motion of the mass is therefore given by,

$$
u(t) = u_s(t) + u_r(t)
\tag{1}
$$

We'll dispense with the $(t)$ as these quantities are clearly a function of time. Consider again the free-body diagram of our SDoF system from two lectures ago, Fig 1. repeated below.

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/simulating-ground-motion-response/img1.jpg"
	alt=" Free-body diagram of the mass. | EngineeringSkills.com"
	width="450"
	height="227"
	caption="Fig 1. Free-body diagram of the mass."
/>

We can construct the same equation of motion we've seen before,

$$
m\ddot{u} + c\dot{u_r} + ku_r = 0
\tag{2}
$$

where $u_r=u-u_s$ and $\dot{u_r} = \dot{u}-\dot{u}_s$. If we now substitute $u=\ddot{u}_s+\ddot{u}_r$ into the equation of motion and rearrange, we obtain,

$$
m\ddot{u}_r + c\dot{u}_r + ku_r = P_{\text{eff}}
\tag{3}
$$

where,

$$
P_{\text{eff}} = -m\ddot{u}_s
\tag{4}
$$

Here $P_{\text{eff}}$ is the effective inertia force due to base motion. Note also that all motion terms on the left side of the equation are relative to the support structure. This form of the equation of motion is very useful because it allows us to simply scale the support motion by the mass and apply this directly as a force onto the structure.

If the base motion is sinusiodal, then the solution has already been demonstrated in the previous lecture. The real utility of this formulation, as we'll see below, is the ability to apply prescribed support motion, e.g. ground acceleration records. What we have presented thus far can be easily captured graphically in Fig. 2 below.

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/simulating-ground-motion-response/img2.jpg"
	alt="Fig 2. (a) SDoF shear building subject to support acceleration u_s, (b) dynamic influence of support acceleration modelled as an effective inertia force due to ground motion, P_eff (c) equivalent spring-mass-damper system. | EngineeringSkills.com"
	width="550"
	height="324"
/>
<Caption>
  Fig 2. (a) SDoF shear building subject to support acceleration $u_s$, (b) dynamic influence of support acceleration modelled as an effective inertia force due to ground motion, $P_{\text{eff}}$ (c) equivalent spring-mass-damper system.
</Caption>

Once the ground motion has been cast as a force applied directly to the SDoF mass, we can obtain the response of the mass relative to the moving support by any of the means discussed in our previous study of fundamental dynamics. For maximum flexibility, it will generally be convenient to implement a numerical solution.

### The Piecewise Exact Method

In this lecture we will use the piecewise exact method to determine the response of our SDoF system to $P_{\text{eff}}$. We discussed this technique in detail in the previous [Fundamentals of Structural Dynamics course](/course/fundamentals-of-structural-dynamics-with-python/the-piecewise-exact-method-concept). In summary, this is a numerical approach in which we break the applied force into a series of straight-line segments.

We can then obtain the closed-form exact solution for the response of the system subject to a linearly varying loading segment. In this way the solution is literally, piecewise exact. The system position and velocity at the end of one linear load step become the initial conditions for the next linear load step. In this way we step through the force record and evaluate the response piece-by-piece.

The complete solution is approximate rather than exact because the applied loading is approximated as a series of straight line segments. However, if the step-size is small enough this method provides an excellent way of simulating response to arbitrary loading, and it will fit the bill nicely here. Refer back to the prerequisite material for a refresher if required.

Next I'll define a function that implements the piecewise exact method. Again, note that we wrote this function in the fundamentals course so head back there if you haven't covered this already.

The function defined below takes a time and force vector as inputs and iterates through the force vector, computing the position and velocity of the system. Note that the function relies on the dynamic characteristics of the system being defined globally, i.e. they are not explicitly passed in as arguments to the function although you could certainly write the function that way.

```python
#Define a function that takes in a time and force vector
def calculateResponse(tVector, fVector):
    '''
    NOTE: All our system parameters are defined globally and not explicitly passed
    into our function as arguments.
    '''

    nPoints = len(tVector) #Number of data points
    delta_t = tVector[1]-tVector[0]

    #Constants (ref. Dynamics of Structures by Humar)
    A = math.e**(-xi*wn*delta_t)*((xi/(math.sqrt(1-(xi**2))))*
        np.sin(wd*delta_t)+np.cos(wd*delta_t))

    B = math.e**(-xi*wn*delta_t)*((1/wd*np.sin(wd*delta_t)))

    C = (1/k)*(((2*xi)/(wn*delta_t)) + math.e**(-xi*wn*delta_t)*
        ((((1-2*(xi**2))/(wd*delta_t))-(xi/(math.sqrt(1-(xi**2)))))*
        math.sin(wd*delta_t)-(1+((2*xi)/(wn*delta_t)))*np.cos(wd*delta_t)))

    D = (1/k)*(1-((2*xi)/(wn*delta_t))+math.e**(-xi*wn*delta_t)*
        (((2*(xi**2)-1)/(wd*delta_t))*np.sin(wd*delta_t)+
        ((2*xi)/(wn*delta_t))*np.cos(wd*delta_t)))

    A1 = -math.e**(-xi*wn*delta_t)*((wn/(math.sqrt(1-(xi**2))))*
         np.sin(wd*delta_t))

    B1 = math.e**(-xi*wn*delta_t)*(np.cos(wd*delta_t)-((xi)/
         (math.sqrt(1-(xi**2))))*np.sin(wd*delta_t))

    C1 =(1/k)*(-(1/delta_t)+math.e**(-xi*wn*delta_t)*
        ((((wn)/(math.sqrt(1-(xi**2))))+((xi)/(delta_t*
        math.sqrt(1-(xi**2)))))*np.sin(wd*delta_t)+(1/delta_t)*np.cos(wd*delta_t)))

    D1 =(1/k)*((1/delta_t)-(math.e**(-xi*wn*delta_t)/delta_t)*
               ((xi/(math.sqrt(1-(xi**2))))*np.sin(wd*delta_t)+np.cos(wd*delta_t)))

    #INITIAL CONDITIONS
    x0=0
    v0=0

    pos = [] #Initialise an array to hold calculated position
    vel = [] #Initialise an array to hold calculated velocity

    #MAIN LOOP (to 'step through' solution)
    for n, timestep in enumerate(tVector):
        #note that enumerate adds a counter to an iterable (n)
        timestep = round(timestep,3) #The current time

        #Update the force values
        if n<nPoints-1:
            Fn = fVector[n]
            Fnp1 = fVector[n+1]
        else:
            #Pickup when we're at the end of the force vector (no more +1 values)
            Fn = fVector[n]
            Fnp1 = 0

        #Calculate the position and velocity at the end of the timestep
        curPos = (A*x0) + (B*v0) + (C*Fn) + (D*Fnp1)
        curVel = (A1*x0) + (B1*v0)+ (C1*Fn) + (D1*Fnp1)

        #Append (save) the position and velocity into the vectors defined above
        pos.append(curPos)
        vel.append(curVel)

        #Update initial conditions for the next loop iteration
        x0 = curPos
        v0 = curVel

    return pos, vel #Return the calculated position and velocity out of the function
```

### Numerical solution to harmonic excitation

Let's start by solving the equation of motion with a harmonic support excitation. We've already established the closed-form (steady-state) solution to this problem two lectures ago. So this is a great starting point for implementing a numerical solution as we'll have something to compare against.

We first need to establish the effective inertia force,

$$
P_{\text{eff}} = -m \:\ddot{u}_s(t)
\tag{5}
$$

Since $u_s(t) = U_0\sin(\bar{\omega}t)$, we can obtain the acceleration by double differentiating the displacement,

$$
\ddot{u}_s = \frac{\mathrm{d}^2}{\mathrm{d}t^2}u_s = -\bar{\omega}^2U_0\sin(\bar{\omega}t)
\tag{6}
$$

Therefore the effective force is given by,

$$
P_{\text{eff}} = m\bar{\omega}^2U_0\sin{\bar{\omega}t}
\tag{7}
$$

$$
P_{\text{eff}} = m\bar{\omega}^2\:u_s(t)
\tag{8}
$$

Now let's define some parameters and a time vector.

```python
# System parameters
m = 1000 #(kg) Mass
xi = 0.05 #Damping ratio
fn = 1 #(Hz) Natural frequency
wn = 2*math.pi*fn #(rad/s) Angular natural frequency
wd = wn*math.sqrt(1-xi**2) #Damped natural frequency
k = m*wn**2 #(N/m) Stiffness

# Time vector
tmax = 20 #(sec) The max time
delta_t = 0.01 #(sec) The timestep
nPoints = tmax/delta_t #Number of equally spaced data points
t = np.linspace(0,tmax, int(nPoints)) # Time vector

# Support motion parameters
U0 = 0.1 #(m) Support motion amplitude
fb = 1 #(Hz) Support motion frequency
wb = 2*math.pi*fb #(rad/s) Angular support motion frequency
us = U0*np.sin(wb*t)
```

Next we can define the effective inertia force and pass it to our `calculateResponse()` function.

```python
Peff = m*wb**2*us #(N) The effective inertia force
position, velocity = calculateResponse(t, Peff) #Call the function
```

Now that we have the numerical solution calculated, let's generate the closed-form solution so we have something to compare it to. We can re-implement the solution we coded two lectures ago. Remember, this closed form solution gives us the absolute position of the mass. We will need to subtract the support motion, $u_s$ in order directly compare it with our numerical solution.

```python
#Generate the closed form response for comparison
beta = wb/wn #Frequency ratio
D = U0*k*math.sqrt(1+(2*xi*beta)**2)
gamma = math.atan(2*xi*beta)

uConst = (D/k)*(1/((1-beta**2)**2+(2*xi*beta)**2))
uVar = (1-beta**2)*np.sin(wb*t+gamma) - 2*xi*beta*np.cos(wb*t+gamma)
u = uConst * uVar #(m) Absolute displacement of the SDoF mass
```

Now we can verify the accuracy of our numerical solution by visually comparing it to the closed-form solution we derived previously.

```python
#Plotting
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,u-us, '-', label='Relative steady-state response')
axes.plot(t,position, label='Relative numerical response')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Support and mass motion')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/multi-degree-of-freedom-dynamics-modal-analysis-and-seismic-response-simulation-in-python/simulating-ground-motion-response/img3.jpg"
	alt="Plot of analytical steady-state and numerical response. | EngineeringSkills.com"
	width="1000"
	height="525"
	caption="Fig 3. Plot of analytical steady-state and numerical response."
/>

We can see that once the transient behaviour of the system dies away, both solutions are in close agreement. This gives us confidence to deploy our numerical solution to something a little more practical.


Please log in or enroll to continue

Please log in or enroll to continue