Skip to content
This repository was archived by the owner on Oct 30, 2022. It is now read-only.

Issue #12: team-01 - Improving current fin flutter method #105

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Issue #12: team-01 - Improving current fin flutter method #105

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
56e2d6b
working on optimizing fin flutter velocity calculation using Apogee N…
aasitvora99 Aug 28, 2022
30f1c7a
changed flutter velocity logic for fin flutter simulation: more detai…
aasitvora99 Aug 28, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/notebooks/dispersion_analysis/RocketPy
Open in desktop
Submodule RocketPy added at 256f90
1,439 changes: 0 additions & 1,439 deletions docs/notebooks/getting_started.ipynb
View file
Open in desktop

This file was deleted.

1,037 changes: 1,037 additions & 0 deletions getting_started.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

123 changes: 64 additions & 59 deletions rocketpy/Flight.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -836,7 +836,7 @@ def __init__(
self.y[0] ** 2
+ self.y[1] ** 2
+ (self.y[2] - self.env.elevation) ** 2
>= self.effective1RL**2
>= self.effective1RL ** 2
):
# Rocket is out of rail
# Check exactly when it went out using root finding
Expand All @@ -850,7 +850,7 @@ def __init__(
# Get points
y0 = (
sum([self.solution[-2][i] ** 2 for i in [1, 2, 3]])
- self.effective1RL**2
- self.effective1RL ** 2
)
yp0 = 2 * sum(
[
Expand All @@ -861,7 +861,7 @@ def __init__(
t1 = self.solution[-1][0] - self.solution[-2][0]
y1 = (
sum([self.solution[-1][i] ** 2 for i in [1, 2, 3]])
- self.effective1RL**2
- self.effective1RL ** 2
)
yp1 = 2 * sum(
[
Expand All @@ -876,15 +876,15 @@ def __init__(
D = float(phase.solver.step_size)
d = float(y0)
c = float(yp0)
b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2))
a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3)) + 1e-5
b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D ** 2))
a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D ** 3)) + 1e-5
# Find roots
d0 = b**2 - 3 * a * c
d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2
c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3)
d0 = b ** 2 - 3 * a * c
d1 = 2 * b ** 3 - 9 * a * b * c + 27 * d * a ** 2
c1 = ((d1 + (d1 ** 2 - 4 * d0 ** 3) ** (0.5)) / 2) ** (1 / 3)
t_roots = []
for k in [0, 1, 2]:
c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k
c2 = c1 * (-1 / 2 + 1j * (3 ** 0.5) / 2) ** k
t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2))
# Find correct root
valid_t_root = []
Expand Down Expand Up @@ -971,15 +971,15 @@ def __init__(
D = float(phase.solver.step_size)
d = float(y0)
c = float(yp0)
b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D**2))
a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D**3))
b = float((3 * y1 - yp1 * D - 2 * c * D - 3 * d) / (D ** 2))
a = float(-(2 * y1 - yp1 * D - c * D - 2 * d) / (D ** 3))
# Find roots
d0 = b**2 - 3 * a * c
d1 = 2 * b**3 - 9 * a * b * c + 27 * d * a**2
c1 = ((d1 + (d1**2 - 4 * d0**3) ** (0.5)) / 2) ** (1 / 3)
d0 = b ** 2 - 3 * a * c
d1 = 2 * b ** 3 - 9 * a * b * c + 27 * d * a ** 2
c1 = ((d1 + (d1 ** 2 - 4 * d0 ** 3) ** (0.5)) / 2) ** (1 / 3)
t_roots = []
for k in [0, 1, 2]:
c2 = c1 * (-1 / 2 + 1j * (3**0.5) / 2) ** k
c2 = c1 * (-1 / 2 + 1j * (3 ** 0.5) / 2) ** k
t_roots.append(-(1 / (3 * a)) * (b + c2 + d0 / c2))
# Find correct root
valid_t_root = []
Expand Down Expand Up @@ -1237,14 +1237,14 @@ def uDotRail1(self, t, u, postProcessing=False):
# Calculate Forces
Thrust = self.rocket.motor.thrust.getValueOpt(t)
rho = self.env.density.getValueOpt(z)
R3 = -0.5 * rho * (freestreamSpeed**2) * self.rocket.area * (dragCoeff)
R3 = -0.5 * rho * (freestreamSpeed ** 2) * self.rocket.area * (dragCoeff)

