diff --git a/CHANGELOG.md b/CHANGELOG.md index 3c2babbf..a5656a49 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/docs/src/private_functions.md b/docs/src/private_functions.md index 5f7eec56..eb94971d 100644 --- a/docs/src/private_functions.md +++ b/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) diff --git a/examples/obj_to_yaml_kite.jl b/examples/obj_to_yaml_kite.jl index 91eff2f9..e756b2e2 100644 --- a/examples/obj_to_yaml_kite.jl +++ b/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, diff --git a/ext/VortexStepMethodMakieExt.jl b/ext/VortexStepMethodMakieExt.jl index f1227f8c..3cd2569a 100644 --- a/ext/VortexStepMethodMakieExt.jl +++ b/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 diff --git a/src/VortexStepMethod.jl b/src/VortexStepMethod.jl index 0521ddfc..f28a73f5 100644 --- a/src/VortexStepMethod.jl +++ b/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`) diff --git a/test/plotting/test_plotting.jl b/test/plotting/test_plotting.jl index d3ffc9cd..669bf42f 100644 --- a/test/plotting/test_plotting.jl +++ b/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)