plot_section_polars draws cl, cd and cm against alpha per panel, for every aero model

CHANGELOG.md

@@ -8,6 +8,9 @@
   doi:10.1088/1742-6596/2767/2/022068): each section gets a drag increment and a force
   along its span from the flow across it, in `solve!`, `solve` and `linearize`. Opt-in
   via `is_with_viscous_drag_correction` (default `false`) on the solver settings.
+- `plot_section_polars(body_aero; panels, alphas, delta)` draws cl, cd and cm against α
+  per panel through `calculate_cl`/`calculate_cd`/`calculate_cm`, for every aero model
+  and at flap deflection `delta`, in one figure instead of one coefficient per call.
 - `linearize` takes a `BodyAerodynamics` with more than one wing; `theta_idxs` and
   `delta_idxs` then run over the unrefined sections of all wings in order.
 

docs/src/private_functions.md

@@ -276,6 +276,7 @@ panel_contour
 panel_normal
 plate_hinge_local
 panel_plate_geometry
+panel_polar_curves
 PLATE_FACES
 Makie.plot!(ax, panel::VortexStepMethod.Panel)
 Makie.plot!(ax, body::VortexStepMethod.BodyAerodynamics)

examples/obj_to_yaml_kite.jl

@@ -80,10 +80,9 @@ if PLOT
     plot_geometry(body_aero, "Ram air kite (converted from .obj)"; is_show=true,
         view_elevation=15, view_azimuth=-120, use_tex=USE_TEX)
 
-    # Airfoils and per-section polars recovered from the converted geometry
+    # Airfoils and panel polars recovered from the converted geometry
     plot_airfoils(geometry_file; symmetric=true, is_show=true)
