Merge pull request #110 from GodotMisogi/develop

Added docs and more bugfixes in profile drag

Author: Arjit Seth

Project.toml

@@ -1,7 +1,7 @@
 name = "AeroMDAO"
 uuid = "ae1513eb-aba2-415b-9b88-80f66c7c7c76"
 authors = ["GodotMisogi <[email protected]>"]
-version = "0.3.8"
+version = "0.3.9"
 
 [deps]
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"

README.md

@@ -60,7 +60,7 @@ If you use AeroMDAO in your research, please cite the following until any releva
   author  = {Arjit Seth, Stephane Redonnet, Rhea P. Liem},
   title   = {AeroMDAO},
   url     = {https://github.com/GodotMisogi/AeroMDAO},
-  version = {0.3.8},
+  version = {0.3.9},
   date    = {2022-4-06},
 }
 ```

docs/src/index.md

@@ -52,7 +52,7 @@ If you use AeroMDAO in your research, please cite the following until any releva
   author  = {Arjit Seth, Rhea P. Liem, Stephane Redonnet},
   title   = {AeroMDAO},
   url     = {https://github.com/GodotMisogi/AeroMDAO},
-  version = {0.3.8},
+  version = {0.3.9},
   date    = {2021-04-06},
 }
 ```

examples/vortex_lattice_method/vlm_aircraft.jl

@@ -93,36 +93,36 @@ end;
 ## Viscous drag prediction
 
 # Equivalent flat-plate skin friction estimation
-CDv_wing  = profile_drag_coefficient(wing_mesh, [0.8, 0.8], system.reference)
-CDv_htail = profile_drag_coefficient(htail_mesh, [0.6, 0.6], system.reference)
-CDv_vtail = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], system.reference)
+CDv_wing  = profile_drag_coefficient(wing,  [0.8, 0.8], system.reference)
+CDv_htail = profile_drag_coefficient(htail, [0.6, 0.6], system.reference)
+CDv_vtail = profile_drag_coefficient(vtail, [0.6, 0.6], system.reference)
 
 CDv_plate = CDv_wing + CDv_htail + CDv_vtail
 
-## Local dissipation form factor friction estimation
+## Local dissipation form factor friction estimation (WRONG???)
 import LinearAlgebra: norm
 
 edge_speeds = norm.(surface_velocities(system)); # Inviscid speeds on the surfaces
 
 # Drag coefficients
-CDvd_wing   = profile_drag_coefficient(wing_mesh,  [0.8, 0.8], edge_speeds.wing,  refs)
-CDvd_htail  = profile_drag_coefficient(htail_mesh, [0.6, 0.6], edge_speeds.htail, refs)
-CDvd_vtail  = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], edge_speeds.vtail, refs)
+CDvd_wing  = profile_drag_coefficient(wing_mesh,  [0.8, 0.8], edge_speeds.wing,  ref)
+CDvd_htail = profile_drag_coefficient(htail_mesh, [0.6, 0.6], edge_speeds.htail, ref)
+CDvd_vtail = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], edge_speeds.vtail, ref)
 
-CDv_diss    = CDvd_wing + CDvd_htail + CDvd_vtail
+CDv_diss   = CDvd_wing + CDvd_htail + CDvd_vtail
 
 ## Viscous drag coefficient
 CDv = CDv_plate
 
-## Total force coefficients with viscous drag prediction
+## Total force coefficients with empirical viscous drag prediction
 CDi_nf, CY_nf, CL_nf, Cl, Cm, Cn = nf = nearfield(system) 
 CDi_ff, CY_ff, CL_ff = ff = farfield(system)
 
 nf_v = (CD = CDi_nf + CDv, CDv = CDv, nf...)
 ff_v = (CD = CDi_ff + CDv, CDv = CDv, ff...)
 
 ## Spanwise forces/lifting line loads
-wing_ll  = spanwise_loading(chord_panels(wing_mesh), CFs.wing, S)
+wing_ll  = spanwise_loading(chord_panels(wing_mesh),  CFs.wing,  S)
 htail_ll = spanwise_loading(chord_panels(htail_mesh), CFs.htail, S)
 vtail_ll = spanwise_loading(chord_panels(vtail_mesh), CFs.vtail, S);
 

