main

 avatar
unknown
plain_text
a year ago
22 kB
3
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
from openpyxl import Workbook
import datetime


# Функции для решения
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
    Mach = closestRoot
    return Mach  # Мах
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 ReynoldsNumberAdjustment(x: float) -> float:
    return 1 - 16 * x - 412.5 * x**2 - 35000 * x**3 - 375000 * x**4
def SpecificHeat_Cp(T: float) -> float:  # Функция поиска удельной теплоемкости
    return Cp * (T/T0[0]) ** fi  # Удельная теплоемкость
def DynamicViscosity_mu(T: float) -> float:  # Функция поиска динамической вязкости
    return mu * (T/T0[0]) ** n  # Динамическая вязкость
def CoeffThermalConductivity_Lambda(T: float) -> float:  # Функция поиска коэффициента теплопроводности
    return Lambda * (T/T0[1]) ** z  # Коэффициент теплопроводности

# Параметры потока в скачках уплотнения
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(T: float, 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(T,p,ro,p1,ro1)
T2 = findFlowParameterTemperature(T,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)
my1 = np.arcsin(1/M1)
M2 = findFlowParameterMach(V2, a2)
my2 = np.arcsin(1/M2)
# Определяем параметры торможения для участка
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)
my3 = np.arcsin(1/M3)
M4 = findMachThroughFictitiousAngle(w4)
my4 = np.arcsin(1/M4)
# Поскольку поток является невязким, давление, плотность и температура торможения в областях 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 области \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

t1, t2, u = 0.3, 0.3, 0.17 # Начальные приближенные значение
tettac5,tettac6,delta = fsolve(findDeflectionAngle, (t1, t2, u))
betta5 = betta + delta
betta6 = betta - delta
# Определяем параметры потока в скачках уплотнения
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(ro4,M4,tettac6)
T5 = findFlowParameterTemperature(T3,p3,ro3,p5,ro5)
T6 = findFlowParameterTemperature(T4,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)

# \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Пограничный слой \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

fi = d.fi
n = d.n
z = d.z
T0 = d.T0
Cp = d.Cp
mu = d.mu
Lambda = d.Lambda

# Удельная теплоемкость

Cp0 = SpecificHeat_Cp(T)
Cp1 = SpecificHeat_Cp(T1)
Cp2 = SpecificHeat_Cp(T2)
Cp3 = SpecificHeat_Cp(T3)
Cp4 = SpecificHeat_Cp(T4)
Cp5 = SpecificHeat_Cp(T5)
Cp6 = SpecificHeat_Cp(T6)

mu0 = DynamicViscosity_mu(T)
mu1 = DynamicViscosity_mu(T1)
mu2 = DynamicViscosity_mu(T2)
mu3 = DynamicViscosity_mu(T3)
mu4 = DynamicViscosity_mu(T4)
mu5 = DynamicViscosity_mu(T5)
mu6 = DynamicViscosity_mu(T6)

Lambda0 = CoeffThermalConductivity_Lambda(T)
Lambda1 = CoeffThermalConductivity_Lambda(T1)
Lambda2 = CoeffThermalConductivity_Lambda(T2)
Lambda3 = CoeffThermalConductivity_Lambda(T3)
Lambda4 = CoeffThermalConductivity_Lambda(T4)
Lambda5 = CoeffThermalConductivity_Lambda(T5)
Lambda6 = CoeffThermalConductivity_Lambda(T6)

print(Cp1,Cp2,Cp3,Cp4)
print(mu1,mu2,mu3,mu4)
print(Lambda1,Lambda2,Lambda3,Lambda4)


#Вывод в эксель

wb = Workbook()
ws = wb.active
ws.append(['1','betta, град','tettac, град','M, [-]','p, 10^4 Па','p0, 10^7  Па','ro, 10^-1 кг/м^3','ro0, 10^1 кг/м^3','T, 10^2 K','T0, 10^3 K', 'V, 10^3 м/с', 'a, 10^2 м/с'])

ws['A1'] = datetime.datetime.now()
ws.append(['1',round(betta1*180/np.pi,3),round(tettac1*180/np.pi,3),round(M1,3),round(p1/10000,3),round(p01/10000000,3),round(ro1*10,3),round(ro01/10,3),round(T1/100,3),round(T01/1000,3),round(V1/1000,3),round(a1/100,3)])
ws['A2'] = datetime.datetime.now()
ws.append(['1',betta2*180/np.pi,tettac2*180/np.pi,M2,p2,p02,ro2,ro02,T2,T02,V2,a2])
ws['A3'] = datetime.datetime.now()
ws.append(['1',betta3*180/np.pi,"-",M3,p3,p03,ro3,ro03,T3,T03,V3,a3])
ws['A4'] = datetime.datetime.now()
ws.append(['1',betta4*180/np.pi,"-",M4,p4,p04,ro4,ro04,T4,T04,V4,a4])
ws['A5'] = datetime.datetime.now()
ws.append(['1',betta5*180/np.pi,tettac5*180/np.pi,M5,p5,p05,ro5,ro05,T5,T05,V5,a5])
ws['A6'] = datetime.datetime.now()
ws.append(['1',betta6*180/np.pi,tettac6*180/np.pi,M6,p6,p06,ro6,ro06,T6,T06,V6,a6])
wb.save("sample.xlsx")


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,2100))+ ' град')
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')
print("Угол мю1: " + str(round(my1,3))+ ' рад; '+ str(round((my1*180)/np.pi,2))+ ' град')
print("Угол мю2: " + str(round(my2,3))+ ' рад; '+ str(round((my2*180)/np.pi,2))+ ' град')
print("Угол мю3: " + str(round(my3,3))+ ' рад; '+ str(round((my3*180)/np.pi,2))+ ' град')
print("Угол мю4: " + str(round(my4,3))+ ' рад; '+ str(round((my4*180)/np.pi,2))+ ' град')

print([betta1*180/np.pi,tettac1*180/np.pi,M1,p1,p01,ro1,ro01,T1,T01,V1,a1])
print([betta2*180/np.pi,tettac2*180/np.pi,M2,p2,p02,ro2,ro02,T2,T02,V2,a2])
print([betta3*180/np.pi,"-",M3,p3,p03,ro3,ro03,T3,T03,V3,a3])
print([betta4*180/np.pi,"-",M4,p4,p04,ro4,ro04,T4,T04,V4,a4])
print([betta5*180/np.pi,tettac5*180/np.pi,M5,p5,p05,ro5,ro05,T5,T05,V5,a5])
print([betta6*180/np.pi,tettac6*180/np.pi,M6,p6,p06,ro6,ro06,T6,T06,V6,a6])