From 86d8317c02c6108fd3277f420b7ffe37ad84366b Mon Sep 17 00:00:00 2001 From: JNewton Date: Sat, 15 Jan 2022 09:49:00 -0600 Subject: [PATCH] Brought the work that has been done to conical flows to this branch --- CompAero/ConicalFlows.py | 230 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 230 insertions(+) create mode 100644 CompAero/ConicalFlows.py diff --git a/CompAero/ConicalFlows.py b/CompAero/ConicalFlows.py new file mode 100644 index 0000000..68c45b0 --- /dev/null +++ b/CompAero/ConicalFlows.py @@ -0,0 +1,230 @@ +from numpy import ndarray +import numpy as np +from scipy.optimize import brenth +from colorama import Back, Style, Fore +from CompAero.internal import checkValue +from CompAero.ObliqueShockRelations import ObliqueShockRelations +from GreekLetters.greekLetters import LowerCaseGreek as lcg, Misc +from CompAero.IsentropecRelations import IsentropicRelations + +from typing import List, Sequence, Union +from math import isnan, sqrt, nan, pow, radians, degrees, tan, cos, sin +from dataclasses import dataclass + + +@dataclass +class ConicalRayRatios: + mach: float + tempRay_tempInit: float + pressRay_pressInit: float + rhoRay_rhoInit: float + velRay_velInit: float + cpRay: float + rayAngle: float + + +class ConicalFlowRelations: + def __init__( + self, + gamma: float, + mach: float, + increment: float = -0.01, + coneAngle: float = nan, + shockAngle: float = nan, + inDegrees: bool = True, + maxIterations: int = 10000, + ) -> None: + self.__gamma = gamma + self.__mach = mach + self.__coneAngle = coneAngle + self.__shockAnlge = shockAngle + self.__increment = increment + self.__inDegrees = inDegrees + self.__maxIter = maxIterations + self.__velcoties = np.empty((1,)) # Store the velocites of the flow field while solving + self.__flowAngles = np.empty((1,)) # Store the angle at which the velocites were stored + + if not self.__inDegrees: + self.__coneAngle = degrees(self.__coneAngle) + self.__shockAnlge = degrees(self.__shockAnlge) + else: + self.__increment = radians(self.__increment) + + if not checkValue(self.__gamma): + err = "Concical flow cannot analysis cannot be completed with current value of gamma. Val: {}".format( + self.__gamma + ) + raise ValueError(Fore.BLACK + Back.RED + err + Style.RESET_ALL) + + if isnan(self.__increment) or self.__increment >= 0.0: + err = "Theta increment must be a negative value. Val: {}".format(self.__increment) + raise ValueError(Fore.BLACK + Back.RED + err + Style.RESET_ALL) + + if checkValue(self.__mach) and checkValue(self.__shockAnlge): + self.__coneAngle = self.__calculateConeAngle(self.__shockAnlge) + print(self.__coneAngle) + + elif checkValue(self.__mach) and checkValue(self.__coneAngle): + self.__shockAnlge = self.__calculateShockAngle() + print(self.__shockAnlge) + + else: + err = "Not enough given arguments to determine conical flow field" + raise ValueError(Fore.BLACK + Back.RED + err + Style.RESET_ALL) + + if not self.__inDegrees: + self.__coneAngle = radians(self.__coneAngle) + self.__shockAnlge = radians(self.__shockAnlge) + else: + self.__increment = degrees(self.__increment) + + @property + def gamma(self) -> float: + """ Ratio of specific heats """ + return self.__gamma + + @property + def theteIncrement(self) -> float: + """ Angle increment used when solveing Taylor-Macoll ODE""" + return self.__increment + + @property + def shockAngle(self) -> float: + """ Angle of Cone shaped shock wave""" + return self.__shockAnlge + + @property + def coneAngle(self) -> float: + """ Angle of the flow deflection angle (Cone) """ + return self.__coneAngle + + @property + def mach(self) -> float: + """ The free stream mach number seen by the cone """ + return self.__mach + + @property + def shockJumpConditions(self) -> ObliqueShockRelations: + return ObliqueShockRelations( + gamma=self.__gamma, mach=self.__mach, shockAngle=self.__shockAnlge, degrees=self.__inDegrees + ) + + def __OdeFunction1(self, theta: float, Vr: float, Vtheta: float, Vr2: float, Vtheta2: float) -> float: + return (self.__gamma - 1) / 2 * (1 - Vr2 - Vtheta2) * (2 * Vr + Vtheta / tan(theta)) - Vtheta2 * Vr + + def __OdeFunction2(self, theta: float, Vr: float, Vtheta: float) -> float: + Vtheta2 = pow(Vtheta, 2) + Vr2 = pow(Vr, 2) + return ( + -1 + * self.__OdeFunction1(theta, Vr, Vtheta, Vr2, Vtheta2) + / ((self.__gamma - 1) / 2 * (1 - Vr2 - Vtheta2) - Vtheta2) + ) + + def __initializeVelocites(self, angle: float) -> Sequence[float]: + obliqueShockFlow = ObliqueShockRelations(self.__gamma, mach=self.__mach, shockAngle=angle) + velocity = 1 / sqrt(2 / ((self.__gamma - 1) * pow(obliqueShockFlow.mach2, 2)) + 1) + Vr = velocity * cos(radians(angle - obliqueShockFlow.wedgeAngle)) + Vtheta = -1 * velocity * sin(radians(angle - obliqueShockFlow.wedgeAngle)) + return [Vr, Vtheta] + + def __calculateConeAngle( + self, shockAngle: float, offset: float = 0.0, anglesToSave: ndarray = np.array([]) + ) -> float: + """Calculates Cone angle given shock angle and upstream mach number""" + vr, vtheta = self.__initializeVelocites(shockAngle) + shockAngle = radians(shockAngle) + + coneAngle = shockAngle + angleSaveIter = 0 + + currentIter = 0 + while vtheta < 0.0 and currentIter <= self.__maxIter: + theta2 = coneAngle + self.__increment / 2 + k11 = self.__increment * vtheta + k12 = self.__increment * self.__OdeFunction2(shockAngle, vr, vtheta) + + k21 = self.__increment * (vtheta + k12 / 2) + k22 = self.__increment * self.__OdeFunction2(theta2, (vr + k11 / 2), (vtheta + k12 / 2)) + + k31 = self.__increment * (vtheta + k22 / 2) + k32 = self.__increment * self.__OdeFunction2(theta2, (vr + k21 / 2), (vtheta + k22 / 2)) + + coneAngle += self.__increment + k41 = self.__increment * (vtheta + k32) + k42 = self.__increment * self.__OdeFunction2(coneAngle, (vr + k31), (vtheta + k32)) + + vr += (k11 + 2 * k21 + 2 * k31 + k41) / 6 + vtheta += (k12 + 2 * k22 + 2 * k32 + k42) / 6 + + currentIter += 1 + + if anglesToSave.size > 0: + if ( + abs(np.min(np.abs(anglesToSave - round(coneAngle, 4)))) <= 1e-4 + and angleSaveIter < self.__velcoties.shape[0] + ): + self.__velcoties[angleSaveIter] = sqrt(pow(vr, 2) + pow(vtheta, 2)) + self.__flowAngles[angleSaveIter] = coneAngle + angleSaveIter += 1 + + return degrees(coneAngle) - offset + + def __calculateShockAngle(self) -> float: + machWaveAngle = degrees(ObliqueShockRelations.calcMachWaveAngle(self.__mach)) + return brenth(self.__calculateConeAngle, machWaveAngle, 80, args=(self.__coneAngle,)) + + def calculateConeFlowParameters( + self, tempInit: float, gasConstR: float, angles: Union[float, ndarray, list], inDegrees: bool = True + ) -> Union[ConicalRayRatios, dict]: + if isinstance(angles, np.ndarray): + self.__velcoties = np.zeros(angles.shape) + self.__flowAngles = np.zeros(angles.shape) + + elif isinstance(angles, float) or isinstance(angles, int): + self.__velcoties = np.array([0]) + self.__flowAngles = np.array([0]) + angles = np.array([angles]) + + elif isinstance(angles, List): + self.__velcoties = np.zeros((len(angles),)) + self.__flowAngles = np.zeros((len(angles),)) + angles = np.array(angles) + + else: + err = "Parameter of Type {} not supported!".format(type(angles)) + raise ValueError(Fore.BLACK + Back.RED + err + Style.RESET_ALL) + + if inDegrees: + angles = np.radians(angles) + + self.__increment = radians(self.__increment) + self.__calculateConeAngle(self.__shockAnlge, anglesToSave=angles) + + flowParams = {} + aInit = sqrt(self.__gamma * gasConstR * tempInit) + shockJumpParams = ObliqueShockRelations( + self.__gamma, mach=self.__mach, shockAngle=self.__shockAnlge, useDegrees=self.__inDegrees + ) + conditionsInit = IsentropicRelations(self.__gamma, mach=self.__mach) + + for velocity, angle in zip(self.__velcoties, self.__flowAngles): + params = ConicalRayRatios + vel = sqrt(2 * self.__gamma * gasConstR / (self.__gamma - 1) * tempInit) * velocity + mach = sqrt( + pow(vel, 2) / (self.__gamma * gasConstR * tempInit - (self.__gamma - 1) / 2 * pow(vel, 2)) + ) + params.mach = mach + isentropRay = IsentropicRelations(self.__gamma, mach) + + # Temperatures + params.tempRay_tempInit = conditionsInit.t0_t / isentropRay.t0_t + params.pressRay_pressInit = conditionsInit.p0_p * shockJumpParams.po2_po1 / isentropRay.p0_p + params.rayAngle = degrees(angle) + params.rhoRay_rhoInit = params.pressRay_pressInit / params.tempRay_tempInit + params.cpRay = 2 / self.__gamma / pow(self.__mach, 2) * (params.pressRay_pressInit - 1) + params.velRay_velInit = velocity / aInit / self.__mach + flowParams[params.rayAngle] = params + + return flowParams +