Hi Guys
Please refer to the below code. I am not going to break down the mathematics, but this system essentially presents a cart with a vertical pendulum that has a rotational stiffness, the cart also has a spring attached to it however this is not always active.
In the “statespace!” function there are a few “if” statements which are used to develop different behavior. The system presents the cart travelling without any spring attached to it for 2.5 seconds (this is the T_to_D variable). When the total displacement of the cart is subtracted by D (the distance traveled when the spring is activated) is smaller than 1e-2 the kb variable which is the stiffness of the spring is set to 0 again. This part oft he code does not seem to be working well.
Refer to the attached plot, this is plot_5. This plots the linear displacement however x[3] has D taken away from it so at the 0 displacement is where the cart hits the spring. The system should be decoupled when “y1” goes past 0 again, which is not occurring.
What could the possible issue be?
using DifferentialEquations
using Plots
global decoupling = []
function statespace!(dx,x,p,t)
l,m,g,I,kc,kb,C1,C2,D,Acc,T_to_D = p
Kb = kb
if t > T_to_D # time at and after impact
kb = Kb
if (x[3]-D) < 1e-2 && t > T_to_D # original position of spring and decoupling of system
kb = 0
push!(decoupling, t) # Append the current time to the decoupling array
end
A = x[2]^2*l*sin(x[1])+(-kb*x[3]-C1*x[4])/m
B = -kc*x[1]-m*g*l*sin(x[1])-C2*x[2]
C = (1-(m*l^2*cos(x[1])^2)/(l^2*m+I))
theta_dotdot = ((-m*A*l*cos(x[1]+B)/(l^2*m+I)))*1/C
x_dotdot = x[2]^2*l*sin(x[1])-theta_dotdot*l*cos(x[1])+(-kb*x[3]-C1*x[4])/m
elseif t < T_to_D # acceleration with cart only
kb = 0
A = x[2]^2*l*sin(x[1])+(-kb*x[3]-C1*x[4])/m
B = -kc*x[1]-m*g*l*sin(x[1])-C2*x[2]
C = (1-(m*l^2*cos(x[1])^2)/(l^2*m+I))
x_dotdot = Acc
theta_dotdot = ((-m*A*l*cos(x[1]+B)/(l^2*m+I)))*1/C
end
dx[1] = x[2]
dx[2] = theta_dotdot
dx[3] = x[4]
dx[4] = x_dotdot
end
m = 2300 # [kg] mass of system
l = 10 # [m] height of structure
g = 9.81 # [m/s] gravity
kc = 795774.71/1000 # [N/rad] rotational stiffness of structure
kb = 236 # [N/m] stiffness of the rubber buffer
acceleration = 0.25 # [m/s] acceleration of the system
time_period = 2.5 # [s] time period which the above acceleration exists, at this time point imapct will occur
C1 = 0.04*(m+kb) # damping
C2 = 0.04*(m+kc) # damping
I = m*(l/2)^2 # mass moment of inertia assuming a point mass
distance_travelled = 0.5*acceleration*time_period^2 # total distance travelled during acceleration
print(distance_travelled)
x0 = [0.0 , 0.0, 0.0, 0] # [rad,rad/s,m,m/s] initial conditions
tspan = (0.0,5) # time for total simulation
p = (m,l,g,I,kc,kb,C1,C2,distance_travelled,acceleration,time_period)
prob = ODEProblem(statespace!, x0, tspan, p)
sol = solve(prob,Rodas5(),saveat = 0.001,reltol=1e-8,abstol=1e-8)
x = sol[3,:]
delta = x.-distance_travelled
plotly()
plt_1 = plot()
plot!(plt_1, sol.t,sol[1,:],xlabel = "Time (sec)",ylabel = "Angular displacement")
plot!(plt_1, [time_period, time_period], [maximum(sol[1,:]),minimum(sol[1,:])],label = "Impact")
plot!(plt_1, [decoupling[1],decoupling[1]], [maximum(sol[1,:]),minimum(sol[1,:])],label = "Decoupling")
display(plot(plt_1))
savefig(plt_1,"Angular displacement.html")
plt_2 = plot()
plot!(plt_2, sol.t,sol[2,:],xlabel = "Time (sec)",ylabel = "Angular velocity")
plot!(plt_2, [time_period, time_period], [maximum(sol[2,:]),minimum(sol[2,:])],label = "Impact")
plot!(plt_2, [decoupling[1],decoupling[1]], [maximum(sol[2,:]),minimum(sol[2,:])],label = "Decoupling")
display(plot(plt_2))
savefig(plt_2,"Angular velocity")
plt_3 = plot()
plot!(plt_3, sol.t,sol[3,:],xlabel = "Time (sec)",ylabel = "Linear displacement")
plot!(plt_3, [time_period, time_period], [maximum(sol[3,:]),minimum(sol[3,:])],label = "Impact")
plot!(plt_3, [decoupling[1],decoupling[1]], [maximum(sol[3,:]),minimum(sol[3,:])],label = "Decoupling")
plot!(plt_3, [0,5], [distance_travelled,distance_travelled],label = "distance travelled to impact")
display(plot(plt_3))
savefig(plt_3,"Linear displacement")
plt_4 = plot()
plot!(plt_4, sol.t,sol[4,:],xlabel = "Time (sec)",ylabel = "Linear velocity")
plot!(plt_4, [time_period, time_period], [maximum(sol[4,:]),minimum(sol[4,:])],label = "Impact")
plot!(plt_4, [decoupling[1],decoupling[1]], [maximum(sol[4,:]),minimum(sol[4,:])],label = "Decoupling")
display(plot(plt_4))
savefig(plt_4,"Linear velocity")
plt_5 = plot()
plot!(plt_5, sol.t,delta,xlabel = "Time (sec)",ylabel = "Linear displacement")
plot!(plt_5, [time_period, time_period], [maximum(sol[3,:]),minimum(sol[3,:])],label = "Impact")
plot!(plt_5, [decoupling[1],decoupling[1]], [maximum(sol[3,:]),minimum(sol[3,:])],label = "Decoupling")
display(plot(plt_5))
savefig(plt_5,"1")