агм дз
unknown
plain_text
2 months ago
9.9 kB
3
Indexable
Never
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Используемые библиотеки \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
import numpy as np
from scipy.optimize import root
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Используемые функции \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
# Функции атмосферы
def geopotencial_v():
return (r * h) / (r + h) # Геопотенциальная высота
def yskorenie_s_p():
return g0 * (r / (r + h)) ** 2 # Ускорение свободного падения
def parametry_t():
H = geopotencial_v()
# Параметры, для определения температуры и давления
if (H >= -2001) and (H < 0): # 3-4 стр. ГОСТа
Hx = -2000
batm = -0.0065
Tx = 301.15
px = 127774
elif (H >= 0) and (H <= 11000): # 5-15 стр. ГОСТа
Hx = 0
batm = -0.0065
Tx = 288.15
px = 101325
elif (H > 11000) and (H < 20000): # 15-24 стр. ГОСТа
Hx = 11000
Tx = 216.65
batm = 0
px = 22632
elif (H >= 20000) and (H < 32000): # 25-36 стр. ГОСТа
Hx = 20000
batm = 0.001
Tx = 216.65
px = 5474.87
elif (H >= 32000) and (H < 47000): # 37-43 стр. ГОСТа
Hx = 32000
batm = 0.0028
Tx = 228.65
px = 868.014
elif (H >= 47000) and (H < 51000): # 43-46 стр. ГОСТа
Hx = 47000
Tx = 270.65
batm = 0
px = 110.906
elif (H >= 51000) and (H < 71000): # 45-52 стр. ГОСТа
Hx = 51000
batm = -0.0028
Tx = 270.650
px = 66.9384
elif (H >= 71000) and (H < 85000): # 51-53 и 159 стр. ГОСТа
Hx = 71000
batm = -0.002
Tx = 214.65
px = 3.95639
return dict(Hx = Hx, batm = batm, Tx = Tx, px = px)
def temperatura():
p = parametry_t()
H = geopotencial_v()
T = p["Tx"] + p["batm"] * (H - p["Hx"]) # Температура
return T
def davlenie():
p = parametry_t()
H = geopotencial_v()
T = temperatura()
if p["batm"] != 0: # Давление
p = 10 ** (np.log10(p["px"]) - (g0 / (p["batm"] * R)) * np.log10(T / p["Tx"]))
else:
p = 10 ** (np.log10(p["px"]) - ((0.434294 * g0 * (H - p["Hx"])) / (T * R)))
return p
def plotnost():
p = davlenie()
T = temperatura()
ro = p / (R * T) # Плотность
return ro
def skorost_z():
T = temperatura()
a = 20.046796 * np.sqrt(T) # Скорость звука
return a
def vivod():
H = geopotencial_v()
g = yskorenie_s_p()
T = temperatura()
p = davlenie()
ro = plotnost()
a = skorost_z()
print("Геометрическая высота: " + "{:.0f}".format(h) + ' м')
print("Геопотенциальная высота: " + "{:.0f}".format(H) + ' м')
print("Ускорение свободного падения: "+str(round(g,3))+' м/с^2')
print("Температура: "+str(round(T/100,3))+' * 10^2 K')
print("Давление: "+str(round(p/10000,3))+' * 10^4'+' Па')
print('Плотность: '+str(round(ro*10,3))+' * 10^-1'+' кг/м^3')
print("Скорость звука: "+str(round(a/100,3))+' * 10^2 м/с')
print(' ')
# Функции для расчета:
# плоского СУ
def findShockWaveAngle(betta: float) -> float: # Функция поиска угла скачка уплотнения
def ugol_SY(tettac, bet): # Трансцендентное уравнение
return np.tan(tettac) / np.tan(tettac - bet) - ((k + 1) * M ** 2 * (np.sin(tettac)) ** 2) / (
2 + (k - 1) * M ** 2 * (np.sin(tettac)) ** 2)
# Генератор листа численных решений функции, с шагом 0.1. Решение находится 50 раз. Округление до 4 знаков после запятой
res = [round(root(ugol_SY, x0=i * 0.1, args=(betta)).x[0], 4) for i in range(50)]
# Удалены повторяющиеся, отсортировано по возрастанию
listRoots = sorted(list(set(res)))
#print(res) # Раскомментировать, чтобы посмотреть вывод всех корней
# Поиск наименьшего положительного корня
closestRoot = None
for listRoot in listRoots:
if listRoot > 0 and (closestRoot is None or listRoot < closestRoot):
closestRoot = listRoot
shockWaweAngle = closestRoot
return shockWaweAngle # Угол СУ
# Касательной скорости
def findVr(θ_c):
return M * skorost_z() * np.cos(θ_c)
# Нормальной скорости
def findVt(θ_c):
return - findVr(θ_c) * np.tan(θ_c) * (
(2 + (k - 1) * M ** 2 * np.sin(θ_c) ** 2) / ((k + 1) * M ** 2 * np.sin(θ_c) ** 2))
# Метода последовательных приближений
def V_r_next(V_r_previous, V_θ_previous, shag):
return V_r_previous + V_θ_previous * shag
def V_tetta_next(θ_c, V_r_previous, V_θ_previous, shag):
T_0 = T * (1 + ((k - 1) / 2) * M ** 2)
a_kr = (k * R * T_0) - ((k - 1) / 2) * (V_r_previous ** 2 - V_θ_previous ** 2)
dV0_dtetta = (-V_θ_previous / np.tan(θ_c) + V_r_previous * ((V_θ_previous ** 2) / a_kr - 2)) / (
1 - (V_θ_previous ** 2) / a_kr)
return V_θ_previous + dV0_dtetta * shag
def calc_block(θ_c, V_r, V_θ, shag):
V_r = V_r_next(V_r, V_θ, shag)
V_θ = V_tetta_next(θ_c, V_r, V_θ, shag)
θ_c += shag
return V_r, V_θ, θ_c
def ConeParametrs():
coef = 0.85
shag = -0.00001
θ_c_plosk = findShockWaveAngle(β_k)
θ_c = coef * θ_c_plosk
while round(np.degrees(θ_c), 2) != np.degrees(β_k):
θ_c = coef * findShockWaveAngle(β_k)
V_r = findVr(θ_c)
V_θ = findVt(θ_c)
#print(θ_c, ' ',V_r, ' ', V_θ)
while V_θ < 0:
V_r = V_r_next(V_r, V_θ, shag)
V_θ = V_tetta_next(θ_c, V_r, V_θ, shag)
θ_c += shag
#print('coef', coef, '|| β_n = ', round(np.degrees(θ_c), 2), '| V_r = ', V_r, '| V_θ = ', V_θ)
if round(θ_c, 2) > β_k:
coef -= 0.0001
elif round(θ_c, 2) < β_k:
coef += 0.0001
θ_k = coef * θ_c_plosk
β_n = θ_c
V_k = V_r
return θ_k, β_n, V_k # Угол θ корнуса, Расчетный угол конуса, # Нормальная скорость
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Исходные данные \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
N = 11
dm = 1.4
P = 0
if (N > 0) and (N < 6):
P += 1
elif (N >= 6) and (N < 11):
P += 2
elif (N >= 11) and (N < 16):
P += 3
elif (N >= 16) and (N < 20):
P += 4
M0 = 2
H0 = 3.5
βk0 = 7.5
x_cm = 0.65
β_k = np.radians(βk0 + 2.5 * (N - 2 * P)) # Угол конуса
lk = dm / (2 * np.tan(np.radians(β_k))) # Длина конуса
Xcm = x_cm * lk # Координата центра масс
l_s_chertoy_k = (2 * lk) / 3 # Длина затупленного конуса
# Данные для атмосферы
h = (H0 + 0.2 * N) * 1000 # Геометрическая высота
M = M0 + (P - 1) # Мах
r = 6356767
g0 = 9.80665
R = 287.05287
k = 1.4
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Определение параметров стандартной атмосферы \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
H = geopotencial_v()
g = yskorenie_s_p()
T = temperatura()
p = davlenie()
ro = plotnost()
a = skorost_z()
V = a * M # Скорость потока
vivod()
print("Мах: " + str(round(M ,3)) )
print("Скорость V: " + str(round(V,3)) + ' м/с')
print("Угол полураствора конуса: " + str(round(β_k, 3)) + ' рад; ' + str(round(np.degrees(β_k), 2)) + ' град')
print("Угол θ плоск: " + str(round(findShockWaveAngle(β_k), 3)) + ' рад; ' + str(round(np.degrees(findShockWaveAngle(β_k)), 2)) + ' град')
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Метод последовательных приближений \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
list_parametrov = ConeParametrs()
θ_c = list_parametrov[0] # Угол θ корнуса
β_n = list_parametrov[1] # Расчетный угол конуса
V_k = list_parametrov[2] # Нормальная скорость
print("Угол θ корнуса: " + str(round(θ_c, 3)) + ' рад; ' + str(round(np.degrees(θ_c), 3)) + ' град')
print("Расчетный угол конуса: " + str(round(β_n, 3)) + ' рад; ' + str(round(np.degrees(β_n), 3)) + ' град')
print("Нормальная скорость V_r: " + str(round(V_k,3)) + ' м/с')
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Расчёт p̅к, Mк, Cx \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
T_0 = T * (1 + ((k - 1) / 2) * M ** 2)
V_max = np.sqrt(2 * k * R * T_0 / (k - 1))
T_k = T_0 * (1 - (V_k / V_max) ** 2)
M_k = V_k / np.sqrt(k * R * T_k)
P_k = (1 + (k - 1) * M_k ** 2 / 2) ** ((- k)/(k - 1))
p̅k = (P_k - p) * 2 / (ro * V ** 2)
print("Коэффициент давления p̅k: " + str(round(p̅k ,3)) )
#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Расчёт по приближенным зависимостям \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
θ_c_emp = np.sqrt(np.arcsin((1 / M ** 2) ))