more efficient polars

Author: 1-Bart-1

examples/ram_air_kite.jl

@@ -35,48 +35,48 @@ vel_app = [
 ] * v_a
 body_aero.va = vel_app
 
-# Plotting geometry
-plot_geometry(
-    body_aero,
-    "";
-    data_type=".svg",
-    save_path="",
-    is_save=false,
-    is_show=true,
-    view_elevation=15,
-    view_azimuth=-120
-)
+# # Plotting geometry
+# plot_geometry(
+#     body_aero,
+#     "";
+#     data_type=".svg",
+#     save_path="",
+#     is_save=false,
+#     is_show=true,
+#     view_elevation=15,
+#     view_azimuth=-120
+# )
 
 # Solving and plotting distributions
 @time results = solve(VSM, body_aero)
 @time results = solve(VSM, body_aero)
 
 CAD_y_coordinates = [panel.aerodynamic_center[2] for panel in body_aero.panels]
 
-plot_distribution(
-    [CAD_y_coordinates],
-    [results],
-    ["VSM"];
-    title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))",
-    data_type=".pdf",
-    is_save=false,
-    is_show=true
-)
+# plot_distribution(
+#     [CAD_y_coordinates],
+#     [results],
+#     ["VSM"];
+#     title="CAD_spanwise_distributions_alpha_$(round(aoa, digits=1))_beta_$(round(side_slip, digits=1))_yaw_$(round(yaw_rate, digits=1))_v_a_$(round(v_a, digits=1))",
+#     data_type=".pdf",
+#     is_save=false,
+#     is_show=true
+# )
 
-plot_polars(
-    [VSM],
-    [body_aero],
-    [
-        "VSM from Ram Air Kite OBJ and DAT file",
-    ];
-    angle_range=range(0, 20, length=20),
-    angle_type="angle_of_attack",
-    angle_of_attack=0,
-    side_slip=0,
-    v_a=10,
-    title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)",
-    data_type=".pdf",
-    is_save=false,
-    is_show=true
-)
+# plot_polars(
+#     [VSM],
+#     [body_aero],
+#     [
+#         "VSM from Ram Air Kite OBJ and DAT file",
+#     ];
+#     angle_range=range(0, 20, length=20),
+#     angle_type="angle_of_attack",
+#     angle_of_attack=0,
+#     side_slip=0,
+#     v_a=10,
+#     title="ram_kite_panels_$(wing.n_panels)_distribution_$(wing.spanwise_panel_distribution)",
+#     data_type=".pdf",
+#     is_save=false,
+#     is_show=true
+# )
 nothing

src/body_aerodynamics.jl

@@ -226,7 +226,7 @@ function calculate_panel_properties(section_list::Vector{Section}, n_panels::Int
 end
 
 """
-    calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String, 
+    calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String, 
                          core_radius_fraction::Float64,
                          va_norm_array::Vector{Float64}, 
                          va_unit_array::Matrix{Float64})
@@ -236,7 +236,7 @@ Calculate Aerodynamic Influence Coefficient matrices.
 Returns:
     Tuple of (AIC_x, AIC_y, AIC_z) matrices
 """
-function calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String,
+function calculate_AIC_matrices!(body_aero::BodyAerodynamics, model::String,
                               core_radius_fraction::Float64,
                               va_norm_array::Vector{Float64}, 
                               va_unit_array::Matrix{Float64})
@@ -249,9 +249,9 @@ function calculate_AIC_matrices(body_aero::BodyAerodynamics, model::String,
     velocity_induced, tempvel, va_unit, U_2D = zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3)
 
     # Calculate influence coefficients
-    for icp in 1:length(body_aero.panels)
+    for icp in eachindex(body_aero.panels)
         ep = getproperty(body_aero.panels[icp], evaluation_point)
-        for jring in 1:body_aero.n_panels
+        for jring in eachindex(body_aero.panels)
             va_unit .= @views va_unit_array[jring, :]
             filaments = body_aero.panels[jring].filaments
             va_norm = va_norm_array[jring]

src/panel.jl

@@ -36,7 +36,9 @@ mutable struct Panel
     cl_coefficients::Union{Nothing,Vector{Float64}}
     cd_coefficients::Union{Nothing,Vector{Float64}}
     cm_coefficients::Union{Nothing,Vector{Float64}}
-    polar_data::Union{Nothing, Tuple{Extrapolation}}
+    cl_interp::Function
+    cd_interp::Function
+    cm_interp::Function
     aerodynamic_center::MVec3
     control_point::MVec3
     bound_point_1::MVec3
@@ -78,10 +80,8 @@ mutable struct Panel
         aero_model = isa(section_1.aero_input, String) ? section_1.aero_input : section_1.aero_input[1]
 
         # Initialize aerodynamic properties
-        cl_coeffs = nothing
-        cd_coeffs = nothing
-        cm_coeffs = nothing
-        polar_data = nothing
+        cl_coeffs, cd_coeffs, cm_coeffs = nothing, nothing, nothing
+        cl_interp, cd_interp, cm_interp = ()->nothing, ()->nothing, ()->nothing
 
         if aero_model == "lei_airfoil_breukels"
             cl_coeffs, cd_coeffs, cm_coeffs = compute_lei_coefficients(section_1, section_2)
@@ -92,29 +92,21 @@ mutable struct Panel
                 throw(ArgumentError("Polar data must have same shape"))
             end
             polar_data = (aero_1 + aero_2) / 2
-            cl_interp = linear_interpolation(
-                polar_data[:,1],
-                polar_data[:,2];
-                extrapolation_bc=NaN
-            )
-            cd_interp = linear_interpolation(
-                polar_data[:,1],
-                polar_data[:,3];
-                extrapolation_bc=NaN
-            )
-            cm_interp = linear_interpolation(
-                polar_data[:,1],
-                polar_data[:,4];
-                extrapolation_bc=NaN
-            )
-            polar_data = (cl_interp, cd_interp, cm_interp)
+            cl_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,2]; extrapolation_bc=NaN)
+            cd_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,3]; extrapolation_bc=NaN)
+            cm_interp_ = linear_interpolation(polar_data[:,1], polar_data[:,4]; extrapolation_bc=NaN)
+            cl_interp = (α) -> cl_interp_(α)
+            cd_interp = (α) -> cd_interp_(α)
+            cm_interp = (α) -> cm_interp_(α)
         elseif aero_model == "interpolations"
             cl_left, cd_left, cm_left = section_1.aero_input[2]
             cl_right, cd_right, cm_right = section_2.aero_input[2]
