Merge pull request #133 from GodotMisogi/develop

Version 0.4.6: New fuselage geometry, stability tutorial, and bugfixes

Project.toml

@@ -1,7 +1,7 @@
 name = "AeroFuse"
 uuid = "477c59f4-51f5-487f-bf1e-8db39645b227"
 authors = ["GodotMisogi <[email protected]>"]
-version = "0.4.5"
+version = "0.4.6"
 
 [deps]
 Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"

README.md

@@ -5,7 +5,7 @@
 [![CI](https://github.com/GodotMisogi/AeroFuse.jl/workflows/CI/badge.svg?branch=main)](https://github.com/GodotMisogi/AeroFuse.jl/actions?query=workflow%3ACI+branch%3Amain)
 
 <p align="center">
-<img width="60%", src="https://raw.githubusercontent.com/GodotMisogi/AeroFuse.jl/main/docs/src/assets/logo.svg">
+<img width="60%", src="docs/src/assets/logo.svg">
 </p>
 
 AeroFuse is meant to be a toolbox for aircraft design analyses. It currently provides convenient methods for developing studies in aircraft geometry, aerodynamics and structures, with implementations in other relevant fields such as flight dynamics and propulsion in progress.
@@ -31,8 +31,16 @@ The current focus is to enable tutorials in computation in an aerospace educatio
 - Automatic differentiation support primarily via [`ForwardDiff.jl`](https://github.com/JuliaDiff/ForwardDiff.jl).
 - Optimization capabilities with these features.
 
+**Aerodynamic Analysis:**
+
 ![](plots/VortexLattice.svg)
-<!-- ![](plots/LinearVortex.svg) -->
+
+**Stability Analysis:**
+
+![](plots/boeing777.svg)
+
+**Optimization:**
+
 ![](plots/SciMLWingOptimization.svg)
 
 **In progress:**
@@ -61,7 +69,7 @@ If you use AeroFuse in your research, please cite the following until any releva
   author  = {Arjit Seth, Rhea P. Liem},
   title   = {AeroFuse},
   url     = {https://github.com/GodotMisogi/AeroFuse},
-  version = {0.4.5},
-  date    = {2023-03-20},
+  version = {0.4.6},
+  date    = {2023-03-27},
 }
 ```

docs/in_progress/flight-dynamics.jl

@@ -0,0 +1,38 @@
+
+# ## Flight Dynamics
+
+# ### Longitudinal Motion
+
+# A standard three degrees-of-freedom rigid body model is used for performing flight dynamics analyses in the 2-dimensional longitudinal plane formulated as an initial-value problem. The coupled differential equations are shown in a canonical state-space representation subject to initial conditions $\mathbf x_0$, with the time-evolution of the state vector $\mathbf x$ driven by the forcing function $\mathbf f$.
+
+# ```math
+# \begin{aligned}
+#     \frac{d\mathbf x}{dt} & = \mathbf f(\mathbf x, t) \\
+#     \frac{d}{dt}
+#     \begin{bmatrix}
+#         u_b \\
+#         w_b \\
+#         Q \\
+#         x_e \\
+#         y_e \\
+#         \Theta \\
+#         \delta_e \\
+#         m
+#     \end{bmatrix} & =
+#     \begin{bmatrix}
+#         -Qw_b + (T - D(V_\infty)\cos\alpha + L(V_\infty, \delta_e)\sin\alpha - W\sin\Theta) / m \\
+#         Qu_b + (  - D(V_\infty)\sin\alpha - L(V_\infty, \delta_e)\cos\alpha + W\cos\Theta) / m \\
+#         (M_A(\alpha, \delta_e, \hat Q) - T\Delta_{zT}) / I_{yy} \\
+#         u_b \cos\Theta - w_b \sin\Theta \\
+#         u_b \sin\Theta + w_b \cos\Theta \\
+#         Q \\ 
+#         g(\mathbf x, t) \\
+#         -c_T T
+#     \end{bmatrix}
+# \end{aligned}
+# ```
+# The fuel burn over time (viz. reduction of mass $m$) is computed using a specific fuel consumption value with a linear dependence on the thrust. A manual controller law for the elevator deflection angle $\delta_e$ can also be implemented by providing the function $g(\mathbf x, t)$.
+
+# ### Full Space
+
+# A standard six degrees-of-freedom rigid body model is used for performing flight dynamics analyses in 3 dimensions.

docs/in_progress/tutorials-structures.jl

@@ -0,0 +1,45 @@
+
+# # ## Euler-Bernoulli Beam Structural Analysis
+# # The tubular beam's relevant properties, viz. the Young's (elastic) modulus $E$, shear modulus $G$, and torsional moment of inertia $J$ must be specified to define the stiffness matrix for its discretization with $n$ elements.
+
+# ## Material properties
+# E = 1. # Elastic modulus
+# G = 1. # Shear modulus
+# J = 2. # Torsional moment of inertia
+# n = 2  # Number of sections
+
+# ## Stiffness matrix
+# K = bending_stiffness_matrix(
+#     fill(E, 2), 
+#     fill(G, 2),
+#     fill(J, 2), 
+#     :z         # Direction of deflection
+# )
+
+# # Fixed, hinged beam subjected to force and moment at the center.
+
+# ## Stiffness matrix
+# A = K[[3,4,6],[3,4,6]]  # v2, φ2, φ3
+
+# ## Load vector
+# b = [-1000, 1000, 0]    # F2, M2, M3
+
+# ## Solution
+# x = A \ b
+
+# ## Forces
+# F1 = K * [ 0.; 0.; x[1:2]; 0.; x[3] ]
+
+# # Propped cantilever beam subjected to force at free end.
+
+# ## Stiffness matrix
+# A = K[[1,2,4],[1,2,4]] # v1, φ1, φ2
+
+# ## Load vector
+# b = [10, 0, 0]
+
+# ## Solution
+# x = A \ b
+
+# ## Forces
+# F2 = K * [ x[1:2]; 0.; x[3]; 0.; 0. ]

docs/lit/howto.jl

@@ -26,11 +26,9 @@ plot(xlabel = L"x", ylabel = L"y", aspect_ratio = 1)
 plot!(my_foil)
 
 # !!! tip
-#     You can also import airfoil coordinates from the internet by using `FileIO`.
+#     You can also import airfoil coordinates from the internet.
 
 ## Read coordinates from an internet address.
-using FileIO
-
 online_foil = read_foil(download("https://m-selig.ae.illinois.edu/ads/coord/s1210.dat"))
 
 plot!(online_foil)
@@ -93,7 +91,11 @@ locs     = panel_location.(panels);   # Upper or lower surface
 # ## Wing Geometry
 # How to work with wing geometry.
 # 
-# To define a wing, AeroFuse provides a `Wing` constructor.
+# To define a wing, AeroFuse provides a `Wing` constructor based on the following parametrization. The named arguments correspond to the foil shapes, chord and span lengths, twist, dihedral and sweep angles.
+# !!! info
+#    A **wing section** consists of two foil profiles and their chord lengths and twist angles. Between them is their span length with associated _leading-edge_ dihedral and sweep angles. So a general half-wing consisting of ``n`` sections will have ``n`` entries for spans $b$, dihedrals $\delta$, sweeps $\Lambda$, and ``n+1`` entries for foils, chords $c$, and twists $\iota$, for some ``n \in \mathbb N``.
+# 
+# ![](https://godot-bloggy.xyz/post/diagrams/WingGeometry.svg)
 airfoil = naca4((2,4,1,2))
 wing = Wing(
     foils     = [ airfoil for i in 1:3 ],   # Foils
@@ -105,10 +107,10 @@ wing = Wing(
     w_sweep   = 0.25,             # Sweep angle location w.r.t. 
                                   ## normalized chord lengths ∈ [0,1]
     symmetry  = true,             # Whether wing is symmetric
-    ## flip      = false           # Whether wing is reflected
+    ## flip      = false           # Whether wing is flipped in x-z plane
 )
 
-# The `symmetry` Boolean argument determines whether the geometry should be reflected in the ``x``-``z`` plane.
+# The `symmetry` Boolean argument specifies whether the geometry should be reflected in the ``x``-``z`` plane. The `flip` Boolean argument specifies whether the coordinates should be flipped in the ``x``-``z`` plane. The `w_sweep` argument specifies the chordwise-ratio of the sweep angles, e.g. 0. = leading edge sweep angle (default), 1. = trailing edge, 0.25 = quarter-chord.
 
 # The following "getter" functions provide quantities of interest such as chord lengths, spans, twist, dihedral, and sweep angles.
 
@@ -118,22 +120,43 @@ wing = Wing(
 # spans(wing)
 # twists(wing)
 # dihedrals(wing)
-# sweeps(wing)
+# sweeps(wing) # Leading-edge sweep angles
+# sweeps(wing, 0.25) # Quarter-chord sweep angles
+# sweeps(wing, 1.0) # Trailing-edge sweep angles
 # ```
 
-# You can evaluate commonly defined properties such as the aspect ratio, projected area, etc. with the following functions.
+# You can evaluate commonly defined properties such as the aspect ratio, projected area, etc., with the following functions.
 
 # ```@repl howto
 # aspect_ratio(wing)
 # projected_area(wing)
 # span(wing)
 # taper_ratio(wing)
+# mean_aerodynamic_chord(wing)
 # ```
 
+# The mean aerodynamic center calculation merits discussion. Its calculation procedure is depicted in the following diagram:
+# 
+# ![](https://raw.githubusercontent.com/HKUST-OCTAD-LAB/MECH3620Materials/main/pics/WingParams.svg)
+mac25_w = mean_aerodynamic_center(wing) # 25% MAC by default
+mac40_w = mean_aerodynamic_center(wing, 0.4) # at 40% MAC
+
 # You can plot the wing with `Plots.jl` quite simply.
-plot(wing, zlim = (-0.5, 0.5) .* span(wing), aspect_ratio = 1, label = "Wing")
+plot(wing, 
+    0.25; # Aerodynamic center location in terms of chord ratio plot 
+    mac = true, # Optional named argument to disable aerodynamic center plot
+    zlim = (-0.5, 0.5) .* span(wing), 
+    aspect_ratio = 1, 
+    label = "Wing",
+)
 
-# There is also a convenient function for pretty-printing information, in which the first argument takes the `Wing` type and the second takes a name (as a `String` or `Symbol`).
+# At 40% MAC
+plot!(wing, 
+    0.40; # 40% MAC
+    label = "40% MAC"
+)
+
+# There is also a convenient function for pretty-printing information (using [PrettyTables.jl](https://ronisbr.github.io/PrettyTables.jl/stable/)), in which the first argument takes the `Wing` type and the second takes a name (as a `String` or `Symbol`).
 print_info(wing, "Wing")
 
 #md # !!! tip
@@ -143,17 +166,15 @@ print_info(wing, "Wing")
 using Accessors
 
 ## Set only chords with other properties remaining identical.
-wing = @set wing.chords = [0.4, 0.1, 0.05]
-
-print_info(wing, "My Wing")
+new_wing = @set wing.chords = [0.4, 0.1, 0.05]
 
 # ## Vortex Lattice Aerodynamic Analyses
 # 
 # How to run a generic aerodynamic analysis on a conventional aircraft configuration.
 #
 # ### Geometry 
 # First we define the lifting surfaces. These can be a combination of `Wing` types constructed using the various methods available.
-# !!! warning "Alert"
+# !!! info
 #     Support for fuselages and control surfaces is in progress.
 
 # Horizontal tail
@@ -196,7 +217,8 @@ vtail = WingSection(
 ## Wing meshes
 wing_mesh = WingMesh(wing, [20,8], 10)
 
-# Optionally the type of spanwise spacing can be specified by the keyword `span_spacing` and providing types `Sine(), Cosine(), Uniform()` or a vector of the combination with length corresponding to the number of sections.
+# !!! tip 
+#     Optionally ,the type of spanwise spacing can be specified by the keyword `span_spacing` and providing types `Sine(), Cosine(), Uniform()` or a vector of the combination with length corresponding to the number of sections. For example, if you have two spanwise sections and the wing is symmetric, then the total number of spanwise sections is four. So the spanwise spacing argument would be, for example, `spanwise_spacing = [Uniform(), Cosine(), Cosine(), Uniform()]`
 
 ## Horizontal tail
 htail_mesh = WingMesh(htail, [10], 5;
@@ -211,7 +233,7 @@ vtail_mesh = WingMesh(vtail, [10], 5); # Vertical tail
 wing_horsies = make_horseshoes(wing_mesh)
 
 # !!! note
-#     You can also generate vortex ring elements by calling `make_vortex rings`.
+#     You can also generate vortex ring elements by calling `make_vortex rings`, but the aerodynamic analysis is slightly slower.
 
 # For multiple lifting surfaces, it is most convenient to define a single vector consisting of all the components' horseshoes using [ComponentArrays.jl](https://github.com/jonniedie/ComponentArrays.jl).
 aircraft = ComponentVector(
@@ -258,9 +280,12 @@ system.influence_matrix[:wing, :htail];
 
 system.boundary_vector[:wing]
 
-# ### Dynamics
+# ### Forces and Moments
 # 
-# The analysis computes the aerodynamic forces by surface pressure integration for nearfield forces. You can specify the axis system for the nearfield forces, with choices of geometry, body, wind, or stability axes. The wind axes are used by default.
+# The analysis computes the aerodynamic force coefficients $C_{\mathbf{F}}$ by surface pressure integration for nearfield forces. The aerodynamic moment coefficients $C_{\mathbf M}$ are calculated by computing the cross product of the moment arm of each control point from the reference point $\mathbf r_i - \mathbf r_{ref}$, i.e. $(C_{\mathbf{M}})_i = (\mathbf r_i - \mathbf r_{ref}) \times (C_{\mathbf{F}})_i$
+
+# !!! note
+#     You can specify the axis system for the nearfield forces, with choices of geometry, body, wind, or stability axes. The wind axes are used by default.
 
 ## Compute dynamics
 ax = Wind() # Body(), Stability(), Geometry()
@@ -304,47 +329,13 @@ dv_data = freestream_derivatives(
     print_components = true,  # Prints the results for all components
 );
 
-# ## Euler-Bernoulli Beam Structural Analysis
-# The tubular beam's relevant properties, viz. the Young's (elastic) modulus $E$, shear modulus $G$, and torsional moment of inertia $J$ must be specified to define the stiffness matrix for its discretization with $n$ elements.
-
-## Material properties
-E = 1. # Elastic modulus
-G = 1. # Shear modulus
-J = 2. # Torsional moment of inertia
-n = 2  # Number of sections
-
-## Stiffness matrix
-K = bending_stiffness_matrix(
-    fill(E, 2), 
-    fill(G, 2),
-    fill(J, 2), 
-    :z         # Direction of deflection
-)
-
-# Fixed, hinged beam subjected to force and moment at the center.
-
-## Stiffness matrix
-A = K[[3,4,6],[3,4,6]]  # v2, φ2, φ3
-
-## Load vector
-b = [-1000, 1000, 0]    # F2, M2, M3
-
-## Solution
-x = A \ b
-
-## Forces
-F1 = K * [ 0.; 0.; x[1:2]; 0.; x[3] ]
-
-# Propped cantilever beam subjected to force at free end.
-
-## Stiffness matrix
-A = K[[1,2,4],[1,2,4]] # v1, φ1, φ2
-
-## Load vector
-b = [10, 0, 0]
-
-## Solution
-x = A \ b
-
-## Forces
-F2 = K * [ x[1:2]; 0.; x[3]; 0.; 0. ]
+#md # !!! note
+#md #     For efficiency, instead of calling `solve_case` to compute the forces and then computing the derivatives, you can directly call:
+#md #     ```julia
+#md #     dvs = freestream_derivatives(
+#md #               aircraft, fs, refs,
+#md #               compressible     = true, # Compressibility option
+#md #               print            = true, # Prints the results for only the aircraft
+#md #               print_components = true, # Prints the results for all components
+#md #           )
+#md #     ```

docs/lit/theory.jl

@@ -137,48 +137,6 @@
 
 # ```
 
-# ## Flight Dynamics
-
-# ### Longitudinal Motion
-
-# A standard three degrees-of-freedom rigid body model is used for performing flight dynamics analyses in the 2-dimensional longitudinal plane formulated as an initial-value problem. The coupled differential equations are shown in a canonical state-space representation subject to initial conditions $\mathbf x_0$, with the time-evolution of the state vector $\mathbf x$ driven by the forcing function $\mathbf f$.
-
-# ```math
-# \begin{aligned}
-#     \frac{d\mathbf x}{dt} & = \mathbf f(\mathbf x, t) \\
-#     \frac{d}{dt}
-#     \begin{bmatrix}
-#         u_b \\
-#         w_b \\
-#         Q \\
-#         x_e \\
-#         y_e \\
-#         \Theta \\
-#         \delta_e \\
-#         m
-#     \end{bmatrix} & =
-#     \begin{bmatrix}
-#         -Qw_b + (T - D(V_\infty)\cos\alpha + L(V_\infty, \delta_e)\sin\alpha - W\sin\Theta) / m \\
-#         Qu_b + (  - D(V_\infty)\sin\alpha - L(V_\infty, \delta_e)\cos\alpha + W\cos\Theta) / m \\
-#         (M_A(\alpha, \delta_e, \hat Q) - T\Delta_{zT}) / I_{yy} \\
-#         u_b \cos\Theta - w_b \sin\Theta \\
-#         u_b \sin\Theta + w_b \cos\Theta \\
-#         Q \\ 
-#         g(\mathbf x, t) \\
-#         -c_T T
-#     \end{bmatrix}
-# \end{aligned}
-# ```
-# The fuel burn over time (viz. reduction of mass $m$) is computed using a specific fuel consumption value with a linear dependence on the thrust. A manual controller law for the elevator deflection angle $\delta_e$ can also be implemented by providing the function $g(\mathbf x, t)$.
-
-# ### Full Space
-
-# A standard six degrees-of-freedom rigid body model is used for performing flight dynamics analyses in 3 dimensions.
-
-# ```math
-
-# ```
-
 # ## References
 
 # 1. Mark Drela. _Flight Vehicle Aerodynamics_. MIT Press, 2014.