added alpha_geometric_dist to the solver output, and updated example V3_kite to also test solve! output

Author: jellepoland

examples/V3_kite.jl

@@ -1,6 +1,7 @@
 using LinearAlgebra
 using VortexStepMethod
 using GLMakie
+using DelimitedFiles
 
 PLOT = true
 USE_TEX = false
@@ -57,7 +58,7 @@ USE_TEX = false
 
 # Solve and plot combined analysis
 results = VortexStepMethod.solve(solver, body_aero; log=true)
-PLOT && plot_combined_analysis(
+fig1 = plot_combined_analysis(
     solver,
     body_aero,
     results;
@@ -69,10 +70,212 @@ PLOT && plot_combined_analysis(
     side_slip=sideslip_deg,
     v_a=wind_speed,
     title="TU Delft V3 Kite",
-    is_show=true,
+    is_show=false,
     use_tex=USE_TEX,
     angle_of_attack_for_spanwise_distribution=10.0,
 )
+scr1 = display(fig1)
+wait(scr1)
+
+
+# Polar sweep including CMy, using solve! and calculate_results
+function compute_polar_with_cm(
+    solver,
+    body_aero,
+    angle_range;
+    angle_type::String="angle_of_attack",
+    angle_of_attack::Float64=0.0,
+    side_slip::Float64=0.0,
+    v_a::Float64=10.0
+)
+    n_angles = length(angle_range)
+    cl = zeros(n_angles)
+    cd = zeros(n_angles)
+    cs = zeros(n_angles)
+    cmy = fill(NaN, n_angles)
+    reynolds_number = zeros(n_angles)
+
+    gamma_prev = solver.sol.gamma_distribution
+    for (i, angle_i) in enumerate(angle_range)
+        if angle_type == "angle_of_attack"
+            α = deg2rad(angle_i)
+            β = deg2rad(side_slip)
+        elseif angle_type == "side_slip"
+            α = deg2rad(angle_of_attack)
+            β = deg2rad(angle_i)
+        else
+            throw(ArgumentError("angle_type must be 'angle_of_attack' or 'side_slip'"))
+        end
+
+        set_va!(body_aero, [cos(α) * cos(β), sin(β), sin(α)] * v_a)
+        solve!(solver, body_aero, gamma_prev; log=false)
+        gamma_prev = solver.sol.gamma_distribution
+
+        results_local = calculate_results(
+            body_aero,
+            solver.lr.gamma_new,
+            zeros(MVec3),
+            solver.density,
+            solver.aerodynamic_model_type,
+            solver.core_radius_fraction,
+            solver.mu,
+            solver.lr.alpha_dist,
+            solver.lr.v_a_dist,
+            solver.sol._chord_dist,
+            solver.sol._x_airf_dist,
+            solver.sol._y_airf_dist,
+            solver.sol._z_airf_dist,
+            solver.sol._va_dist,
+            solver.br.va_norm_dist,
+            solver.br.va_unit_dist,
+            body_aero.panels,
+            solver.is_only_f_and_gamma_output;
+            correct_aoa=solver.correct_aoa
+        )
+
+        cl[i] = results_local["cl"]
+        cd[i] = results_local["cd"]
+        cs[i] = results_local["cs"]
+        cmy[i] = get(results_local, "cmy", NaN)
+        reynolds_number[i] = results_local["Rey"]
+    end
+
+    return (angle=angle_range, cl=cl, cd=cd, cs=cs, cmy=cmy, rey=reynolds_number)
+end
+
+function plot_polars_with_cmy(
+    solver_list,
+    body_aero_list,
+    label_list;
+    literature_path_list::Vector{String}=String[],
+    angle_range=range(-10, 40, step=1),
+    angle_type::String="angle_of_attack",
+    angle_of_attack::Float64=0.0,
+    side_slip::Float64=0.0,
+    v_a::Float64=10.0,
+    title::String="polar_with_cm",
+    fig_size::Tuple{Int,Int}=(1200, 800),
+    angle_xlim::Tuple{Real,Real}=(-10, 40)
+)
+    total_cases = length(body_aero_list) + length(literature_path_list)
+    length(label_list) == total_cases || throw(ArgumentError("labels length ($(length(label_list))) must match number of cases ($total_cases)"))
+    length(solver_list) == length(body_aero_list) || throw(ArgumentError("solver_list length must match body_aero_list length"))
+
+    polar_data_list = Vector{Any}()
+    labels_full = String[]
+
+    # Computational cases
+    for (solver_i, body, lbl) in zip(solver_list, body_aero_list, label_list[1:length(solver_list)])
+        pd = compute_polar_with_cm(
+            solver_i,
+            body,
+            angle_range;
+            angle_type=angle_type,
+            angle_of_attack=angle_of_attack,
+            side_slip=side_slip,
+            v_a=v_a
+        )
+        @info "polar sample (solver)" label=lbl first_cl=pd.cl[1] first_cd=pd.cd[1] first_cs=pd.cs[1] first_cmy=pd.cmy[1]
+        push!(polar_data_list, pd)
+        re_tag = round(Int, first(pd.rey) * 1e-5)
+        push!(labels_full, "$(lbl) Re=$(re_tag)e5")
+    end
+
+    # Literature cases
+    for (path, lbl) in zip(literature_path_list, label_list[length(solver_list)+1:end])
+        data = readdlm(path, ',')
+        header_raw = string.(data[1, :])
+        header = lowercase.(strip.(header_raw))
+        alpha_idx = findfirst(x -> occursin("alpha", x) || occursin("aoa", x), header)
+        cl_idx    = findfirst(x -> occursin("cl", x), header)
+        cd_idx    = findfirst(x -> occursin("cd", x), header)
+        cs_idx    = findfirst(x -> occursin("cs", x), header)
+        cmy_idx   = findfirst(x -> occursin("cmy", x), header)
+
+        parse_col(col) = begin
+            vals = Float64[]
+            for v in col
+                if v isa Real
+                    push!(vals, Float64(v))
+                else
+                    s = strip(String(v))
+                    y = tryparse(Float64, s)
+                    push!(vals, isnothing(y) ? NaN : y)
+                end
+            end
+            vals
+        end
+
+        alpha_col = parse_col(data[2:end, alpha_idx])
+        cl_col    = parse_col(data[2:end, cl_idx])
+        cd_col    = parse_col(data[2:end, cd_idx])
+        cs_col    = cs_idx === nothing ? zeros(size(data, 1)-1) : parse_col(data[2:end, cs_idx])
+        cmy_col   = cmy_idx === nothing ? nothing : parse_col(data[2:end, cmy_idx])
+
+        push!(polar_data_list, (angle=alpha_col, cl=cl_col, cd=cd_col, cs=cs_col, cmy=cmy_col, rey=fill(NaN, length(alpha_col))))
+        push!(labels_full, lbl)
+    end
+
+    fig = Figure(size=fig_size)
+    Label(fig[0, :], title; fontsize=20, font=:bold)
+
+    ax_cl    = Axis(fig[1, 1], title="CL vs $angle_type [deg]",  xlabel="$angle_type [deg]", ylabel="CL [-]")
+    ax_cd    = Axis(fig[1, 2], title="CD vs $angle_type [deg]",  xlabel="$angle_type [deg]", ylabel="CD [-]")
+    ax_cmy   = Axis(fig[1, 3], title="CMy vs $angle_type [deg]", xlabel="$angle_type [deg]", ylabel="CMy [-]")
+    ax_cs    = Axis(fig[2, 1], title="CS vs $angle_type [deg]",  xlabel="$angle_type [deg]", ylabel="CS [-]")
+    ax_polar = Axis(fig[2, 2], title="CL vs CD", xlabel="CD [-]", ylabel="CL [-]")
+
+    xlims!(ax_cl, angle_xlim...)
+    xlims!(ax_cd, angle_xlim...)
+    xlims!(ax_cmy, angle_xlim...)
+    xlims!(ax_cs, angle_xlim...)
+
+    colors = Makie.wong_colors()
+
+    for (idx, (pd, lbl)) in enumerate(zip(polar_data_list, labels_full))
+        color = colors[mod1(idx, length(colors))]
+
+        lines!(ax_cl, pd.angle, pd.cl; color, label=lbl)
+        scatter!(ax_cl, pd.angle, pd.cl; color, markersize=6)
+
+        lines!(ax_cd, pd.angle, pd.cd; color, label=lbl)
+        scatter!(ax_cd, pd.angle, pd.cd; color, markersize=6)
+
+        lines!(ax_cs, pd.angle, pd.cs; color, label=lbl)
+        scatter!(ax_cs, pd.angle, pd.cs; color, markersize=6)
+
+        lines!(ax_polar, pd.cd, pd.cl; color, label=lbl)
+        scatter!(ax_polar, pd.cd, pd.cl; color, markersize=6)
+
+        if pd.cmy !== nothing
+            cmy_vals = Float64.(pd.cmy)
+            if !all(isnan, cmy_vals)
+                lines!(ax_cmy, pd.angle, cmy_vals; color, label=lbl)
+                scatter!(ax_cmy, pd.angle, cmy_vals; color, markersize=6)
+            end
+        end
+    end
+
+    Legend(fig[2, 3], ax_cl)
+    return fig
+end
+
+
+fig2 = plot_polars_with_cmy(
+    [solver],
+    [body_aero],
+    labels;
+    literature_path_list=literature_paths,
+    angle_range=range(-5, 40, step=1),
+    angle_type="angle_of_attack",
+    angle_of_attack=angle_of_attack_deg,
+    side_slip=sideslip_deg,
+    v_a=wind_speed,
+    title="Testing solve! on V3 kite",
+    angle_xlim=(-5, 40)
+)
+scr2 = display(fig2)
+wait(scr2)
 
 
 nothing

src/solver.jl

@@ -43,6 +43,7 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all
     ### end of private vectors
     width_dist::Vector{Float64} = zeros(P)
     alpha_dist::Vector{Float64} = zeros(P)
+    alpha_geometric_dist::Vector{Float64} = zeros(P)
     cl_dist::Vector{Float64} = zeros(P)
     cd_dist::Vector{Float64} = zeros(P)
     cm_dist::Vector{Float64} = zeros(P)
@@ -208,6 +209,7 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
     converged = solver.lr.converged
     alpha_dist = solver.lr.alpha_dist
     alpha_corrected = solver.sol.alpha_dist
+    alpha_geometric_dist = solver.sol.alpha_geometric_dist
     v_a_dist = solver.lr.v_a_dist
     panels = body_aero.panels
 
@@ -224,6 +226,28 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
         cl_dist[i] = calculate_cl(panel, alpha_dist[i])
         cd_dist[i], cm_dist[i] = calculate_cd_cm(panel, alpha_dist[i])
         width_dist[i] = panel.width
+
+        # Geometric AoA using panel-local axes and prescribed freestream
+        # @views makes slices like solver.sol._va_dist[i, :] return views instead of copies. 
+        #   Without it, each [...] would allocate a new array; with it, you reuse a lightweight 
+        #   window into the original array, which cuts allocations in this tight loop.
+        @views begin
+            va_panel = solver.sol._va_dist[i, :]
+            x_airf = solver.sol._x_airf_dist[i, :]
+            z_airf = solver.sol._z_airf_dist[i, :]
+            va_norm = norm(va_panel)
+            x_norm = norm(x_airf)
+            z_norm = norm(z_airf)
+            if va_norm == 0 || x_norm == 0 || z_norm == 0
+                alpha_geometric_dist[i] = NaN
+            else
+                v_unit = va_panel / va_norm
+                v_tangential = dot(x_airf, -v_unit) / x_norm
+                v_normal = dot(z_airf, -v_unit) / z_norm
+                alpha_geometric_dist[i] = pi + atan(v_normal, v_tangential)
+            end
+        end
+
     end
 
     # create an alias for the three vertical output vectors
@@ -465,6 +489,10 @@ function solve(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=n
         solver.is_only_f_and_gamma_output;
         correct_aoa=solver.correct_aoa
     )
+    # Attach geometric AoA (already computed in calculate_results) to solver.sol
+    if haskey(results, "alpha_geometric")
+        solver.sol.alpha_geometric_dist .= results["alpha_geometric"]
+    end
     return results
 end
 
@@ -885,4 +913,3 @@ function linearize(solver::Solver, body_aero::BodyAerodynamics, y::Vector{T};
     calc_results!(results, y)
     return jac, results
 end
-