-            cl_interp = (α, β) -> 0.5cl_left(α, β) + 0.5cl_right(α, β)
-            cd_interp = (α, β) -> 0.5cd_left(α, β) + 0.5cd_right(α, β)
-            cm_interp = (α, β) -> 0.5cm_left(α, β) + 0.5cm_right(α, β)
-            polar_data = (cl_interp, cd_interp, cm_interp)
+            cl_interp = cl_left === cl_right ? cl_left :
+                (α, β) -> 0.5cl_left(α, β) + 0.5cl_right(α, β)
+            cd_interp = cd_left === cd_right ? cd_left :
+                (α, β) -> 0.5cd_left(α, β) + 0.5cd_right(α, β)
+            cm_interp = cm_left === cm_right ? cm_left :
+                (α, β) -> 0.5cm_left(α, β) + 0.5cm_right(α, β)
         elseif aero_model != "inviscid"
             throw(ArgumentError("Unsupported aero model: $aero_model"))
         end
@@ -132,7 +124,8 @@ mutable struct Panel
         new(
             TE_point_1, LE_point_1, TE_point_2, LE_point_2,
             chord, nothing, corner_points, aero_model,
-            cl_coeffs, cd_coeffs, cm_coeffs, polar_data,
+            cl_coeffs, cd_coeffs, cm_coeffs,
+            cl_interp, cd_interp, cm_interp,
             aerodynamic_center, control_point,
             bound_point_1, bound_point_2,
             x_airf, y_airf, z_airf,
@@ -269,9 +262,9 @@ function calculate_cl(panel::Panel, alpha::Float64)
     elseif panel.aero_model == "inviscid"
         return 2π * alpha
     elseif panel.aero_model == "polar_data"
-        return panel.polar_data[1](alpha)
+        return panel.cl_interp(alpha)
     elseif panel.aero_model == "interpolations"
-        return panel.polar_data[1](alpha, 0.0)
+        return panel.cl_interp(alpha, 0.0)
     else
         throw(ArgumentError("Unsupported aero model: $(panel.aero_model)"))
     end
@@ -281,6 +274,14 @@ function calculate_cl(polar_data::Nothing, model_type::String, alpha::Float64)
     return 2π * alpha
 end
 
+function calculate_cl(polar_data::Tuple{Function, Function, Function}, model_type::String, alpha::Float64)
+    if model_type == "polar_data"
+        return polar_data[1](alpha)
+    elseif model_type == "interpolations"
+        return polar_data[1](alpha, 0.0)
+    end
+end
+
 """
     calculate_cd_cm(panel::Panel, alpha::Float64)
 
@@ -297,12 +298,12 @@ function calculate_cd_cm(panel::Panel, alpha::Float64)
     elseif panel.aero_model == "inviscid"
         return 0.0, 0.0
     elseif panel.aero_model == "polar_data"
-        cd = panel.polar_data[2](alpha)
-        cm = panel.polar_data[3](alpha)
+        cd = panel.cd_interp(alpha)
+        cm = panel.cm_interp(alpha)
         return cd, cm
     elseif panel.aero_model == "interpolations"
-        cd = panel.polar_data[2](alpha, 0.0)
-        cm = panel.polar_data[3](alpha, 0.0)
+        cd = panel.cd_interp(alpha, 0.0)
+        cm = panel.cm_interp(alpha, 0.0)
         return cd, cm
     else
         throw(ArgumentError("Unsupported aero model: $(panel.aero_model)"))

src/solver.jl

@@ -188,7 +188,7 @@ function gamma_loop(
     relaxation_factor::Float64
 )
     converged = false
-    n_panels = body_aero.n_panels
+    n_panels = length(body_aero.panels)
     alpha_array = body_aero.alpha_array
     v_a_array = body_aero.v_a_array
     Umagw_array = similar(v_a_array)
@@ -244,7 +244,7 @@ function gamma_loop(
         end
 
         for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array))
-            cl_array[i] = calculate_cl(panel.polar_data, panel.aero_model, alpha)
+            cl_array[i] = calculate_cl(panel, alpha)
         end
         gamma_new .= 0.5 .* v_a_array.^2 ./ Umagw_array .* cl_array .* chord_array
 

src/wing_geometry.jl

@@ -239,23 +239,14 @@ function calculate_new_aero_input(aero_input,
         cl_left, cd_left, cm_left = aero_input[section_index][2]
         cl_right, cd_right, cm_right = aero_input[section_index + 1][2]
 
-        @info "Interpolation object comparison" begin
-            objects_equal = (
-                cl_left === cl_right,
-                cd_left === cd_right,
-                cm_left === cm_right
-            )
-            (; objects_equal)
-        end
-
-        CL_interp = cl_left === cl_right ? cl_left :
+        cl_interp = cl_left === cl_right ? (α, β) -> cl_left[α, β] :
             (α, β) -> left_weight * cl_left[α, β] + right_weight * cl_right[α, β]
-        CD_interp = cd_left === cd_right ? cd_left :
+        cd_interp = cd_left === cd_right ? (α, β) -> cd_left[α, β] :
             (α, β) -> left_weight * cd_left[α, β] + right_weight * cd_right[α, β]
-        CM_interp = cm_left === cm_right ? cm_left :
+        cm_interp = cm_left === cm_right ? (α, β) -> cm_left[α, β] :
             (α, β) -> left_weight * cm_left[α, β] + right_weight * cm_right[α, β]
 
-        return ("interpolations", (CL_interp, CD_interp, CM_interp))
+        return ("interpolations", (cl_interp, cd_interp, cm_interp))
 
     elseif model_type == "lei_airfoil_breukels"
         tube_diameter_left = aero_input[section_index][2][1]

test/test_body_aerodynamics.jl

@@ -195,13 +195,13 @@ end
         va_norm_array = fill(norm(Uinf), length(coord))
         va_unit_array = repeat(reshape(Uinf ./ norm(Uinf), 1, 3), length(coord))
         calculate_AIC_matrices!(
-            wing_aero,
+            body_aero,
             "VSM",
             core_radius_fraction,
             va_norm_array,
             va_unit_array
         )
-        AIC_x, AIC_y, AIC_z = wing_aero.AIC[1, :, :], wing_aero.AIC[2, :, :], wing_aero.AIC[3, :, :]
+        AIC_x, AIC_y, AIC_z = body_aero.AIC[1, :, :], body_aero.AIC[2, :, :], body_aero.AIC[3, :, :]
 
         # Compare matrices with higher precision for VSM
         @test isapprox(MatrixU, AIC_x, atol=1e-8)

test/test_panel.jl

@@ -134,17 +134,6 @@ end
     # Create panel
     panel = create_panel(section1, section2)
 
-    @testset "Panel Polar Data Initialization" begin
-        # Test if panel has polar data
-        @test hasproperty(panel, :polar_data)
-        @test !isnothing(panel.polar_data)
-        @test size(panel.polar_data) == size(polar_data)
-        
-        # Test if panel polar data is correctly averaged
-        expected_data = (polar_data + big_polar_data) / 2
-        @test isapprox(panel.polar_data, expected_data, rtol=1e-5)
-    end
-    
     @testset "Coefficient Interpolation" begin
         test_alphas = deg2rad.([-5.0, 0.0, 5.0, 10.0])
         for alpha in test_alphas