Untitled
unknown
plain_text
a month ago
7.6 kB
2
Indexable
Never
# Библиотеки
import numpy as np
# Константы
R_luna = 1738000 # [м] Луна - сфера радиуса R_m
mu_luna = 4.903 * (10 ** 12) # м^3/с^2 - гравитационный параметр Луны
m_0 = 3015 # [кг] начальная масса СВ
delta_m_t_max = 1550 # [кг] предельная масса рабочего запаса топлива ДУ СВ
m_const = 615 # [кг] масса конструкции СВ
t_v = 14 # [с] продолжительность вертикального участка выведения
delta_t = 0.1 # [с] шаг интегрирования
# Начальные данные
W_ist = 3420 # [м/с] эффективная скорость истечения топлива
P_1 = 9100 # [Н] тяга СВ 1
P_2 = 10600 # [Н] тяга СВ 2
h_isl1_1 = 125000 # [м]
h_isl1_2 = 200000 # [м]
h_isl2_1 = 200000 # [м]
h_isl2_2 = 355000 # [м]
# Задаем текущую высоту орбиты и тягу
h = h_isl1_1
P = P_1
# н.у.
Vx = 0
Vy = 0
x = 0
y = 0
t = 0
# углы для управления вектором тяги ДУ - радианы
theta_1_pr = -0.005 # radians
theta_2 = -0.05 # radians
# значения конечного времени для 1 и 2 участка
t_1 = 360.738
t_2 = 700.245
# проекции ускорения от грав. поля Луны
def g_x(x, y):
return - (mu_luna * x) / ((x ** 2 + (R_luna + y) ** 2) ** (3 / 2))
def g_y(x, y):
return - (mu_luna * (R_luna + y)) / ((x ** 2 + (R_luna + y) ** 2) ** (3 / 2))
# Закон изменения угла тангажа на активных участках
def theta(t, thata_1_pr, theta_2):
if t <= t_v:
return np.pi / 2
elif (t_v < t <= t_1):
return np.pi / 2 + thata_1_pr * (t - t_v)
elif (t_2 <= t):
return theta_2
else:
return 0
# мат.модель
def dVx_dt(t, P, m, x, y, theta_1_dot, theta_2):
return P * np.cos(theta(t, theta_1_dot, theta_2)) / m + g_x(x, y)
def dVy_dt(t, P, m, x, y, theta_1_dot, theta_2):
return P * np.sin(theta(t, theta_1_dot, theta_2)) / m + g_y(x, y)
def dm_dt(P):
return - (P / W_ist)
#РунгеКутта4
def RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt):
# расчет коэффициентов: k1, k2, k3, k4 через \n\n соответственно
Vx_k1 = dt * dVx_dt(t, P, m, x, y, theta_1_dot, theta_2)
Vy_k1 = dt * dVy_dt(t, P, m, x, y, theta_1_dot, theta_2)
x_k1 = dt * Vx
y_k1 = dt * Vy
m_k1 = dt * dm_dt(P)
Vx_k2 = dt * dVx_dt(t + dt / 2, P, m + m_k1 / 2, x + x_k1 / 2, y + y_k1 / 2, theta_1_dot, theta_2)
Vy_k2 = dt * dVy_dt(t + dt / 2, P, m + m_k1 / 2, x + x_k1 / 2, y + y_k1 / 2, theta_1_dot, theta_2)
x_k2 = dt * (Vx + Vx_k1 / 2)
y_k2 = dt * (Vy + Vy_k1 / 2)
m_k2 = dt * dm_dt(P)
Vx_k3 = dt * dVx_dt(t + dt / 2, P, m + m_k2 / 2, x + x_k2 / 2, y + y_k2 / 2, theta_1_dot, theta_2)
Vy_k3 = dt * dVy_dt(t + dt / 2, P, m + m_k2 / 2, x + x_k2 / 2, y + y_k2 / 2, theta_1_dot, theta_2)
x_k3 = dt * (Vx + Vx_k2 / 2)
y_k3 = dt * (Vy + Vy_k2 / 2)
m_k3 = dt * dm_dt(P)
Vx_k4 = dt * dVx_dt(t + dt, P, m + m_k3, x + x_k3, y + y_k3, theta_1_dot, theta_2)
Vy_k4 = dt * dVy_dt(t + dt, P, m + m_k3, x + x_k3, y + y_k3, theta_1_dot, theta_2)
x_k4 = dt * (Vx + Vx_k3)
y_k4 = dt * (Vy + Vy_k3)
m_k4 = dt * dm_dt(P)
# расчет новых величин в конце участка, соответствующего времени t + dt
Vx += (Vx_k1 + Vx_k2 * 2 + Vx_k3 * 2 + Vx_k4) / 6
Vy += (Vy_k1 + Vy_k2 * 2 + Vy_k3 * 2 + Vy_k4) / 6
x += (x_k1 + x_k2 * 2 + x_k3 * 2 + x_k4) / 6
y += (y_k1 + y_k2 * 2 + y_k3 * 2 + y_k4) / 6
m += (m_k1 + m_k2 * 2 + m_k3 * 2 + m_k4) / 6
return Vx, Vy, x, y, m
# Функция основного расчета параметров траектории
def integration(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt, t_1, t_2):
print('\nДанные выводятся в формате: Vx [м/с], Vy [м/с], x [м], y [м], m [кг], t [с]')
print('\nРАСЧЕТ ТОЧКИ t_v\n')
while t < t_v:
Vx, Vy, x, y, m = RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt)
t += dt
if round(t, 1) == t_v:
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
print('\nРАСЧЕТ ТОЧКИ t_1\n')
print(t_1 - (t_1 % 0.1))
while t < t_1:
if round(t, 1) == round(t_1 - (t_1 % 0.1),3):
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
# первый шаг с тягой до критической точки t_1
Vx, Vy, x, y, m = RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, t_1 % 0.1)
t += t_1 % 0.1
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
# второй шаг без тяги от критической точки t_1
Vx, Vy, x, y, m = RK4(t, 0, m, x, y, Vx, Vy, theta_1_dot, theta_2, delta_t - t_1 % 0.1)
t += delta_t - t_1 % 0.1
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
else:
Vx, Vy, x, y, m = RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt)
t += dt
#print(t,round(t, 1))
print('\nРАСЧЕТ ТОЧКИ t_2\n')
while t < t_2:
if round(t, 1) == t_2 - (t_2 % 0.1):
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
# первый шаг без тяги до критической точки t_2
Vx, Vy, x, y, m = RK4(t, 0, m, x, y, Vx, Vy, theta_1_dot, theta_2, t_2 % 0.1)
t += t_2 % 0.1
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
# второй шаг с тягой от критической точки t_2
Vx, Vy, x, y, m = RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, delta_t - t_2 % 0.1)
t += delta_t - t_2 % 0.1
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', round(t, 3))
else:
Vx, Vy, x, y, m = RK4(t, 0, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt)
t += dt
print('\nРАСЧЕТ 2АУТ И КОНЕЧНОЙ ТОЧКИ\n')
while abs(np.sqrt(Vx ** 2 + Vy ** 2) - (np.sqrt(mu_luna / (R_luna + h)))) > (10 ** (-7)):
Vx_do = Vx
Vy_do = Vy
x_do = x
y_do = y
m_do = m
Vx, Vy, x, y, m = RK4(t, P, m, x, y, Vx, Vy, theta_1_dot, theta_2, dt)
if np.sqrt(Vx ** 2 + Vy ** 2) > (np.sqrt(mu_luna / (R_luna + h))):
Vx = Vx_do
Vy = Vy_do
x = x_do
y = y_do
m = m_do
dt /= 10
else:
t += dt
if round(m, ) < m_0 - delta_m_t_max:
Vx = Vx_do
Vy = Vy_do
x = x_do
y = y_do
m = m_do
t -= dt
break
print(Vx, ',', Vy, ',', x, ',', y, ',', m, ',', t)
print('Полученная скорость:', np.sqrt(Vx ** 2 + Vy ** 2))
print('Действительная скорость:', np.sqrt(mu_luna / (R_luna + h)))
return Vx, Vy, x, y, m, t
Vx, Vy, x, y, m, t = integration(t, P, m_0, x, y, Vx, Vy, theta_1_pr, theta_2, delta_t, t_1, t_2)