Untitled

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 sympy import diff, symbols, sin, tan

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 LongitudinalAndNormalForces(pi: float) -> float:  # Функция поиска продольной и нормальной сил
    X = (pi - p) * (c / 2)
    Y = (pi - p) * (b / 2)
    return X, Y # Продольная и нормальная силы

# ПС
def CriticalReynoldsNumber(Tr, M):
    Tst_otn = Tst / Tr  # Относительная температура стенки
    x = (Tst_otn - 1) / M ** 2
    y = 1 - 16 * x - 412.5 * x**2 - 35000 * x**3 - 375000 * x**4
    return criticalReynoldsNumber * y
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  findPressureCoefficient (pi: float) -> float: # Функция поиска коэффициента давления
    return (pi - p) / q # Коэффициент давления
def findDefiningTemperature (lam_or_turb: float,T: float, M: float) -> float: # Функция поика температуры восстановления и определяющей температуры
    rl = 0.83
    rt = 0.88
    r = rl
    if (lam_or_turb == 1):
        r = rt
    Tr = T * (1 + r * (k - 1) * 0.5 * M ** 2)
    T_z = 0.5 * (Tst + T) + 0.22 * (Tr - T)
    T_z_previous = 0
    while abs(T_z - T_z_previous) > (T_z * 0.001):
        T_z_previous = T_z
        Cpi = SpecificHeat_Cp(T_z)
        mui = DynamicViscosity_mu(T_z)
        Lambdai = CoeffThermalConductivity_Lambda(T_z)
        Prl = (Cpi * mui) / Lambdai
        r = np.sqrt(Prl)
        if (lam_or_turb == 0):
            r = Prl ** (1/3)
    return Tr,T_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 findFlowParameterSpeedOfSound(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  # Температура торможения

# Производная Vθ по θ
def findFlowParameterАcceleration(Vr: float, Vθ: float) -> float:
    return ((k + 1) / 2) * a_kr ** 2 - ((k - 1) / 2) * (Vr ** 2 - Vθ ** 2)
def dVθ_na_dθ(Vr: float, Vθ: float, θ: float) -> float:
    return (-Vθ * (np.tan(θ) ** (-1)) + Vr * ((Vθ ** 2)/(findFlowParameterАcceleration(Vr, Vθ)) - 2)) / (1 - (Vθ ** 2)/(findFlowParameterАcceleration(Vr, Vθ)))

# Касательная и нормальная скорости
def Vr(θ: float) -> float:
    return V * np.cos(θ)
def Vθ(Vr: float,θc: float) -> float:
    return (-1) * Vr * (2 + (k - 1) * M ** 2 * (np.sin(θc) ** 2))/((k + 1) * M ** 2 * (np.sin(θc)) ** 2)

# Получение данных из data.py
variant = d.variant
dm = d.dm
P = d.P
βk = np.radians(d.bettak)
print("Угол полураствора конуса: " + str(round(βk,3))+ ' рад; '+ str(round((βk*180)/np.pi,2))+ ' град')
lk = d.lk
Xcm = d.Xcm
l_s_chertoy_k = d.l_s_chertoy_k
h = d.h
M = d.M_inf
R = d.R

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Определение параметров стандартной атмосферы \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

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  # Скорость потока
atm.vivod()

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Плоский СУ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
# alfa = 0
k = 1.4
θc_plosk = findShockWaveAngle(βk)
print("Угол плоского скачка уплотнения: " + str(round(θc_plosk,3))+ ' рад; '+ str(round((θc_plosk*180)/np.pi,2))+ ' град')

#\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ Метод последовательных приближений \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

T0 = T * (1 + ((k - 1) * M ** 2) / 2)
a_kr = np.sqrt((2 * k * R * T0) / (k + 1))

coef = 0.85
θc = coef * θc_plosk
print("Угол плоского скачка уплотнения умноженый на 0.85: " + str(round(θc,3))+ ' рад; '+ str(round((θc*180)/np.pi,2))+ ' град')
θkonusa = θc
βn = 60
dθ = 0.0001
dθn = 0.000001
lol, lol2 = 0, 0
while abs((βn-βk)/βk) > 0.001:
    θc -= dθ
    Vr_c = Vr(θc)
    Vθ_c = Vθ(Vr_c,θc)
    Vr_next = Vr_c
    Vθ_next = Vθ_c
    θ = θc
    while Vθ_next < 0:
        Vr_c = Vr_next
        Vθ_c = Vθ_next
        Vr_next = Vr_c - Vθ_c * dθ
        Vθ_next = Vθ_c - dVθ_na_dθ(Vr_c, Vθ_c, θ) * dθ
        θ -= dθn
    βn = θc
    if (lol>1) and (lol2<1):
        print('РАСЧЁТ: 1%')
        lol2+=1
    if (lol>130) and (lol2<2):
        print('РАСЧЁТ: 10%')
        lol2+=1
    if (lol>260) and (lol2<3):
        print('РАСЧЁТ: 20%')
        lol2+=1
    if (lol>390) and (lol2<4):
        print('РАСЧЁТ: 30%')
        lol2+=1
    if (lol>520) and (lol2<5):
        print('РАСЧЁТ: 40%')
        lol2+=1
    if (lol>650) and (lol2<6):
        print('РАСЧЁТ: 50%')
        lol2+=1
    if (lol>1200) and (lol2<7):
        print('РАСЧЁТ: 95%')
        lol2+=1
    lol += 1

print(abs((βn-βk)/βk))
print("Расчетный угол СУ конуса: " + str(round(θkonusa,3))+ ' рад; '+ str(round((θkonusa*180)/np.pi,2))+ ' град')
print("Расчетный угол конуса: " + str(round(βn,3))+ ' рад; '+ str(round((βn*180)/np.pi,2))+ ' град')
print("Скорость Vr_с: " + str(round(Vr_c,3)) + ' м/с')
print("Скорость Vr_next: " + str(round(Vr_next,3)) + ' м/с')
print("Скорость Vθ_с: " + str(round(Vθ_c,3)) + ' м/с')
print("Скорость Vθ_next: " + str(round(Vθ_next,3)) + ' м/с')
print(lol)