Skip to content

Commit

Permalink
BUG: remove airfoil implementation
Browse files Browse the repository at this point in the history 
The airfoil implementation from pull request #47 was removed and
will be readded after careful review and reimplementation. The
current implementation is buggy at best. This commit fixes
issue #115.
  • Loading branch information
@Lucas-KB
Lucas-KB committed Feb 3, 2022
1 parent 13fdc8e commit 1e36730
Show file tree
Hide file tree
Showing 2 changed files with 29 additions and 131 deletions.
9 changes: 7 additions & 2 deletions rocketpy/Flight.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1298,6 +1298,7 @@ def uDot(self, t, u, postProcessing=False):
# Calculate lift and moment for each component of the rocket
for aerodynamicSurface in self.rocket.aerodynamicSurfaces:
compCp = aerodynamicSurface[0][2]
clalpha = aerodynamicSurface[1]
# Component absolute velocity in body frame
compVxB = vxB + compCp * omega2
compVyB = vyB - compCp * omega1
Expand All @@ -1324,10 +1325,14 @@ def uDot(self, t, u, postProcessing=False):
compStreamVzBn = compStreamVzB / compStreamSpeed
if -1 * compStreamVzBn < 1:
compAttackAngle = np.arccos(-compStreamVzBn)
cLift = abs(aerodynamicSurface[1](compAttackAngle))
# Component lift force magnitude
compLift = (
0.5 * rho * (compStreamSpeed ** 2) * self.rocket.area * cLift
0.5
* rho
* (compStreamSpeed ** 2)
* self.rocket.area
* clalpha
* compAttackAngle
)
# Component lift force components
liftDirNorm = (compStreamVxB ** 2 + compStreamVyB ** 2) ** 0.5
Expand Down
151 changes: 22 additions & 129 deletions rocketpy/Rocket.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -370,13 +370,8 @@ def evaluateStaticMargin(self):
# Calculate total lift coefficient derivative and center of pressure
if len(self.aerodynamicSurfaces) > 0:
for aerodynamicSurface in self.aerodynamicSurfaces:
self.totalLiftCoeffDer += aerodynamicSurface[1].differentiate(
x=1e-2, dx=1e-3
)
self.cpPosition += (
aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3)
* aerodynamicSurface[0][2]
)
self.totalLiftCoeffDer += aerodynamicSurface[1]
self.cpPosition += aerodynamicSurface[1] * aerodynamicSurface[0][2]
self.cpPosition /= self.totalLiftCoeffDer

# Calculate static margin
Expand Down Expand Up @@ -416,8 +411,6 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM):
Returns
-------
cldata : Function
Object of the Function class. Contains tail's lift data.
self : Rocket
Object of the Rocket class.
"""
Expand All @@ -435,16 +428,9 @@ def addTail(self, topRadius, bottomRadius, length, distanceToCM):

# Calculate clalpha
clalpha = -2 * (1 - r ** (-2)) * (topRadius / rref) ** 2
cldata = Function(
lambda x: clalpha * x,
"Alpha (rad)",
"Cl",
interpolation="linear",
extrapolation="natural",
)

# Store values as new aerodynamic surface
tail = [(0, 0, cpz), cldata, "Tail"]
tail = [(0, 0, cpz), clalpha, "Tail"]
self.aerodynamicSurfaces.append(tail)

# Refresh static margin calculation
Expand Down Expand Up @@ -476,8 +462,6 @@ def addNose(self, length, kind, distanceToCM):
Returns
-------
cldata : Function
Object of the Function class. Contains nose's lift data.
self : Rocket
Object of the Rocket class.
"""
Expand All @@ -499,16 +483,9 @@ def addNose(self, length, kind, distanceToCM):

# Calculate clalpha
clalpha = 2
cldata = Function(
lambda x: clalpha * x,
"Alpha (rad)",
"Cl",
interpolation="linear",
extrapolation="natural",
)

# Store values
nose = [(0, 0, cpz), cldata, "Nose Cone"]
nose = [(0, 0, cpz), clalpha, "Nose Cone"]
self.aerodynamicSurfaces.append(nose)

# Refresh static margin calculation
Expand All @@ -526,7 +503,6 @@ def addFins(
distanceToCM,
radius=0,
cantAngle=0,
airfoil=None,
):
"""Create a fin set, storing its parameters as part of the
aerodynamicSurfaces list. Its parameters are the axial position
Expand Down Expand Up @@ -556,16 +532,9 @@ def addFins(
cantAngle : int, float, optional
Fins cant angle with respect to the rocket centerline. Must
be given in degrees.
airfoil : string
Fin's lift curve. It must be a .csv file. The .csv file shall
contain no headers and the first column must specify time in
seconds, while the second column specifies lift coefficient. Lift
coefficient is dimensionaless.
Returns
-------
cldata : Function
Object of the Function class. Contains fin's lift data.
self : Rocket
Object of the Rocket class.
"""
Expand Down Expand Up @@ -605,102 +574,26 @@ def addFins(
+ (1 / 6) * (Cr + Ct - Cr * Ct / (Cr + Ct))
)