# Calculate Linear acceleration
a3 = (R3 + Thrust) / M - (e0**2 - e1**2 - e2**2 + e3**2) * self.env.g
a3 = (R3 + Thrust) / M - (e0 ** 2 - e1 ** 2 - e2 ** 2 + e3 ** 2) * self.env.g
if a3 > 0:
ax = 2 * (e1 * e3 + e0 * e2) * a3
ay = 2 * (e2 * e3 - e0 * e1) * a3
az = (1 - 2 * (e1**2 + e2**2)) * a3
az = (1 - 2 * (e1 ** 2 + e2 ** 2)) * a3
else:
ax, ay, az = 0, 0, 0

Expand Down Expand Up @@ -1343,15 +1343,15 @@ def uDot(self, t, u, postProcessing=False):
a = b * Mt / M
rN = self.rocket.motor.nozzleRadius
# Prepare transformation matrix
a11 = 1 - 2 * (e2**2 + e3**2)
a11 = 1 - 2 * (e2 ** 2 + e3 ** 2)
a12 = 2 * (e1 * e2 - e0 * e3)
a13 = 2 * (e1 * e3 + e0 * e2)
a21 = 2 * (e1 * e2 + e0 * e3)
a22 = 1 - 2 * (e1**2 + e3**2)
a22 = 1 - 2 * (e1 ** 2 + e3 ** 2)
a23 = 2 * (e2 * e3 - e0 * e1)
a31 = 2 * (e1 * e3 - e0 * e2)
a32 = 2 * (e2 * e3 + e0 * e1)
a33 = 1 - 2 * (e1**2 + e2**2)
a33 = 1 - 2 * (e1 ** 2 + e2 ** 2)
# Transformation matrix: (123) -> (XYZ)
K = [[a11, a12, a13], [a21, a22, a23], [a31, a32, a33]]
# Transformation matrix: (XYZ) -> (123) or K transpose
Expand All @@ -1373,7 +1373,7 @@ def uDot(self, t, u, postProcessing=False):
else:
dragCoeff = self.rocket.powerOffDrag.getValueOpt(freestreamMach)
rho = self.env.density.getValueOpt(z)
R3 = -0.5 * rho * (freestreamSpeed**2) * self.rocket.area * (dragCoeff)
R3 = -0.5 * rho * (freestreamSpeed ** 2) * self.rocket.area * (dragCoeff)
# Off center moment
M1 += self.rocket.cpEccentricityY * R3
M2 -= self.rocket.cpEccentricityX * R3
Expand All @@ -1400,23 +1400,23 @@ def uDot(self, t, u, postProcessing=False):
compStreamVyB = compWindVyB - compVyB
compStreamVzB = compWindVzB - compVzB
compStreamSpeed = (
compStreamVxB**2 + compStreamVyB**2 + compStreamVzB**2
compStreamVxB ** 2 + compStreamVyB ** 2 + compStreamVzB ** 2
) ** 0.5
# Component attack angle and lift force
compAttackAngle = 0
compLift, compLiftXB, compLiftYB = 0, 0, 0
if compStreamVxB**2 + compStreamVyB**2 != 0:
if compStreamVxB ** 2 + compStreamVyB ** 2 != 0:
# Normalize component stream velocity in body frame
compStreamVzBn = compStreamVzB / compStreamSpeed
if -1 * compStreamVzBn < 1:
compAttackAngle = np.arccos(-compStreamVzBn)
cLift = aerodynamicSurface["cl"](compAttackAngle, freestreamMach)
# Component lift force magnitude
compLift = (
0.5 * rho * (compStreamSpeed**2) * self.rocket.area * cLift
0.5 * rho * (compStreamSpeed ** 2) * self.rocket.area * cLift
)
# Component lift force components
liftDirNorm = (compStreamVxB**2 + compStreamVyB**2) ** 0.5
liftDirNorm = (compStreamVxB ** 2 + compStreamVyB ** 2) ** 0.5
compLiftXB = compLift * (compStreamVxB / liftDirNorm)
compLiftYB = compLift * (compStreamVyB / liftDirNorm)
# Add to total lift force
Expand All @@ -1429,7 +1429,7 @@ def uDot(self, t, u, postProcessing=False):
if aerodynamicSurface["name"] == "Fins":
Clfdelta, Cldomega, cantAngleRad = aerodynamicSurface["roll parameters"]
M3f = (
(1 / 2 * rho * freestreamSpeed**2)
(1 / 2 * rho * freestreamSpeed ** 2)
* self.rocket.area
* 2
* self.rocket.radius
Expand All @@ -1450,26 +1450,26 @@ def uDot(self, t, u, postProcessing=False):
alpha1 = (
M1
- (
omega2 * omega3 * (Rz + Tz - Ri - Ti - mu * b**2)
omega2 * omega3 * (Rz + Tz - Ri - Ti - mu * b ** 2)
+ omega1
* (
(TiDot + MtDot * (Mr - 1) * (b / M) ** 2)
- MtDot * ((rN / 2) ** 2 + (c - b * mu / Mr) ** 2)
)
)
) / (Ri + Ti + mu * b**2)
) / (Ri + Ti + mu * b ** 2)
alpha2 = (
M2
- (
omega1 * omega3 * (Ri + Ti + mu * b**2 - Rz - Tz)
omega1 * omega3 * (Ri + Ti + mu * b ** 2 - Rz - Tz)
+ omega2
* (
(TiDot + MtDot * (Mr - 1) * (b / M) ** 2)
- MtDot * ((rN / 2) ** 2 + (c - b * mu / Mr) ** 2)
)
)
) / (Ri + Ti + mu * b**2)
alpha3 = (M3 - omega3 * (TzDot - MtDot * (rN**2) / 2)) / (Rz + Tz)
) / (Ri + Ti + mu * b ** 2)
alpha3 = (M3 - omega3 * (TzDot - MtDot * (rN ** 2) / 2)) / (Rz + Tz)
# Euler parameters derivative
e0Dot = 0.5 * (-omega1 * e1 - omega2 * e2 - omega3 * e3)
e1Dot = 0.5 * (omega1 * e0 + omega3 * e2 - omega2 * e3)
Expand All @@ -1478,7 +1478,7 @@ def uDot(self, t, u, postProcessing=False):