examples/vortex_lattice_method/vlm_wing.jl

@@ -60,11 +60,10 @@ FFs = farfield_coefficients(system)
 comp_coeffs = mapreduce(name -> [ sum(CFs[name]); sum(CMs[name]); FFs[name] ], hcat, keys(system.vortices))
 
 ## Equivalent flat-plate skin-friction estimation
-x_tr        = [0.98, 0.98]              # Transition locations over sections
-CDv_plate   = profile_drag_coefficient(wing, x_tr, refs)
+x_tr        = fill(0.98, 4)              # Transition locations over sections
+CDv_plate   = profile_drag_coefficient(wing_mesh, x_tr, refs)
 
 ## Local-dissipation drag estimation (WRONG???)
-x_tr        = fill(0.98, 4)
 cam_panels  = camber_panels(wing_mesh)
 edge_speeds = surface_velocities(system).wing
 CDv_diss    = profile_drag_coefficient(wing, x_tr, edge_speeds, cam_panels, refs)

src/AeroMDAO.jl

@@ -95,9 +95,9 @@ import .AircraftGeometry: Fuselage, projected_area, length, cosine_interpolation
 export Fuselage, projected_area, length, cosine_interpolation
 
 # Wing
-import .AircraftGeometry: HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties
+import .AircraftGeometry: HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties, wetted_area_ratio
 
-export HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties
+export HalfWing, HalfWingSection, Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, aspect_ratio, projected_area, taper_ratio, leading_edge, trailing_edge, chop_leading_edge, chop_trailing_edge, chop_wing, chop_sections, chop_coordinates, chop_spanwise_sections, chop_chords, chop_spans, wing_bounds, make_panels, mesh_chords, mesh_wing, mesh_cambers, max_thickness_to_chord_ratio_sweeps, mean_aerodynamic_center, panel_wing, number_of_spanwise_panels, symmetric_spacing, coordinates, chord_coordinates, camber_coordinates, surface_coordinates, foils, chords, twists, spans, dihedrals, sweeps, position, orientation, WingMesh, chord_panels, camber_panels, surface_panels, AbstractSpacing, Sine, Cosine, Uniform, properties, wetted_area_ratio
 
 # Surfaces
 import .AircraftGeometry: HorizontalTail, VerticalTail

src/Aerodynamics/profile_drag.jl

@@ -10,7 +10,7 @@ function form_factor(wing :: HalfWing, M)
     Kf  = form_factor_wing.(xcs, tcs, sweeps, M)
 end
 
-form_factor(wing :: Wing, M) = (form_factor(wing.left, M) + form_factor(wing.right, M)) / 2
+form_factor(wing :: Wing, M) = [ form_factor(wing.left, M); form_factor(wing.right, M) ]
 
 # Schlichting averaged skin-friction coefficients
 cf_lam(Re_c, k_lam = 1.) = 1.328 / √(Re_c * k_lam)
@@ -33,14 +33,16 @@ function wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
 
     # Profile drag
     Dp_by_q = sum(@. cfs * S_wets * K_fs)
+
+    return Dp_by_q
 end
 
-function wetted_area_drag(wing :: HalfWing, x_tr, V, ρ, a_ref = 330., μ = 1.5e-5)
+function wetted_area_drag(wing :: AbstractWing, x_tr, V, ρ, a_ref, μ)
     # Chord processing
     mean_chords = (forward_sum ∘ chords)(wing) / 2
 
     # Wetted areas
-    S_wets  = @. mean_chords * wing.spans / cos(wing.dihedrals)
+    S_wets  = mean_chords .* spans(wing) ./ cos.(dihedrals(wing))
 
     # Form factors
     M       = V / a_ref
