
In this lecture, we'll develop a simple algorithm to implement the piecewise-exact method. Let's start by defining a force-time history to work with. For convenience we'll assume that the force applied to our structure is a simple sine wave; but remember, the process we're about to implement will work for **any** force-time history - not just simple sine waves.

```python
#Define a time vector
tmax = 30 #(sec) The max time
delta_t = 0.01 #(sec) The time-step
nPoints = tmax/delta_t #Number of equally spaced data points
t = np.linspace(0,tmax, int(nPoints)) # Time vector

# Define a force-time history to apply to our structure
P = 750 # (N) Force magnitude
f = 1 # (Hz) Forcing frequency
omega = 2*pi*f #(rads/sec) Angular forcing frequency
F_sine = P*np.sin(omega*t) #Force vector
```

Next we can plot the force; we'll refer to this force as an _experimental_ force just to emphasise that it needn't necessarily be a sine wave.

```python
#Plot the 'experimental' force and sampling points
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,F_sine,label='Experimental force')
axes.plot(t,F_sine,'*',label='Sample points')
axes.set_xlim([0,1])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Force (N)')
axes.set_title('Force-time history')
axes.grid()
axes.legend()
plt.show()
```

<Image
	src="/images/courses/fundamentals-of-structural-dynamics-with-python/implementing-the-piecewise-exact-method/img1.jpg"
	alt=" Force-time history. | EngineeringSkills.com"
	width="1000"
	height="528"
	caption="Fig 1. Force-time history."
/>

We've also plotted the individual data points that make up the sine wave. These will serve as our force sample points within the algorithm. Our task now is to build a function that calculates the response of a SDoF system subject to the force indicated by the orange data points above.

We'll define a function that accepts a time vector and force vector as input arguments and returns the structural response.

```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

    #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 #Return the calculated position out of the function
```

Within the function above, we are:

- defining the constants $A$ to $D_1$ identified in the previous lecture.

- setting the initial position and velocity of the system to zero, although these could also be passed into the function as arguments allowing us to handle non-zero initial conditions.

- initialising two data structures or containers to hold the calculated position and velocity values.

- setting up a `for` loop to cycle through each time-step.

- within the `for` loop, on each iteration we are calculating the current value of the force, the position and velocity, saving the values calculated and updating the initial conditions for the next iteration of the loop.

- finally, once the loop is complete, we return the array of calculated positions from the function.

Logically, it's all quite straightforward, even if some of the equations along the way were quite long and unwieldy. Now we can call the function and pass it our time and force data. But first we need to redefine the SDoF system parameters as our function relies on these being globally defined.

```python
# SYSTEM CONSTANTS
m = 1000 #(kg) Mass
k = 20000 #(N/m) Stiffness
xi = 0.05 # Damping ratio
omega_n = math.sqrt(k/m) #(rads/sec) Undamped natural frequency
omega_d = omega_n*math.sqrt(1-xi**2) #Damped natural frequency
f = 1 # (Hz) Forcing frequency
beta = omega/omega_n # Frequency ratio
```

```python
position = calculateResponse(t, F_sine) #Calling our function
```

And plotting the response,

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,1.5])
axes.plot(t,position,'-')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Displacement-time history')
axes.grid()
plt.show()
```

<Image
	src="/images/courses/fundamentals-of-structural-dynamics-with-python/implementing-the-piecewise-exact-method/img2.jpg"
	alt="Displacement-time history. | EngineeringSkills.com"
	width="1000"
	height="527"
	caption="Fig 2. Displacement-time history."
/>

**Remember...**

- This approach will work for **any** force vector - we just happen to define a sinusoid here because it was convenient for demonstration purposes.

- This **numerical solution** is well suited to evaluating the impact of experimentally observed forces (represented as a list of sampled force values).

- This method is also well suited to modelling dynamic behaviour of structures that exhibit non-linear behaviour because we can alter our algorithm to change dynamic properties between iterations.

- We can also alter the applied force between iterations even making it dependent on the structural response in the previous iteration - modelling feedback between force and response.

For completeness, let's compare our numerical solution to the closed-form solution for harmonic excitation we derived previously,

$$
u(t) = e^{-\xi \omega_nt}\:\left[A\sin(\omega_d t) +B\cos(\omega_d t)\right] + \frac{P_0}{k}\left[\frac{1}{(1-\beta^2)^2+(2\xi\beta)^2}\right]\left[(1-\beta^2)\sin(\omega t)-2\xi\beta\cos(\omega t)\right]
\tag{1}
$$

Following the same procedure we worked through in Worked Example #3 to calculate the response:

```python
#Constants for calculation of closed-form solution (refer to Worked Example #3)
O = (P/k)*(1/((1-beta**2)**2 + (2*xi*beta)**2)) #Capital Omega defined above
C1 = O*(1-beta**2) #Defined above
C2 = -O*2*xi*beta #Defined above
omega_d = omega_n*math.sqrt(1-xi**2) #Damped natural frequency
A = (-1/omega_d)*(C2*xi*omega_n + C1*omega) #Constant of integration
B = -C2 #Constant of integration

ut = math.e**(-xi*omega_n*t)*(A*np.sin(omega_d*t) + B*np.cos(omega_d*t)) #Transient response
us = O*((1-beta**2)*np.sin(omega*t) - 2*xi*beta*np.cos(omega*t)) #Steady-state response
```

Now we can compare the solutions generated by our numerical and closed-form solutions.

```python
fig = plt.figure()
axes = fig.add_axes([0.1,0.1,2,2])
axes.plot(t,(ut+us), label='Closed-form solution')
axes.plot(t,position, label= 'Numerical solution')
axes.set_xlim([0,tmax])
axes.set_xlabel('Time (sec)')
axes.set_ylabel('Displacement (m)')
axes.set_title('Displacement-time history')
axes.grid()
axes.legend(loc='lower right')
plt.show()
```

<Image
	src="/images/courses/fundamentals-of-structural-dynamics-with-python/implementing-the-piecewise-exact-method/img3.jpg"
	alt="Displacement-time history comparison between numerical and analytical solution. | EngineeringSkills.com"
	width="1000"
	height="681"
	caption="Fig 3. Displacement-time history comparison between numerical and analytical solution."
/>

As we might expect there is pretty much perfect agreement between both responses. In the next lecture we'll consider some of the factors that influence the accuracy of our numerical response.


Please log in or enroll to continue

Please log in or enroll to continue