# Linear acceleration
L = [
(R1 - b * Mt * (omega2**2 + omega3**2) - 2 * c * MtDot * omega2) / M,
(R1 - b * Mt * (omega2 ** 2 + omega3 ** 2) - 2 * c * MtDot * omega2) / M,
(R2 + b * Mt * (alpha3 + omega1 * omega2) + 2 * c * MtDot * omega1) / M,
(R3 - b * Mt * (alpha2 - omega1 * omega3) + Thrust) / M,
]
Expand Down Expand Up @@ -1549,11 +1549,11 @@ def uDotParachute(self, t, u, postProcessing=False):
R = 1.5
rho = self.env.density.getValueOpt(u[2])
to = 1.2
ma = ka * rho * (4 / 3) * np.pi * R**3
ma = ka * rho * (4 / 3) * np.pi * R ** 3
mp = self.rocket.mass
eta = 1
Rdot = (6 * R * (1 - eta) / (1.2**6)) * (
(1 - eta) * t**5 + eta * (to**3) * (t**2)
Rdot = (6 * R * (1 - eta) / (1.2 ** 6)) * (
(1 - eta) * t ** 5 + eta * (to ** 3) * (t ** 2)
)
Rdot = 0
# Get relevant state data
Expand All @@ -1569,7 +1569,7 @@ def uDotParachute(self, t, u, postProcessing=False):
freestreamZ = vz
# Determine drag force
pseudoD = (
-0.5 * rho * CdS * freestreamSpeed - ka * rho * 4 * np.pi * (R**2) * Rdot
-0.5 * rho * CdS * freestreamSpeed - ka * rho * 4 * np.pi * (R ** 2) * Rdot
)
Dx = pseudoD * freestreamX
Dy = pseudoD * freestreamY
Expand Down Expand Up @@ -1749,19 +1749,19 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):