@@ -49,20 +51,7 @@ function wetted_area_drag(wing :: HalfWing, x_tr, V, ρ, a_ref = 330., μ = 1.5e
     wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
 end
 
-function wetted_area_drag(wing :: WingMesh, x_tr, V, ρ, a_ref = 330., μ = 1.5e-5)
-    # Chord processing
-    mean_chords = (forward_sum ∘ chords)(wing.surface) / 2
-
-    # Wetted areas
-    surf_pans = camber_panels(wing)
-    S_wets    = sum(panel_area, surf_pans, dims = 1)
-
-    # Form factors
-    M       = V / a_ref
-    K_fs    = form_factor(wing.surface, M)
-
-    wetted_area_drag(mean_chords, S_wets, K_fs, x_tr, V, ρ, M, μ)
-end
+wetted_area_drag(wing :: WingMesh, x_tr, V, ρ, a_ref, μ) = wetted_area_drag(wing.surface, x_tr, V, ρ, a_ref, μ)
 
 # Sato's local-friction and local-dissipation based on power balance method from Mark Drela, Flight Vehicle Aerodynamics, eq. 4.115.
 function local_dissipation_drag(wing :: Wing, wetted_areas, ρ_es, u_es, x_tr, V, ρ, M, μ)
@@ -79,8 +68,10 @@ function local_dissipation_drag(wing :: WingMesh, ρ_es, u_es, x_tr, V, ρ, M, 
     # Chord processing
     mean_chords = (forward_sum ∘ chords)(wing.surface) / 2
 
-    # Compute weighted wetted areas based on inviscid edge velocity distribution.
+    # Compute wetted area of camber (for now).
     S_wets = panel_area.(camber_panels(wing))
+
+    # Calculate weighted wetted areas based on inviscid edge velocity distribution.
     weighted_S_wets = sum(@. ρ_es * u_es^3 * S_wets; dims = 1) ./ (ρ * V^3)
 
     wetted_area_drag(mean_chords, weighted_S_wets, 1., x_tr, V, ρ, M, μ)
@@ -91,18 +82,26 @@ abstract type AbstractProfileDrag end
 struct FormFactor <: AbstractProfileDrag end
 struct Dissipation <: AbstractProfileDrag end
 
-profile_drag_coefficient(wing :: HalfWing, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
-profile_drag_coefficient(wing :: WingMesh, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
-profile_drag_coefficient(wing :: Wing, x_tr, V, rho_ref, a_ref, area_ref, μ) = profile_drag_coefficient(wing.left, x_tr, V, rho_ref, a_ref, area_ref, μ) + profile_drag_coefficient(wing.right, x_tr, V, rho_ref, a_ref, area_ref, μ)
+profile_drag_coefficient(wing :: AbstractWing, x_tr, V, rho_ref, a_ref, area_ref, μ) = wetted_area_drag(wing, x_tr, V, rho_ref, a_ref, μ) / area_ref
+
+"""
+    profile_drag_coefficient(wing :: AbstractWing, x_tr, refs :: References)
 
+Estimate the profile drag coefficient of a `Wing`` using the **wetted-area method** based on Schlichting's skin-friction coefficient formula.
+"""
 profile_drag_coefficient(wing :: AbstractWing, x_tr, refs :: References) = profile_drag_coefficient(wing, x_tr, refs.speed, refs.density, refs.sound_speed, refs.area, refs.viscosity)
 
 profile_drag_coefficient(wing :: AbstractWing, x_tr, edge_speeds, panels, refs :: References) = local_dissipation_drag(wing, panel_area.(panels), refs.density, edge_speeds, x_tr, refs.speed, refs.density, mach_number(refs), refs.viscosity) / refs.area
 
+"""
+    profile_drag_coefficient(wing :: WingMesh, x_tr, edge_speeds, panels, refs :: References)
+
+Estimate the profile drag coefficient of a `Wing` using the **local-friction and local-dissipation method** based on Schlichting's skin-friction coefficient formula.
+"""
 profile_drag_coefficient(wing :: WingMesh, x_tr, edge_speeds, refs :: References) = local_dissipation_drag(wing, refs.density, edge_speeds, x_tr, refs.speed, refs.density, mach_number(refs), refs.viscosity) / refs.area
 
 function wave_drag(M, Λ, t_by_c, Cl, κ_A)
     M_drag_divergence = κ_A / cos(Λ) - t_by_c / cos(Λ)^2 - Cl / (10 * cos(Λ)^3) # Drag divergence Mach number
-    M_crit   = M_drag_divergnece - (0.1 / 80)^(1/3) # Critical Mach number
+    M_crit  = M_drag_divergence - (0.1 / 80)^(1/3) # Critical Mach number
     CD_wave = ifelse(M > M_crit, 20 * (Mach - M_crit)^4, 0.)
 end

src/Geometry/AircraftGeometry/AircraftGeometry.jl

@@ -17,7 +17,7 @@ using Interpolations
 import ..MathTools: uniform_spacing, linear_spacing, sine_spacing, cosine_spacing, cosine_interp, splitat, adj3, slope, columns, forward_sum, forward_division, forward_difference, weighted_vector, vectarray, extend_yz
 
 # Panel geometry
-import ..PanelGeometry: Panel2D, Panel3D, panel_area, normal_vector, transform, make_panels
+import ..PanelGeometry: Panel2D, Panel3D, panel_area, normal_vector, transform, make_panels, wetted_area
 
 import ..AeroMDAO: properties
 

src/Geometry/AircraftGeometry/Wings/mesh_wing.jl

@@ -86,6 +86,8 @@ chord_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spaci
 
 camber_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = camber_coordinates(wing, number_of_spanwise_panels(wing, span_num), chord_num; spacings = spacings)
 
+surface_coordinates(wing :: Wing, span_num :: Integer, chord_num :: Integer; spacings = symmetric_spacing(wing)) = surface_coordinates(wing, number_of_spanwise_panels(wing, span_num), chord_num; spacings = spacings)
+
 ## Panelling
 #==========================================================================================#
 
@@ -164,11 +166,37 @@ check_definition(surf :: HalfWing, n_span) = @assert length(n_span) == length(su
 check_definition(surf :: Wing, n_span) = @assert length(n_span) == length(surf.right.spans) == length(surf.left.spans) "The spanwise number vector's length must be the same as the number of sections of the surface."
 
 ##
+"""
+    chord_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)
+
+Generate the chord coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
+"""
 chord_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = chord_coordinates(wing.surface, n_span, n_chord)
+
+"""
+    camber_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)
+
+Generate the camber coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
+"""
 camber_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = camber_coordinates(wing.surface, n_span, n_chord)
+
+"""
+    surface_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)
+
+Generate the surface coordinates of a `WingMesh` with default spanwise ``n_s`` and chordwise ``n_c`` panel distributions from the mesh.
+"""
 surface_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) = surface_coordinates(wing.surface, n_span, n_chord)
 
-surface_panels(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord)  = (make_panels ∘ surface_coordinates)(wing, n_span, n_chord)
+"""
+    surface_panels(wing_mesh :: WingMesh, 
+                   n_s = wing_mesh.num_span, 
+                   n_c = length(first(foils(wing_mesh.surface))).x)
+
+Generate the surface panel distribution from a `WingMesh` with the default spanwise ``n_s`` panel distribution from the mesh and the chordwise panel ``n_c`` distribution from the airfoil.
+
+In case of strange results, provide a higher number of chordwise panels to represent the airfoils more accurately
+""" 
+surface_panels(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x))  = (make_panels ∘ surface_coordinates)(wing, n_span, n_chord)
 
 """
     chord_panels(wing_mesh :: WingMesh)
@@ -184,6 +212,24 @@ Generate the camber panel distribution from a `WingMesh`.
 """
 camber_panels(wing :: WingMesh) = make_panels(wing.camber_mesh)
 
+"""
+    wetted_area(wing_mesh :: WingMesh, 
+                n_s = wing_mesh.num_span, 
+                n_c = length(first(foils(wing_mesh.surface))).x)
+
+Determine the wetted area ``S_{wet}`` of a `WingMesh` by calculating the total area of the surface panels.
+"""
+wetted_area(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x)) = wetted_area(surface_panels(wing, n_span, n_chord))
+
+"""
+    wetted_area_ratio(wing_mesh :: WingMesh, 
+                      n_s = wing_mesh.num_span, 
+                      n_c = length(first(foils(wing_mesh.surface))).x)
+
+Determine the wetted area ratio ``S_{wet}/S`` of a `WingMesh` by calculating the ratio of the total area of the surface panels to the projected area of the `Wing`.
+"""
+wetted_area_ratio(wing :: WingMesh, n_span = wing.num_span, n_chord = length(first(foils(wing.surface)).x)) = wetted_area(wing, n_span, n_chord) / projected_area(wing.surface)
+
 function Base.show(io :: IO, mesh :: WingMesh)
     n_c, n_s = size(mesh.chord_mesh) .- 1
     println(io, "WingMesh —")