HW2 Question

[apolicow]

Cris,

I hope your work trip is going well. We had a question regarding HW 2 we'd like to ask. We used the "SpringMassRandom" code to model the homework and are getting a little confused with our results. As shown in the attached screenshot, the graph of force vs time does not exactly look random. Even though we are passing the random function "rand()" through this code, it still seems to have a steady state and transient response plot. Based on the video, we counted cycles over time and did the math to get wn = 8.9 rad/sec. Since the pole is moving at its natural frequency, we know the angular velocity of the wind (w) is very close to equaling wn. Therefore, we input w=8 rad/sec into the Julia code. What is the best way to model the forcing function to truly make it random? My code is included below. Thank you

#define Julia ODE package that we will use as solver
using OrdinaryDiffEq

#define constants
m=100 #kg
g=9.8 #m/sec^2
Tn=0.706 #seconds
fn=1/Tn
ωn=2pifn #radians/second

k=ωn^2*m #N/m

ζ=0.005 #viscous damping ratio
ccr=2mωn #critical viscous damping coefficient
c=ζ*ccr #viscous damping coefficient, units of kg/sec.

μ=0 #friction coefficient
N=mg #normal force
F=μN #friction force

#calculate frequency and period
ωn=sqrt(k/m) #radians per second
fn=ωn/(2*pi) #cycles per second
Tn=1/fn #period, seconds

#forcing function
a=1278.9 #N
ω=8 #radians/sec.

#write equation of motion
function SpringMassDamperEOM(ddu,du,u,p,t)
m, k, c, a, ω, F = p

#make the forcing function random
r=rand()   #generates a random number from a uniform distribution [0,1]

#harmonic forcing function
pt=a*r*sin(ω*t[1])

#friction damping
if du[1]>0
    SignF=1
else
    SignF=-1
end

ddu[1] = -c/m*du[1] - k/m*u[1] + pt/m -F*SignF/m

end

#solver controls
dt=0.001 #seconds
totalt=36 #seconds

#define initial conditions
uo=0.1 #meters
duo=0.0 #m/sec.

#solve the ODE
prob = SecondOrderODEProblem(SpringMassDamperEOM, [duo], [uo], (0.,totalt), (m, k, c, a, ω, F))
sol = solve(prob, DPRKN6(),tstops=0:dt:totalt)

#separate out u(t),du(t), and t from solution
du=first.(sol.u)
u=last.(sol.u) #angular velocity
t=sol.t

#plot results
using Plots #start Julia plot package
plot(t,u,linewidth=1,
xaxis="t [sec.]",yaxis="u [meters]", legend=false)

[cristophermoen]

Hi Aaron,

Everything is looking reasonable to me. You are seeing vibration primarily at the natural frequency. Slight changes in the damping ratio, or in omega (the forcing function) change the response. If I set w=wn for example, you see resonance and deformation growth to a constant, even with the random forcing function. If I make w=1 rad/sec, you can see the transient and steady state response. Also, if I change the damping ratio, the response can vary significantly.

I think your specific plot is showing the influence of the random load, not a transient response riding on a steady state response. You can become more comfortable with your results by studying some benchmarking problems like those in Chopra Fig. 3.2.5.

Cris