# Calculate lift parameters for planar fins
if not airfoil:
# Calculate clalpha
clalpha = (4 * n * (s / d) ** 2) / (1 + np.sqrt(1 + (2 * Lf / Yr) ** 2))
clalpha *= 1 + radius / (s + radius)

# # Create a function of lift values by attack angle
cldata = Function(
lambda x: clalpha * x, "Alpha (rad)", "Cl", interpolation="linear"
)
# Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
clfDelta = n * (Ymac + radius) * clalpha / d
cldOmega = (
n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
)
rollParameters = (
[clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
)

# Store values
fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
self.aerodynamicSurfaces.append(fin)

# Refresh static margin calculation
self.evaluateStaticMargin()

# Return self
return self.aerodynamicSurfaces[-1]

else:

def cnalfa1(cn):
"""Calculates the normal force coefficient derivative of a 3D
airfoil for a given Cnalfa0
Parameters
----------
cn : int
Normal force coefficient derivative of a 2D airfoil.
Returns
-------
Cnalfa1 : int
Normal force coefficient derivative of a 3D airfoil.
"""

# Retrieve parameters for calculations
Af = (Cr + Ct) * span / 2
# fin area
AR = 2 * (span ** 2) / Af # Aspect ratio
gamac = np.arctan((Cr - Ct) / (2 * span))
# mid chord angle
FD = 2 * np.pi * AR / (cn * np.cos(gamac))
Cnalfa1 = (
cn
* FD
* (Af / self.area)
* np.cos(gamac)
/ (2 + FD * (1 + (4 / FD ** 2)) ** 0.5)
)
return Cnalfa1

# Import the lift curve as a function of lift values by attack angle
read = genfromtxt(airfoil, delimiter=",")

# Applies number of fins to lift coefficient data
data = [[cl[0], (n / 2) * cnalfa1(cl[1])] for cl in read]
cldata = Function(
data,
"Alpha (rad)",
"Cl",
interpolation="linear",
extrapolation="natural",
)

# Takes an approximation to an angular coefficient
clalpha = cldata.differentiate(x=0, dx=1e-2)

# Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
clfDelta = n * (Ymac + radius) * clalpha / d
cldOmega = (
n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
)
rollParameters = (
[clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
)
# Calculate clalpha
clalpha = (4 * n * (s / d) ** 2) / (1 + np.sqrt(1 + (2 * Lf / Yr) ** 2))
clalpha *= 1 + radius / (s + radius)

# Parameters for Roll Moment. Documented at: https://github.com/Projeto-Jupiter/RocketPy/blob/develop/docs/technical/aerodynamics/Roll_Equations.pdf
clfDelta = n * (Ymac + radius) * clalpha / d
cldOmega = n * clalpha * np.cos(cantAngleRad) * trapezoidalConstant / (Af * d)
rollParameters = (
[clfDelta, cldOmega, cantAngleRad] if cantAngleRad != 0 else [0, 0, 0]
)

# Store values
fin = [(0, 0, cpz), cldata, rollParameters, "Fins"]
self.aerodynamicSurfaces.append(fin)
# Store values
fin = [(0, 0, cpz), clalpha, rollParameters, "Fins"]
self.aerodynamicSurfaces.append(fin)

# Refresh static margin calculation
self.evaluateStaticMargin()
# Refresh static margin calculation
self.evaluateStaticMargin()

# Return self
return self.aerodynamicSurfaces[-1]
# Return self
return self.aerodynamicSurfaces[-1]

def addParachute(
self, name, CdS, trigger, samplingRate=100, lag=0, noise=(0, 0, 0)
Expand Down Expand Up @@ -987,7 +880,7 @@ def allInfo(self):
print("\nAerodynamics Lift Coefficient Derivatives")
for aerodynamicSurface in self.aerodynamicSurfaces:
name = aerodynamicSurface[-1]
clalpha = aerodynamicSurface[1].differentiate(x=1e-2, dx=1e-3)
clalpha = aerodynamicSurface[1]
print(
name + " Lift Coefficient Derivative: {:.3f}".format(clalpha) + "/rad"
)
Expand Down

0 comments on commit 1e36730

Please sign in to comment.