Skip to content

Project: Plotting rocket motion for 3c #36

Open
dhk15102 opened this issue Apr 3, 2020 · 1 comment
Open

Project: Plotting rocket motion for 3c #36

dhk15102 opened this issue Apr 3, 2020 · 1 comment

Comments

@dhk15102
Copy link

dhk15102 commented Apr 3, 2020

I got my f_m function to work and tested it with my rocket function from earlier, but my plots with the new dmdt that the mod_secant function gives me are weird.

I copied my rocket function from above and just changed the name, but when I replace the mass flow rate in the function with the one I get for 3c the plot caps out around 350m instead of 300m. When I use the original mass flow rate of 0.05 kg/s in the function but use the calculated dmdt for the period and time step, the plot looks right. Any idea why this is happening?

def rocket2(state,dmdt=0.05, u=250,c=0.18e-3):
    g = 9.81
    dstate = np.zeros(np.shape(state))
    dstate[0] = state[1]
    dstate[1] = (u/state[2])*dmdt - g - (c/state[2])*state[1]**2
    dstate[2] = -dmdt
    return dstate
t = np.linspace(0,(m0-mf)/dmdt1,N)
dt = (m0 - mf)/(dmdt1*N)
num_sol_3c = np.zeros([N,3])
num_sol_3c[0,0] = y0
num_sol_3c[0,1] = v0
num_sol_3c[0,2] = m0
for i in range(N-1):
    num_sol_3c[i+1] = heun_step(num_sol_3c[i], rocket2, dt)


fig = plt.figure(figsize=(15,10))
plt.plot(t,num_sol_3c[:,0],linewidth=2, color='g', linestyle='-', label='Rocket')
plt.xlabel('time (s)')
plt.ylabel('y(t) (m)')
plt.legend()
plt.title('Altitude vs. time');

image

@rcc02007
Copy link
Owner

rcc02007 commented Apr 8, 2020

What is dmdt1 in the lines

dt = (m0 - mf)/(dmdt1*N)
and
dt = (m0 - mf)/(dmdt1*N)

are the values set to what you intended? If you want to change the dmdt in the ode you need

heun_step(num_sol_3c[i], lambda state: rocket2(state,dmdt=dmdt1), dt)

Sign in to join this conversation on GitHub.
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants