main.py
Основной код программы
unknown
plain_text
a year ago
18 kB
4
Indexable
import numpy as np
from scipy.optimize import fsolve
from scipy.optimize import root
from functions import atmosphere as atm
import data as d
# Функции для решения
def findShockWaveAngle(betta: float) -> float: # Функция поиска угла скачка уплотнения
def ugol_SY(tettac, betta): # Трансцендентное уравнение
return np.tan(tettac) / np.tan(tettac - betta) - ((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 findFictitiousAngle (M: float) -> float: # Функция поиска фиктивного угла
FictitiousAngle = np.sqrt((k + 1) / (k - 1)) * np.arctan(np.sqrt(((k - 1) * (M ** 2 - 1)) / (k + 1))) - np.arctan(np.sqrt(M ** 2 - 1))
return FictitiousAngle # Фиктивный угол
def findMachThroughFictitiousAngle(w: float) -> float: # Функция поиска Маха, через фиктивный угол
def phiktivnyi_u(M, w) -> float: # Трансцендентное уравнение
return np.sqrt((k + 1) / (k - 1)) * np.arctan(np.sqrt(((k - 1) * (M ** 2 - 1)) / (k + 1))) - np.arctan(np.sqrt(M ** 2 - 1)) - w
# Генератор листа численных решений функции, с шагом 0.1. Решение находится 50 раз. Округление до 4 знаков после запятой
res = [round(root(phiktivnyi_u, x0=i * 1, args=(w)).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 findDeflectionAngle(p): # Функция поиска углов СУ 5 и 6, и угла отклонения
tettac5,tettac6,delta = p
return (
np.tan(tettac5) / np.tan(tettac5 - betta - delta) - ((k + 1) * M3 ** 2 * (np.sin(tettac5)) ** 2) / (2 + (k - 1) * M3 ** 2 * (np.sin(tettac5)) ** 2),
np.tan(tettac6) / np.tan(tettac6 - betta + delta) - ((k + 1) * M4 ** 2 * (np.sin(tettac6)) ** 2) / (2 + (k - 1) * M4 ** 2 * (np.sin(tettac6)) ** 2),
p3 * ((2 * k * M3**2 * np.sin(tettac5)**2 - (k-1))/(k+1)) - p4 * ((2 * k * M4**2 * np.sin(tettac6)**2 - (k-1))/(k+1))
)
# Параметры потока в скачках уплотнения
def findFlowParameterPressure(p: float, M: float, tettac: float) -> float: # Функция поиска давления
FlowParameterPressure = p * ((2 * k * M**2 * np.sin(tettac)**2 - (k-1))/(k+1))
return FlowParameterPressure # Давление области
def findFlowParameterDensity(ro: float, M: float,tettac: float) -> float: # Функция поиска плотности
FlowParameterDensity = (ro * (k+1) * M**2 * np.sin(tettac)**2)/(2 + ((k-1) * M**2 * np.sin(tettac)**2))
return FlowParameterDensity # Плотность области
def findFlowParameterTemperature(p: float, ro: float, pi: float, roi: float) -> float: # Функция поиска температуры
FlowParameterTemperature = (T * pi * ro)/(p * roi)
return FlowParameterTemperature # Температура области
def findFlowParameterVelocity(V: float, betta: float, tettac: float) -> float: # Функция поиска скорости
FlowParameterVelocity = (V * np.cos(tettac))/(np.cos(tettac-betta))
return FlowParameterVelocity # Скорость области
def findFlowParameterАcceleration(T: float) -> float: # Функция поиска скорости
FlowParameterАcceleration = np.sqrt(k * R * T)
return FlowParameterАcceleration # Скорость области
def findFlowParameterMach(V: float, a: float) -> float: # Функция поиска маха
FlowParameterMach = V / a
return FlowParameterMach # Число маха области
# Параметры торможения
def findBrakingParameterPressure(p: float,M: float) -> float: # Функция поиска давления торможения
BrakingParameterPressure = p * (((k - 1) * M ** 2) / 2 + 1) ** (k / (k - 1))
return BrakingParameterPressure # Давление торможения
def findBrakingParameterDensity(ro: float, M: float) -> float: # Функция поиска плотности торможения
BrakingParameterDensity = ro * (((k - 1) * M ** 2) / 2 + 1) ** (1 / (k - 1))
return BrakingParameterDensity # Плотность торможения
def findBrakingParameterTemperature(T: float, M: float) -> float: # Функция поиска температуры торможения
BrakingParameterTemperature = T * (((k - 1) * M ** 2) / 2 + 1)
return BrakingParameterTemperature # Температура торможения
# Получение данных из data.py
M = d.M
alfa = (d.alfaDegrees * np.pi) / 180
betta = d.betta
chordLength = d.chordLength
h = d.h
TempCooledWall = d.TempCooledWall
criticalReynoldsNumber = d.criticalReynoldsNumber
R = d.R
# Получение данных из atmosphere.py
H = atm.geopotencial_v()
g = atm.yskorenie_s_p()
T = atm.temperatura()
p = atm.davlenie()
ro = atm.plotnost()
a = atm.skorost_z()
V = a * M # Скорость потока
# 1 и 2 области
# Углы отклонения потока
betta1 = betta - alfa
betta2 = betta + alfa
# Поиск углов СУ
k = 1.4
tettac1 = findShockWaveAngle(betta1)
tettac2 = findShockWaveAngle(betta2)
# Определяем параметры потока в скачках уплотнения
p1 = findFlowParameterPressure(p,M,tettac1)
p2 = findFlowParameterPressure(p,M,tettac2)
ro1 = findFlowParameterDensity(ro,M,tettac1)
ro2 = findFlowParameterDensity(ro,M,tettac2)
T1 = findFlowParameterTemperature(p,ro,p1,ro1)
T2 = findFlowParameterTemperature(p,ro,p2,ro2)
V1 = findFlowParameterVelocity(V,betta1,tettac1)
V2 = findFlowParameterVelocity(V,betta2,tettac2)
a1 = findFlowParameterАcceleration(T1)
a2 = findFlowParameterАcceleration(T2)
M1 = findFlowParameterMach(V1, a1)
M2 = findFlowParameterMach(V2, a2)
# Определяем параметры торможения для участка
p01 = findBrakingParameterPressure(p1,M1)
p02 = findBrakingParameterPressure(p2,M2)
ro01 = findBrakingParameterDensity(ro1,M1)
ro02 = findBrakingParameterDensity(ro2,M2)
T01 = findBrakingParameterTemperature(T1,M1)
T02 = findBrakingParameterTemperature(T2,M2)
# 3 и 4 области
# Прямая задача течения Прандтля – Майера
betta3 = 2 * betta
betta4 = 2 * betta
# Фиктивные углы и Мах
w1 = findFictitiousAngle(M1)
w3 = w1 + betta3
w2 = findFictitiousAngle(M2)
w4 = w2 + betta4
M3 = findMachThroughFictitiousAngle(w3)
M4 = findMachThroughFictitiousAngle(w4)
# Поскольку поток является невязким, давление, плотность и температура торможения в областях 1 и 3 (2 и 4) равны
p03 = p01
p04 = p02
ro03 = ro01
ro04 = ro02
T03 = T01
T04 = T02
# Найдем параметры потока
p3 = 1 / findBrakingParameterPressure(1/p03,M3) # 1/х используется, для поиска давления, но по формуле давления торможения
p4 = 1 / findBrakingParameterPressure(1/p04,M4)
ro3 = 1 / findBrakingParameterDensity(1/ro03,M3)
ro4 = 1 / findBrakingParameterDensity(1/ro04,M4)
T3 = 1 / findBrakingParameterTemperature(1/T01,M3)
T4 = 1 / findBrakingParameterTemperature(1/T02,M4)
a3 = findFlowParameterАcceleration(T3)
a4 = findFlowParameterАcceleration(T4)
V3 = findFlowParameterMach(M3, 1/a3) # используя формулу Маха, найдем скорость
V4 = findFlowParameterMach(M4, 1/a4)
# 5 и 6 области
betta5 = betta
betta6 = betta
t1, t2, u = 0.3, 0.3, 0.17 # Начальные приближенные значение
tettac5,tettac6,delta = fsolve(findDeflectionAngle, (t1, t2, u))
# Определяем параметры потока в скачках уплотнения
p5 = findFlowParameterPressure(p3,M3,tettac5)
p6 = findFlowParameterPressure(p4,M4,tettac6)
V5 = findFlowParameterVelocity(V3,betta5,tettac5)
V6 = findFlowParameterVelocity(V4,betta6,tettac6)
ro5 = findFlowParameterDensity(ro3,M3,tettac5)
ro6 = findFlowParameterDensity(ro3,M4,tettac6)
T5 = findFlowParameterTemperature(p3,ro3,p5,ro5)
T6 = findFlowParameterTemperature(p4,ro4,p6,ro6)
a5 = findFlowParameterАcceleration(T5)
a6 = findFlowParameterАcceleration(T6)
M5 = findFlowParameterMach(V5, a5)
M6 = findFlowParameterMach(V6, a6)
# Определяем параметры торможения для участка
p05 = findBrakingParameterPressure(p5,M5)
p06 = findBrakingParameterPressure(p6,M6)
ro05 = findBrakingParameterDensity(ro5,M5)
ro06 = findBrakingParameterDensity(ro6,M6)
T05 = findBrakingParameterTemperature(T5,M5)
T06 = findBrakingParameterTemperature(T6,M6)
#Проверочный вывод
atm.vivod()
print("Угол бетта: " + str(round(betta,3))+ ' рад; '+ str(round((betta*180)/np.pi,2))+ ' град')
print("Угол альфа: " + str(round(alfa,3))+ ' рад; '+ str(round((alfa*180)/np.pi,2))+ ' град')
print("Угол бетта1: " + str(round(betta1,3))+ ' рад; '+ str(round((betta1*180)/np.pi,2))+ ' град')
print("Угол бетта2: " + str(round(betta2,3))+ ' рад; '+ str(round((betta2*180)/np.pi,2))+ ' град')
print("Угол теттац1: " + str(round(tettac1,3))+ ' рад; '+ str(round((tettac1*180)/np.pi,2))+ ' град')
print("Угол теттац2: " + str(round(tettac2,3))+ ' рад; '+ str(round((tettac2*180)/np.pi,2))+ ' град')
print("Давление 1 области: " + str(round(p1/10000,3)) + ' * 10^4'+' Па')
print("Давление 2 области: " + str(round(p2/10000,3)) + ' * 10^4'+' Па')
print("Плотность 1 области: " + str(round(ro1*10,3))+' * 10^-1'+' кг/м^3')
print("Плотность 2 области: " + str(round(ro2*10,3))+' * 10^-1'+' кг/м^3')
print("Температура 1 области: " + str(round(T1/100,3)) + ' * 10^2 K')
print("Температура 2 области: " + str(round(T2/100,3)) + ' * 10^2 K')
print("Скорость 1 области: " + str(round(V1/1000,3)) + ' * 10^3 м/с')
print("Скорость 2 области: " + str(round(V2/1000,3)) + ' * 10^3 м/с')
print("Скорость звука 1 области: " + str(round(a1/100,3))+' * 10^2 м/с')
print("Скорость звука 2 области: " + str(round(a2/100,3))+' * 10^2 м/с')
print("Число Маха 1 области: " + str(round(M1,3)))
print("Число Маха 2 области: " + str(round(M2,3)))
print("Давление торможения 1 области: " + str(round(p01/10000000,3)) + ' * 10^7'+' Па')
print("Давление торможения 2 области: " + str(round(p02/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 1 области: " + str(round(ro01/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 2 области: " + str(round(ro02/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 1 области: " + str(round(T01/1000,3)) + ' * 10^3 K')
print("Температура торможения 2 области: " + str(round(T02/1000,3)) + ' * 10^3 K')
print("Угол бетта3: " + str(round(betta3,3))+ ' рад; '+ str(round((betta3*180)/np.pi,2))+ ' град')
print("Угол бетта4: " + str(round(betta4,3))+ ' рад; '+ str(round((betta4*180)/np.pi,2))+ ' град')
print("Фиктивный угол 1: " + str(round(w1,3))+ ' рад; '+ str(round((w1*180)/np.pi,2))+ ' град')
print("Фиктивный угол 2: " + str(round(w2,3))+ ' рад; '+ str(round((w2*180)/np.pi,2))+ ' град')
print("Фиктивный угол 3: " + str(round(w3,3))+ ' рад; '+ str(round((w3*180)/np.pi,2))+ ' град')
print("Фиктивный угол 4: " + str(round(w4,3))+ ' рад; '+ str(round((w4*180)/np.pi,2))+ ' град')
print("Число Маха 3 области: " + str(round(M3,3)))
print("Число Маха 4 области: " + str(round(M4,3)))
print("Давление торможения 3 области: " + str(round(p03/10000000,3)) + ' * 10^7'+' Па')
print("Давление торможения 4 области: " + str(round(p04/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 3 области: " + str(round(ro03/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 4 области: " + str(round(ro04/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 3 области: " + str(round(T03/1000,3)) + ' * 10^3 K')
print("Температура торможения 4 области: " + str(round(T04/1000,3)) + ' * 10^3 K')
print("Давление 3 области: " + str(round(p3/1000,3)) + ' * 10^3'+' Па')
print("Давление 4 области: " + str(round(p4/1000,3)) + ' * 10^3'+' Па')
print("Плотность 3 области: " + str(round(ro3*100,3))+' * 10^-2'+' кг/м^3')
print("Плотность 4 области: " + str(round(ro4*10,3))+' * 10^-1'+' кг/м^3')
print("Температура 3 области: " + str(round(T3/100,3)) + ' * 10^2 K')
print("Температура 4 области: " + str(round(T4/100,3)) + ' * 10^2 K')
print("Скорость звука 3 области: " + str(round(a3/100,3))+' * 10^2 м/с')
print("Скорость звука 4 области: " + str(round(a4/100,3))+' * 10^2 м/с')
print("Скорость 3 области: " + str(round(V3/1000,3)) + ' * 10^3 м/с')
print("Скорость 4 области: " + str(round(V4/1000,3)) + ' * 10^3 м/с')
print("Давление 5 области: " + str(round(p5 / 10000, 3)) + ' * 10^4' + ' Па')
print("Давление 6 области: " + str(round(p6 / 10000, 3)) + ' * 10^4' + ' Па')
print("Угол теттац5: " + str(round(tettac5,3))+ ' рад; '+ str(round((tettac5*180)/np.pi,2))+ ' град')
print("Угол теттац6: " + str(round(tettac6,3))+ ' рад; '+ str(round((tettac6*180)/np.pi,2))+ ' град')
print("Угол отклонения: " + str(round(delta,3))+ ' рад; '+ str(round((delta*180)/np.pi,2))+ ' град')
print("Угол бетта5: " + str(round(betta5,3))+ ' рад; '+ str(round((betta5*180)/np.pi,2))+ ' град')
print("Угол бетта6: " + str(round(betta6,3))+ ' рад; '+ str(round((betta6*180)/np.pi,2))+ ' град')
print("Плотность 5 области: " + str(round(ro5*10,3))+' * 10^-1'+' кг/м^3')
print("Плотность 6 области: " + str(round(ro6*10,3))+' * 10^-1'+' кг/м^3')
print("Температура 5 области: " + str(round(T5/100,3)) + ' * 10^2 K')
print("Температура 6 области: " + str(round(T6/100,3)) + ' * 10^2 K')
print("Скорость 5 области: " + str(round(V5/1000,3)) + ' * 10^3 м/с')
print("Скорость 6 области: " + str(round(V6/1000,3)) + ' * 10^3 м/с')
print("Скорость звука 5 области: " + str(round(a5/100,3))+' * 10^2 м/с')
print("Скорость звука 6 области: " + str(round(a6/100,3))+' * 10^2 м/с')
print("Число Маха 5 области: " + str(round(M5,3)))
print("Число Маха 6 области: " + str(round(M6,3)))
print("Давление торможения 5 области: " + str(round(p05/10000000,3)) + ' * 10^7'+' Па')
print("Давление торможения 6 области: " + str(round(p06/10000000,3)) + ' * 10^7'+' Па')
print("Плотность торможения 5 области: " + str(round(ro05/10,3))+' * 10^1'+' кг/м^3')
print("Плотность торможения 6 области: " + str(round(ro06/10,3))+' * 10^1'+' кг/м^3')
print("Температура торможения 5 области: " + str(round(T05/1000,3)) + ' * 10^3 K')
print("Температура торможения 6 области: " + str(round(T06/1000,3)) + ' * 10^3 K')Editor is loading...