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