агм дз
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) ))