Equation of motion

The following is a discussion about solving non-homogeneous second order partial differential equation of the following form

  • The above equation of motion defines the single degree of freedom system subjected to arbitrary force.
  • The arbitrary force can be periodic or non-periodic function w.r.t time.
  • The variable u(t) is the displacement response to the excitation force.


The following are the numerical integration techniques discussed below

      • Central difference method
      • Linear acceleration method
      • Neuman β method
      • Hilber Hughes Taylor (HHT) – α method

Central difference method

Python code for central difference method

def numericalintegration_CDM(initial_disp,initial_velocity,mass,damping,stiffness,natural_period,t_list0, force_inpt):
    # force input
    t_list = t_list0.tolist() # convert to list from np.linespace (time list)

    displacement = []
    velocity = []
    acceleration = []
    #selection of delta t
    # the force delta_t is 0.001
    delta_t =(0.0001*natural_period)/np.pi # delta_t has to be less than Tn/pi stability condition


    #initial calculation
    inital_accl = (force_inpt[0] - (damping*initial_velocity) - (stiffness*initial_disp))/mass
    disp_minus1 = initial_disp - (delta_t*initial_velocity )+ ((delta_t**2)*0.5*inital_accl)
    disp_0 = initial_disp
    k_hat = (mass/(delta_t**2)) + (damping/(2*delta_t))
    a = (mass/(delta_t**2)) - (damping/(2*delta_t))
    b = stiffness - ((2*mass)/(delta_t**2))
    u_static = max(force_inpt)/stiffness # static displacement p0/k
    vel_0 = initial_velocity # initial velocity
    accl_0 = inital_accl # initial acceleration

    time_tdtn = [] #time t/Tn

    for t0 in t_list:
        t_index0 = t_list.index(t0)  # index of time
        p_minus1 = force_inpt[t_index0]  # force at time t_index

        displacement.append(disp_0/u_static)
        time_tdtn.append(t0/natural_period)
        velocity.append(vel_0)
        acceleration.append(accl_0)
        if t_index0 == (len(t_list)-1):
            break # check the end of time
        t1 = t_list[t_index0 +1]
        p0 = force_inpt[t_index0+1]
        for t_x in np.arange(t0,t1,delta_t):
            t_int = (t_x - t0) / (t1-t0)
            p_i = (p_minus1 * (1-t_int) )+ (t_int * p0)
            p_hat = p_i - (a*disp_minus1) - (b*disp_0)
            accl_0 = (disp_minus1 - (2*disp_0) + (p_hat/k_hat))/(delta_t**2)
            disp_minus1 = disp_0
            disp_0 = p_hat/k_hat
            vel_0 = (disp_0 - disp_minus1)/(2*delta_t)
    return displacement,time_tdtn

Solution 1:

Dynamic response of un-damped SDF system to half-cycle sine pulse force

For the above solution, case td/Tn = 1.5,2.5,... the mass stays still after the force ends (Page 145: Dynamic of Structures - Anil K Chopra. Contrary to the text book solution my displacement response didn't become zero for td/Tn = 1.5,2.5,... at the end of pulse. The response is very low but not zero after the pulse ends unlike the theoretical solution.

newmark's - Beta method

•Newmark method is named after Nathan M. Newmark, former professor of civil engineering at the University of Illinois at Urbana–Champaign, who developed it in 1959 for use in structural dynamics.

•Interpolation equations (Newmark 1959)

•Parameter γ = ½

β = 0 central difference method

β = ¼ constant acceleration method, equivalent to trapezoidal rule

β = 1/6 linear acceleration method

β = 1/12 Leslie Fox and Goodwin proceddure

Python code for linear acceleration method

def numericalintegration_linearacceleration2(initial_disp,initial_velocity,mass,damping,stiffness,natural_period,t_list0, force_inpt,force_tstep):
    # force input
    t_list = t_list0.tolist()  # convert to list from np.linespace (time list)

    displacement = []
    velocity = []
    acceleration = []
    # selection of delta t
    # the force delta_t is 0.001
    delta_t = (0.01*force_tstep * natural_period) / np.pi  # delta_t has to be less than Tn/pi (which is the requirement for stability)

    gamma = 1/2
    beta = 1/6

    # initial calculation
    inital_accl = (force_inpt[0] - (damping * initial_velocity) - (stiffness * initial_disp)) / mass
    disp_0 = initial_disp  # initial displacement
    a1 = (1/ (beta*(delta_t**2)))*mass  + (gamma/(beta*delta_t))*damping
    a2 = (1/(beta*delta_t))*mass + ((gamma/beta)-1)*damping
    a3 = ((1/(2*beta))-1)*mass  + (delta_t*((gamma/(2*beta))-1))*damping

    k_hat = stiffness + a1
    u_static = max(force_inpt) / stiffness  # static displacement p0/k
    vel_0 = initial_velocity  # initial velocity
    accl_0 = inital_accl  # initial acceleration

    time_tdtn = []  # time t/Tn

    for t0 in t_list:
        t_index0 = t_list.index(t0)  # index of time
        p_index0 = force_inpt[t_index0]  # force at time t_index

        displacement.append(disp_0 / u_static)
        time_tdtn.append(t0 / natural_period)
        velocity.append(vel_0)
        acceleration.append(accl_0)
        if t_index0 == (len(t_list) - 1):
            break  # check the end of time
        t1 = t_list[t_index0 + 1]
        p1 = force_inpt[t_index0 + 1]
        disp_minus1 = disp_0
        vel_minus1 = vel_0
        accl_minus1 = accl_0

        #p_minus1 = p_index0 # delta force

        for t_x in np.arange(t0, t1, delta_t):
            t_int = (t_x - t0) / (t1 - t0)
            p_i = (p_index0 * (1 - t_int)) + (t_int * p1)

            p_hat = p_i + (a1*disp_0) + (a2*vel_0) + (a3*accl_0)

            disp_0 = p_hat/k_hat
            vel_0 = ((gamma/(beta*delta_t))*(disp_0-disp_minus1))+((1-(gamma/beta))*vel_minus1 )+ (delta_t*(1-(gamma/(2*beta)))*accl_minus1)
            accl_0 = ((1/(beta*(delta_t**2)))*(disp_0-disp_minus1))-((1/(beta*delta_t))*vel_minus1) - (((1/(2*beta))-1)*accl_minus1)
            disp_minus1 = disp_0
            vel_minus1 = vel_0
            accl_minus1 = accl_0
    return displacement

Solution 2:

Dynamic response of un-damped SDF system to half-cycle square pulse force

Below is the heat map showing the dynamic response of un-damped SDF system to half-sine pulse force of time 1 sec duration (td = 1sec) for various td/Tn ratio and for 3 sec total duration of response.

When looking only at the absolute value of the heat map above. We can clearly see the no response zone @ td/Tn = 1.5,2.5,3.5 etc after the impulse ended.

Shock spectrum

Plotting the heat map in 3d shows the shock spectrum for various td/Tn

Hilber Hughes Taylor (HHT) – alpha method

in progress