using PyPlot rc("font",family ="Osaka",size=16) using ODE function gravity(altitude) # altitude: meter g0 = 9.80665 radius_earth = 6378137 g = g0 * (radius_earth / (radius_earth + altitude)) ^2 end # gravity関数のテスト test_alt = 1:1000:10000000 test_g = zeros(length(test_alt)) index = 1 for i = test_alt test_g[index] = gravity(i) index += 1 end plot(test_alt/1000, test_g) ylim(0,10) xlabel("高度 km");ylabel("重力加速度 kg/s2");grid() function air(altitude) # 標準大気モデルを用いた、高度による温度、音速、大気圧、空気密度の関数 # 高度は基準ジオポテンシャル高度を元にしている。 # 標準大気の各層ごとの気温減率から定義式を用いて計算している。 # Standard Atmosphere 1976 ISO 2533:1975 # 中間圏高度86kmまでの気温に対応している。それ以上は国際標準大気に当てはまらないので注意。 # cf. http://www.pdas.com/hydro.pdf # @param altitude 高度[m] # @return T 温度[K] # @return a 音速[m/s] # @return P 気圧[Pa] # @return rho 空気密度[kg/m3] h = altitude g = 9.80655 gamma = 1.4 R = 287.0531 # height of atomospheric layer HAL = [0, 11000, 20000, 32000, 47000, 51000, 71000, 84852] # Lapse Rate Kelvin per meter LR = [-0.0065, 0.0, 0.001, 0.0028, 0, -0.0028, -0.002, 0.0] # Tempareture Kelvin T0 = [288.15, 216.65, 216.65, 228.65, 270.65, 270.65, 214.65, 186.95] # Pressure Pa P0 = [101325, 22632, 5474.9, 868.02, 110.91, 66.939, 3.9564, 0.3734] if ( h >= HAL[1] && h < HAL[2]) k = 1 elseif ( h >= HAL[2] && h < HAL[3]) k = 2 elseif ( h >= HAL[3] && h < HAL[4]) k = 3 elseif ( h >= HAL[4] && h < HAL[5]) k = 4 elseif ( h >= HAL[5] && h < HAL[6]) k = 5 elseif ( h >= HAL[6] && h < HAL[7]) k = 6 elseif ( h >= HAL[7] && h < HAL[8]) k = 7 elseif ( h >= HAL[8]) k = 8 else k = 1 end T = T0[k] + LR[k] * (h - HAL[k]) a = sqrt(T * gamma * R) if LR[k] != 0 P = P0[k] * ((T0[k] + LR[k] * (h - HAL[k])) / T0[k]) ^ (g / -LR[k] / R) else P = P0[k] * exp(g / R * (HAL[k] - h) / T0[k]) end rho = P / R / T [T, a, P, rho] end # air関数のテスト test_alt = 1:100:100000 test_air = zeros(length(test_alt)) index = 1 for i = test_alt T, a, P, test_air[index] = air(i) index += 1 end plot(test_alt/1000, test_air) xlabel("高度 km");ylabel("空気密度 kg/m3");grid() # using PyCall # @pyimport scipy.interpolate as interp # 補間が使えなかった。。。 function Cd(mach) # M-Vロケットの抗力係数を参考にした # @param Mach マッハ数[-] # @return self.Cd 抗力係数[-] M = [0, 0.2, 0.4, 0.6, 0.8, 1.0, 1.1, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0]; cd_row = [0.28, 0.28, 0.28, 0.29, 0.35, 0.64, 0.67, 0.69, 0.66, 0.62, 0.58, 0.55, 0.48, 0.42, 0.38, 0.355, 0.33] # 線形補間のアルゴリズム、外挿はせずに for i = 1:length(M)-1 if mach >= M[i] && mach < M[i+1] alpha = (mach - M[i]) / (M[i+1] - M[i]) Cd = cd_row[i] + alpha * (cd_row[i+1] - cd_row[i]) break end end if mach < M[1] Cd = cd_row[1] elseif mach >= M[end] Cd = cd_row[end] end Cd end test_mach = 0:0.1:7 test_Cd = zeros(length(test_mach)) index = 1 for i = test_mach test_Cd[index] = Cd(i) index += 1 end plot(test_mach, test_Cd) ylim(0,1) xlabel("マッハ数 -");ylabel("Cd -");grid() # ロケットTypeを作る。ダイナミクスなどの関数にまとめてパラメータを入れるため type Rocket vy::Float64 y::Float64 Isp::Float64 thrust::Float64 burn_time::Float64 mass_init::Float64 area::Float64 end rocket = Rocket(0.0,0.0,200.0,500.0,3.0,10.0,0.5) function drag(rocket::Rocket) T, a, P, rho = air(rocket.y) mach = rocket.vy / a drag = 0.5 * Cd(mach) * rho * rocket.area * rocket.vy ^2 end test_vel = 0:10:700 test_drag = zeros(length(test_vel)) index = 1 for i = test_vel rocket = Rocket(i, 0, 200, 500, 3, 10, 0.12) test_drag[index] = drag(rocket) index += 1 end plot(test_vel, test_drag, label="面積0.12m2のとき") # ylim(0,1) xlabel("速度 m/s");ylabel("抗力 N");grid() legend(loc="best") function flight_dynamics(t, x, rocket::Rocket) # x = [mass, vy, y] # type Rocket : vy, y, Isp, thrust, burn_time, mass_init, area dx = similar(x) g0 = 9.80655 if t < rocket.burn_time thrust = rocket.thrust else thrust = 0 end m_dot = thrust / rocket.Isp / g0 rocket.vy = x[2]; rocket.y = x[3] if x[2] > 0 D = drag(rocket) else D = - drag(rocket) end dx[1] = -m_dot dx[2] = 1 / x[1] * (thrust - D) - gravity(x[3]) dx[3] = x[2] if x[3] < 0.0 dx[1] = 0; dx[2] = 0; dx[3] = 0 end dx end t = [0., 400] Isp = 200; thrust = 31800; burn_time = 30; mass_init = 722; area = 0.07 # S-310 rocket rocket = Rocket(0.0, 0.0, Isp, thrust, burn_time, mass_init, area) @time t,y = ode23s((t, x) -> flight_dynamics(t,x,rocket), [rocket.mass_init, 0.0, 0.0], t; abstol=1e-5, reltol=1e-5, maxstep=1e11/10, minstep=1e11/1e18) Void # ==== PLOT ==== mass = zeros(length(t)) velocity = zeros(length(t)) altitude = zeros(length(t)) for i = 1:length(t) mass[i] = y[i][1] velocity[i] = y[i][2] altitude[i] = y[i][3] end time_S310 = t mass_S310 = mass vel_S310 = velocity alt_S310 = altitude subplot(311) title("JAXA S-310") plot(t,mass); ylim(0, rocket.mass_init*1.1) ylabel("質量 kg");grid() subplot(312) plot(t, velocity) ylabel("速度 m/s") axhline(330, linewidth=1, linestyle = "--") subplot(313) plot(t, altitude/1000) ylabel("高度 km") xlabel("時間 s");grid() t = [0., 400] Isp = 360; thrust = 420000; burn_time = 120; mass_init = 29000; area = 4.9 # Blue Origin New Shepard rocket = Rocket(0.0, 0.0, Isp, thrust, burn_time, mass_init, area) @time t,y = ode23s((t, x) -> flight_dynamics(t,x,rocket), [rocket.mass_init, 0.0, 0.0], t; abstol=1e-4, reltol=1e-4, maxstep=1e11/10, minstep=1e11/1e18) Void # ==== PLOT ==== mass = zeros(length(t)) velocity = zeros(length(t)) altitude = zeros(length(t)) for i = 1:length(t) mass[i] = y[i][1] velocity[i] = y[i][2] altitude[i] = y[i][3] end time_NS = t mass_NS = mass vel_NS = velocity alt_NS = altitude subplot(311) title("Blue Origin New Shepard") plot(t,mass); ylim(0, rocket.mass_init*1.1) ylabel("質量 kg");grid() subplot(312) plot(t, velocity) ylabel("速度 m/s");grid() axhline(330, linewidth=1, linestyle = "--") subplot(313) plot(t, altitude/1000) ylabel("高度 km") xlabel("時間 s");grid() subplot(311) title("S-310 VS New Shepard") plot(time_S310,mass_S310, label="S-310") plot(time_NS,mass_NS, label="New Shepard") ylim(0,30000);xlim(0,300) legend(loc = "best", fontsize = 10) # plot(time_NS,mass_NS) ylabel("質量 kg") subplot(312) plot(time_S310, vel_S310, label="S-310") plot(time_NS, vel_NS, label="New Shepard") ylabel("速度 m/s") axhline(330, linewidth=1, linestyle = "--") ylim(-500, 2000);grid();xlim(0,300) subplot(313) plot(time_S310, alt_S310/1000, label="S-310") plot(time_NS, alt_NS/1000, label="New Shepard") ylabel("高度 km") xlabel("時間 s") grid();ylim(0,200);xlim(0,300)