-    plot_section_polars(body_aero, :cl; is_show=true)
-    plot_section_polars(body_aero, :cd; is_show=true)
+    plot_section_polars(body_aero; panels=[1, 10], is_show=true)
 
     plot_polars([solver], [body_aero], ["VSM (NeuralFoil polars from .obj)"];
         angle_range=range(-5, 20, length=26), v_a=va,

ext/VortexStepMethodMakieExt.jl

@@ -1,7 +1,8 @@
 module VortexStepMethodMakieExt
 using MakieControlPlots.Makie, VortexStepMethod, LinearAlgebra, Statistics, DelimitedFiles
 import MakieControlPlots
-import VortexStepMethod: calculate_filaments_for_plotting
+import VortexStepMethod: calculate_filaments_for_plotting, calculate_cl, calculate_cd,
+    calculate_cm
 import VortexStepMethod: ObjAdapter, AirfoilAero
 
 export plot_geometry, plot_distribution, plot_polars, save_plot, show_plot,
@@ -1324,13 +1325,12 @@ function VortexStepMethod.plot_combined_analysis(
             xlabel="α [°]",
             ylabel="Cm")
 
-        cl_vals = [first_body.panels[1].cl_interp(a) for a in alphas]
-        cd_vals = [first_body.panels[1].cd_interp(a) for a in alphas]
-        cm_vals = [first_body.panels[1].cm_interp(a) for a in alphas]
+        panel = first_body.panels[1]
+        cl, cd, cm = panel_polar_curves([panel], alphas, panel.delta)
 
-        lines!(ax_cl_curve, alphas_deg, cl_vals; color=:blue, linewidth=2)
-        lines!(ax_cd_curve, alphas_deg, cd_vals; color=:red, linewidth=2)
-        lines!(ax_cm_curve, alphas_deg, cm_vals; color=:green, linewidth=2)
+        lines!(ax_cl_curve, alphas_deg, only(cl); color=:blue, linewidth=2)
+        lines!(ax_cd_curve, alphas_deg, only(cd); color=:red, linewidth=2)
+        lines!(ax_cm_curve, alphas_deg, only(cm); color=:green, linewidth=2)
     end
 
     # [2,1] Spanwise Distributions (3×3 grid)
@@ -1477,48 +1477,43 @@ function VortexStepMethod.plot_combined_analysis(
 end
 
 """
-    plot_section_polars(body_aero, coefficient=:cl; is_show=true,
-                        is_save=false, save_path=nothing, data_type=".png")
+    panel_polar_curves(panels, alphas, deltas) -> (cl, cd, cm)
+
+Lift, drag and moment coefficients of each of `panels` over `alphas` [rad], each panel
+at its flap deflection in `deltas` [rad] (or one shared deflection), as one vector per
+panel per coefficient.
+"""
+function panel_polar_curves(panels, alphas, deltas)
+    cl = collect.(eachrow(calculate_cl.(panels, alphas', deltas)))
+    cd = collect.(eachrow(calculate_cd.(panels, alphas', deltas)))
+    cm = collect.(eachrow(calculate_cm.(panels, alphas', deltas)))
+    return cl, cd, cm
+end
+
+"""
+    plot_section_polars(body_aero; kwargs...)
 
 Implementation of [`plot_section_polars`](@ref); rendered through `MakieControlPlots`.
 """
-function VortexStepMethod.plot_section_polars(body_aero::BodyAerodynamics,
-    coefficient::Symbol=:cl; is_show::Bool=true, is_save::Bool=false,
-    save_path=nothing, data_type::String=".png")
-
-    coefficient in (:cl, :cd, :cm) ||
-        throw(ArgumentError("coefficient must be :cl, :cd, or :cm, got :$coefficient"))
-    idx = coefficient === :cl ? 2 : coefficient === :cd ? 3 : 4
-    label = uppercasefirst(string(coefficient))
-
-    alphas_deg = nothing
-    series = Vector{Float64}[]
-    labels = String[]
-    for wing in body_aero.wings
-        for (s, section) in enumerate(wing.unrefined_sections)
-            section.aero_model == POLAR_VECTORS || continue
-            aero = section.aero_data
-            aero === nothing && continue
-            section_alphas = rad2deg.(aero[1])
-            if isnothing(alphas_deg)
-                alphas_deg = collect(section_alphas)
-            elseif length(section_alphas) != length(alphas_deg)
-                @warn "section $s has a different α grid; plotting against the first section's α"
-            end
-            push!(series, Float64.(aero[idx]))
-            push!(labels, "section $s")
-        end
-    end
-    isempty(series) && error("No POLAR_VECTORS sections found in body")
+function VortexStepMethod.plot_section_polars(body_aero::BodyAerodynamics;
+    panels=eachindex(body_aero.panels), alphas=deg2rad.(-20:0.5:30), delta=nothing,
+    is_show::Bool=true, is_save::Bool=false, save_path=nothing,
+    data_type::String=".png")
+
+    panel_indices = vcat(panels)
+    chosen_panels = body_aero.panels[panel_indices]
+    deltas = something.(delta, getproperty.(chosen_panels, :delta))
+    cl, cd, cm = panel_polar_curves(chosen_panels, alphas, deltas)
+    labels = ["panel $i ($(panel.aero_model))"
+              for (i, panel) in zip(panel_indices, chosen_panels)]
 
-    plt = MakieControlPlots.plot(alphas_deg, series;
-        xlabel="α [deg]", ylabel=label, title="$label per section",
-        labels=labels, disp=(is_show || is_save))
+    plt = MakieControlPlots.plotx(rad2deg.(alphas), cl, cd, cm;
+        xlabel="α [deg]", ylabels=["cl", "cd", "cm"], title="Section polars",
+        labels=[labels], disp=(is_show || is_save))
 
     if is_save && !isnothing(save_path)
         isdir(save_path) || mkpath(save_path)
-        MakieControlPlots.savefig(
-            joinpath(save_path, "section_polars_$(coefficient)$(data_type)"))
+        MakieControlPlots.savefig(joinpath(save_path, "section_polars$(data_type)"))
     end
     return plt
 end

src/VortexStepMethod.jl

@@ -174,17 +174,22 @@ in sequence.
 function plot_combined_analysis end
 
 """
-    plot_section_polars(body_aero::BodyAerodynamics, coefficient=:cl; kwargs...)
+    plot_section_polars(body_aero::BodyAerodynamics; panels=eachindex(body_aero.panels),
+                        alphas=deg2rad.(-20:0.5:30), delta=nothing, kwargs...)
 
-Plot one polar coefficient (`:cl`, `:cd`, or `:cm`) against angle of attack for
-every section of a wing using stored `POLAR_VECTORS` data. Rendered through
-`MakieControlPlots`.
+Plot the lift, drag and moment coefficients against angle of attack for the chosen
+`panels`, one curve per panel, as each panel's [`calculate_cl`](@ref),
+[`calculate_cd`](@ref) and [`calculate_cm`](@ref) evaluate them for its aero model.
+Rendered through `MakieControlPlots`; returns its plot object.
 
 # Arguments
 - `body_aero`: the [`BodyAerodynamics`](@ref) to plot
-- `coefficient`: `:cl`, `:cd`, or `:cm` (default: `:cl`)
 
 # Keyword arguments
+- `panels`: index or indices into `body_aero.panels` (default: all panels)
+- `alphas`: angles of attack [rad] (default: `deg2rad.(-20:0.5:30)`)
+- `delta`: flap deflection [rad] a `POLAR_MATRICES` panel is evaluated at
+  (default: `nothing`, each panel's own `delta`)
 - `is_show`: whether to display (default: `true`)
 - `is_save`: whether to save (default: `false`)
 - `save_path`: directory to save the figure (default: `nothing`)

test/plotting/test_plotting.jl

@@ -20,6 +20,7 @@ using MakieControlPlots
 CairoMakie.activate!()
 
 using VortexStepMethod
+using VortexStepMethod.AirfoilAero: lei_poly_coeffs
 using Test
 
 const makie_ext = Base.get_extension(VortexStepMethod, :VortexStepMethodMakieExt)
@@ -28,7 +29,7 @@ global ram_wing = ram_air_matrix_wing(; n_panels=20, n_sections=4,
                           alpha_range=deg2rad.(-1:1.0:1),
                           delta_range=deg2rad.(-1:1.0:1))
 
-function create_body_aero()
+function create_body_aero(; aero_model=INVISCID, aero_data=nothing)
     n_panels = 20          # Number of panels
     span = 20.0            # Wing span [m]
     chord = 1.0            # Chord length [m]
@@ -41,11 +42,11 @@ function create_body_aero()
     add_section!(wing,
         [0.0, span/2, 0.0],
         [chord, span/2, 0.0],
-        INVISCID)
+        aero_model, aero_data)
     add_section!(wing,
         [0.0, -span/2, 0.0],
         [chord, -span/2, 0.0],
-        INVISCID)
+        aero_model, aero_data)
 
     refine!(wing)
     body_aero = BodyAerodynamics([wing])
@@ -433,6 +434,75 @@ function create_body_aero_with_skin(; n_panels=4)
     return body_aero, n_node
 end
 
+"""
+The `(cl, cd, cm)` curves `plot_section_polars` drew for its `curve`-th panel.
+"""
+section_polar_curves(plt, curve) = Tuple(channel[curve] for channel in plt.Y)
+
+@testset "plot_section_polars draws cl, cd and cm per panel for every aero model" begin
+    alphas = deg2rad.(-2:1.0:2)
+    inviscid = create_body_aero()
+    inviscid_plt = plot_section_polars(inviscid; panels=[1, 20], alphas, is_show=false)
+
+    @testset "one curve per chosen panel, all panels by default" begin
+        @test inviscid_plt.X ≈ rad2deg.(alphas)
+        @test length(inviscid_plt.Y) == 3
+        @test all(length(channel) == 2 for channel in inviscid_plt.Y)
+        all_panels = plot_section_polars(inviscid; is_show=false)
+        @test length(all_panels.Y[1]) == length(inviscid.panels)
+    end
+
+    @testset "INVISCID panel is thin-airfoil lift without drag or moment" begin
+        cl, cd, cm = section_polar_curves(inviscid_plt, 2)
+        @test cl ≈ 2π .* alphas
+        @test all(iszero, cd)
+        @test all(iszero, cm)
+    end
+
+    @testset "POLAR_VECTORS panel follows its polar table" begin
+        vectors, _ = create_body_aero_with_skin()
+        plt = plot_section_polars(vectors; panels=1, alphas, is_show=false)
+        cl, cd, cm = section_polar_curves(plt, 1)
+        @test cl ≈ 0.5 .+ 0.25 .* rad2deg.(alphas)
+        @test cd ≈ fill(0.02, length(alphas))
+        @test cm ≈ fill(-0.05, length(alphas))
+    end
+
+    @testset "POLY panel follows its Breukels polynomials" begin
+        cl_coeffs, cd_coeffs, cm_coeffs = lei_poly_coeffs(2.0, 0.5)
+        poly = create_body_aero(; aero_model=POLY,
+                                aero_data=(cl_coeffs, cd_coeffs, cm_coeffs))
+        plt = plot_section_polars(poly; panels=2, alphas, is_show=false)
+        cl, cd, cm = section_polar_curves(plt, 1)
+        @test cl ≈ evalpoly.(rad2deg.(alphas), Ref(cl_coeffs))
+        @test cd ≈ evalpoly.(rad2deg.(alphas), Ref(cd_coeffs))
+        @test cm ≈ evalpoly.(rad2deg.(alphas), Ref(cm_coeffs))
+    end
+
+    @testset "POLAR_MATRICES panel is evaluated at the passed delta" begin
+        matrices = BodyAerodynamics([ram_wing])
+        panel = matrices.panels[3]
+        flap_alphas = deg2rad.(-1:0.5:1)
+        delta = deg2rad(1.0)
+        deflected = section_polar_curves(plot_section_polars(matrices; panels=3,
+            alphas=flap_alphas, delta, is_show=false), 1)
+        stored = section_polar_curves(plot_section_polars(matrices; panels=3,
+            alphas=flap_alphas, is_show=false), 1)
+        @test deflected[1] ≈ panel.cl_interp.(flap_alphas, delta)
+        @test deflected[2] ≈ panel.cd_interp.(flap_alphas, delta)
+        @test deflected[3] ≈ panel.cm_interp.(flap_alphas, delta)
+        @test stored[2] ≈ panel.cd_interp.(flap_alphas, panel.delta)
+        @test deflected[2] != stored[2]
+    end
+
+    @testset "is_save writes section_polars.png" begin
+        save_dir = mktempdir()
+        plot_section_polars(inviscid; panels=1, alphas, is_show=false, is_save=true,
+                            save_path=save_dir)
+        @test isfile(joinpath(save_dir, "section_polars.png"))
+    end
+end
+
 @testset "Airfoil skin (Makie)" begin
     airfoil_skin_geometry = getfield(makie_ext, :airfoil_skin_geometry)
     skin_observables = getfield(makie_ext, :AIRFOIL_SKIN_OBSERVABLES)