# Kinematics functions and values
# Velocity Magnitude
self.speed = (self.vx**2 + self.vy**2 + self.vz**2) ** 0.5
self.speed = (self.vx ** 2 + self.vy ** 2 + self.vz ** 2) ** 0.5
self.speed.setOutputs("Speed - Velocity Magnitude (m/s)")
maxSpeedTimeIndex = np.argmax(self.speed[:, 1])
self.maxSpeed = self.speed[maxSpeedTimeIndex, 1]
self.maxSpeedTime = self.speed[maxSpeedTimeIndex, 0]
# Acceleration
self.acceleration = (self.ax**2 + self.ay**2 + self.az**2) ** 0.5
self.acceleration = (self.ax ** 2 + self.ay ** 2 + self.az ** 2) ** 0.5
self.acceleration.setOutputs("Acceleration Magnitude (m/s²)")
maxAccelerationTimeIndex = np.argmax(self.acceleration[:, 1])
self.maxAcceleration = self.acceleration[maxAccelerationTimeIndex, 1]
self.maxAccelerationTime = self.acceleration[maxAccelerationTimeIndex, 0]
# Path Angle
self.horizontalSpeed = (self.vx**2 + self.vy**2) ** 0.5
self.horizontalSpeed = (self.vx ** 2 + self.vy ** 2) ** 0.5
pathAngle = (180 / np.pi) * np.arctan2(
self.vz[:, 1], self.horizontalSpeed[:, 1]
)
Expand All @@ -1770,9 +1770,9 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
# Attitude Angle
self.attitudeVectorX = 2 * (self.e1 * self.e3 + self.e0 * self.e2) # a13
self.attitudeVectorY = 2 * (self.e2 * self.e3 - self.e0 * self.e1) # a23
self.attitudeVectorZ = 1 - 2 * (self.e1**2 + self.e2**2) # a33
self.attitudeVectorZ = 1 - 2 * (self.e1 ** 2 + self.e2 ** 2) # a33
horizontalAttitudeProj = (
self.attitudeVectorX**2 + self.attitudeVectorY**2
self.attitudeVectorX ** 2 + self.attitudeVectorY ** 2
) ** 0.5
attitudeAngle = (180 / np.pi) * np.arctan2(
self.attitudeVectorZ[:, 1], horizontalAttitudeProj[:, 1]
Expand All @@ -1793,9 +1793,9 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
attitudeLateralPlaneProjY = self.attitudeVectorY[:, 1] - attitudeLateralProjY
attitudeLateralPlaneProjZ = self.attitudeVectorZ[:, 1]
attitudeLateralPlaneProj = (
attitudeLateralPlaneProjX**2
+ attitudeLateralPlaneProjY**2
+ attitudeLateralPlaneProjZ**2
attitudeLateralPlaneProjX ** 2
+ attitudeLateralPlaneProjY ** 2
+ attitudeLateralPlaneProjZ ** 2
) ** 0.5
lateralAttitudeAngle = (180 / np.pi) * np.arctan2(
attitudeLateralProj, attitudeLateralPlaneProj
Expand Down Expand Up @@ -1887,11 +1887,11 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
self.railButton2ShearForce[:outOfRailTimeIndex]
)
# Aerodynamic Lift and Drag
self.aerodynamicLift = (self.R1**2 + self.R2**2) ** 0.5
self.aerodynamicLift = (self.R1 ** 2 + self.R2 ** 2) ** 0.5
self.aerodynamicLift.setOutputs("Aerodynamic Lift Force (N)")
self.aerodynamicDrag = -1 * self.R3
self.aerodynamicDrag.setOutputs("Aerodynamic Drag Force (N)")
self.aerodynamicBendingMoment = (self.M1**2 + self.M2**2) ** 0.5
self.aerodynamicBendingMoment = (self.M1 ** 2 + self.M2 ** 2) ** 0.5
self.aerodynamicBendingMoment.setOutputs("Aerodynamic Bending Moment (N m)")
self.aerodynamicSpinMoment = self.M3
self.aerodynamicSpinMoment.setOutputs("Aerodynamic Spin Moment (N m)")
Expand All @@ -1903,7 +1903,7 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
Ri = self.rocket.inertiaI
Tz = self.rocket.motor.inertiaZ
Ti = self.rocket.motor.inertiaI
I1, I2, I3 = (Ri + Ti + mu * b**2), (Ri + Ti + mu * b**2), (Rz + Tz)
I1, I2, I3 = (Ri + Ti + mu * b ** 2), (Ri + Ti + mu * b ** 2), (Rz + Tz)
# Redefine I1, I2 and I3 grid
grid = self.vx[:, 0]
I1 = Function(np.column_stack([grid, I1(grid)]), "Time (s)")
Expand All @@ -1919,9 +1919,9 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
vx, vy, vz = self.vx, self.vy, self.vz
w1, w2, w3 = self.w1, self.w2, self.w3
# Kinetic Energy
self.rotationalEnergy = 0.5 * (I1 * w1**2 + I2 * w2**2 + I3 * w3**2)
self.rotationalEnergy = 0.5 * (I1 * w1 ** 2 + I2 * w2 ** 2 + I3 * w3 ** 2)
self.rotationalEnergy.setOutputs("Rotational Kinetic Energy (J)")
self.translationalEnergy = 0.5 * totalMass * (vx**2 + vy**2 + vz**2)
self.translationalEnergy = 0.5 * totalMass * (vx ** 2 + vy ** 2 + vz ** 2)
self.translationalEnergy.setOutputs("Translational Kinetic Energy (J)")
self.kineticEnergy = self.rotationalEnergy + self.translationalEnergy
self.kineticEnergy.setOutputs("Kinetic Energy (J)")
Expand Down Expand Up @@ -2030,9 +2030,9 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
self.streamVelocityZ = -1 * self.vz
self.streamVelocityZ.setOutputs("Freestream Velocity Z (m/s)")
self.freestreamSpeed = (
self.streamVelocityX**2
+ self.streamVelocityY**2
+ self.streamVelocityZ**2
self.streamVelocityX ** 2
+ self.streamVelocityY ** 2
+ self.streamVelocityZ ** 2
) ** 0.5
self.freestreamSpeed.setOutputs("Freestream Speed (m/s)")
# Apogee Freestream speed
Expand All @@ -2052,15 +2052,15 @@ def postProcess(self, interpolation="spline", extrapolation="natural"):
self.maxReynoldsNumberTime = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 0]
self.maxReynoldsNumber = self.ReynoldsNumber[maxReynoldsNumberTimeIndex, 1]
# Dynamic Pressure
self.dynamicPressure = 0.5 * self.density * self.freestreamSpeed**2
self.dynamicPressure = 0.5 * self.density * self.freestreamSpeed ** 2
self.dynamicPressure.setOutputs("Dynamic Pressure (Pa)")
maxDynamicPressureTimeIndex = np.argmax(self.dynamicPressure[:, 1])
self.maxDynamicPressureTime = self.dynamicPressure[
maxDynamicPressureTimeIndex, 0
]
self.maxDynamicPressure = self.dynamicPressure[maxDynamicPressureTimeIndex, 1]
# Total Pressure
self.totalPressure = self.pressure * (1 + 0.2 * self.MachNumber**2) ** (3.5)
self.totalPressure = self.pressure * (1 + 0.2 * self.MachNumber ** 2) ** (3.5)
self.totalPressure.setOutputs("Total Pressure (Pa)")
maxtotalPressureTimeIndex = np.argmax(self.totalPressure[:, 1])
self.maxtotalPressureTime = self.totalPressure[maxtotalPressureTimeIndex, 0]
Expand Down Expand Up @@ -2507,7 +2507,7 @@ def calculateStallWindVelocity(self, stallAngle):
wV = (
2 * vF * math.cos(theta) / c
+ (
4 * vF * vF * math.cos(theta) * math.cos(theta) / (c**2)
4 * vF * vF * math.cos(theta) * math.cos(theta) / (c ** 2)
+ 4 * 1 * vF * vF / c
)
** 0.5
Expand Down Expand Up @@ -3166,7 +3166,7 @@ def plotFluidMechanicsData(self):

return None

def calculateFinFlutterAnalysis(self, finThickness, shearModulus):
def calculateFinFlutterAnalysis(self, finThickness, shearModulus, altMaxV):
"""Calculate, create and plot the Fin Flutter velocity, based on the
pressure profile provided by Atmospheric model selected. It considers the
Flutter Boundary Equation that is based on a calculation published in
Expand All @@ -3182,6 +3182,8 @@ def calculateFinFlutterAnalysis(self, finThickness, shearModulus):
The fin thickness, in meters
shearModulus : float
Shear Modulus of fins' material, must be given in Pascal
altMaxV : float
Altitude at maximum velocity of flight, in feet

Return
------
Expand All @@ -3194,11 +3196,14 @@ def calculateFinFlutterAnalysis(self, finThickness, shearModulus):
s = (self.rocket.tipChord + self.rocket.rootChord) * self.rocket.span / 2
ar = self.rocket.span * self.rocket.span / s
la = self.rocket.tipChord / self.rocket.rootChord

T = 59 - (0.0356 * altMaxV)
self.pressure = (self.pressure * 0) + (
2116 * ((T + 459.7) / 518.6) ** 5.256
) * 6894.7572931783
# Calculate the Fin Flutter Mach Number
self.flutterMachNumber = (
(shearModulus * 2 * (ar + 2) * (finThickness / self.rocket.rootChord) ** 3)
/ (1.337 * (ar**3) * (la + 1) * self.pressure)
/ (1.337 * (ar ** 3) * (la + 1) * self.pressure)
) ** 0.5

# Calculate difference between Fin Flutter Mach Number and the Rocket Speed
Expand Down
Loading