add results against output test

Author: 1-Bart-1

src/panel.jl

@@ -239,12 +239,12 @@ Calculate lift coefficient for given angle of attack.
 function calculate_cl(panel::Panel, alpha::Float64)
     if panel.panel_aero_model == "lei_airfoil_breukels"
         cl = evalpoly(rad2deg(alpha), reverse(panel.cl_coefficients))
-        if abs(alpha) > (π/9)  # Outside ±20 degrees
+        if abs(alpha) > (π/9)
             cl = 2 * cos(alpha) * sin(alpha)^2
         end
         return cl
     elseif panel.panel_aero_model == "inviscid"
-        return 2π * alpha  # Changed: 2 * π to 2π for more idiomatic Julia
+        return 2π * alpha
     elseif panel.panel_aero_model == "polar_data"
         return linear_interpolation(
             panel.panel_polar_data[:,1],

src/solver.jl

@@ -60,7 +60,7 @@ end
 Main solving routine for the aerodynamic model.
 """
 function solve(solver::Solver, wing_aero::WingAerodynamics, gamma_distribution=nothing)
-    isnothing(wing_aero.va) && throw(ArgumentError("Inflow conditions are not set"))
+    isnothing(wing_aero.panels[1].va) && throw(ArgumentError("Inflow conditions are not set, use set_va!(wing_aero, va)"))
 
     # Initialize variables
     panels = wing_aero.panels

test/test_wing_aerodynamics.jl

@@ -0,0 +1,102 @@
+# using VortexStepMethod: Wing, WingAerodynamics, BoundFilament, SemiInfiniteFilament, add_section!, set_va!, solve, calculate_cl
+using VortexStepMethod
+using VortexStepMethod: calculate_cl, calculate_cd_cm, calculate_projected_area
+using LinearAlgebra
+using Test
+using Logging
+
+include("utils.jl")
+
+@testset "Calculate results against output results" begin
+    # Setup
+    density = 1.225
+    N = 40
+    max_chord = 1.0
+    span = 15.709  # AR = 20
+    Umag = 20.0
+    AR = span^2 / (π * span * max_chord / 4)
+    aoa = deg2rad(5)
+    Uinf = [cos(aoa), 0.0, sin(aoa)] .* Umag
+    model = "VSM"
+
+    # Setup wing geometry
+    dist = "cos"
+    core_radius_fraction = 1e-20
+    coord = generate_coordinates_el_wing(max_chord, span, N, dist)
+    coord_left_to_right = flip_created_coord_in_pairs(deepcopy(coord))
+    wing = Wing(N; spanwise_panel_distribution="unchanged")
+    for idx in 1:2:length(coord_left_to_right[:, 1])
+        @debug "coord_left_to_right[$idx] = $(coord_left_to_right[idx,:])"
+        add_section!(
+            wing,
+            coord_left_to_right[idx,:],
+            coord_left_to_right[idx+1,:],
+            "inviscid"
+        )
+    end
+    
+    wing_aero = WingAerodynamics([wing])
+    set_va!(wing_aero, (Uinf, 0.0))
+
+    # Run analysis
+    solver_object = Solver(
+        aerodynamic_model_type=model, 
+        core_radius_fraction=core_radius_fraction
+    )
+    results_NEW = solve(solver_object, wing_aero)
+
+    @test results_NEW isa Dict
+
+    # Calculate forces using uncorrected alpha
+    alpha = results_NEW["alpha_uncorrected"]
+    dyn_visc = 0.5 * density * norm(Uinf)^2
+    n_panels = length(wing_aero.panels)
+    lift = zeros(n_panels)
+    drag = zeros(n_panels)
+    moment = zeros(n_panels)
+    
+    for (i, panel) in enumerate(wing_aero.panels)
+        lift[i] = dyn_visc * calculate_cl(panel, alpha[i]) * panel.chord
+        cd_cm = calculate_cd_cm(panel, alpha[i])
+        drag[i] = dyn_visc * cd_cm[1] * panel.chord
+        moment[i] = dyn_visc * cd_cm[2] * panel.chord^2
+        @info "lift: $lift, drag: $drag, moment: $moment"
+    end
+    Fmag = hcat(lift, drag, moment)
+
+    # Calculate coefficients using corrected alpha
+    alpha = results_NEW["alpha_at_ac"]
+    aero_coeffs = hcat(
+        [alpha[i] for (i, panel) in enumerate(wing_aero.panels)],
+        [calculate_cl(panel, alpha[i]) for (i, panel) in enumerate(wing_aero.panels)],
+        [calculate_cd_cm(panel, alpha[i])[1] for (i, panel) in enumerate(wing_aero.panels)],
+        [calculate_cd_cm(panel, alpha[i])[2] for (i, panel) in enumerate(wing_aero.panels)]
+    )
+    
+    ringvec = [Dict("r0" => panel.width * panel.z_airf) for panel in wing_aero.panels]
+    controlpoints = [Dict("tangential" => panel.y_airf, "normal" => panel.x_airf) 
+                    for panel in wing_aero.panels]
+    Atot = calculate_projected_area(wing)
+
+    F_rel_ref, F_gl_ref, Ltot_ref, Dtot_ref, CL_ref, CD_ref, CS_ref = 
+        output_results(Fmag, aero_coeffs, ringvec, Uinf, controlpoints, Atot)
+
+    # Compare results
+    @info "Comparing results"
+    @info "cl_calculated: $(results_NEW["cl"]), CL_ref: $CL_ref"
+    @info "cd_calculated: $(results_NEW["cd"]), CD_ref: $CD_ref"
+    @info "cs_calculated: $(results_NEW["cs"]), CS_ref: $CS_ref"
+    @info "L_calculated: $(results_NEW["lift"]), Ltot_ref: $Ltot_ref"
+    @info "D_calculated: $(results_NEW["drag"]), Dtot_ref: $Dtot_ref"
+
+    # Assert results
+    @test isapprox(results_NEW["cl"], CL_ref, rtol=1e-4)
+    @test isapprox(results_NEW["cd"], CD_ref, rtol=1e-4)
+    @test isapprox(results_NEW["cs"], CS_ref, rtol=1e-4)
+    @test isapprox(results_NEW["lift"], Ltot_ref, rtol=1e-4)
+    @test isapprox(results_NEW["drag"], Dtot_ref, rtol=1e-4)
+
+    # Check array shapes
+    @test length(results_NEW["cl_distribution"]) == length(wing_aero.panels)
+    @test length(results_NEW["cd_distribution"]) == length(wing_aero.panels)
+end

test/thesis_oriol_cayon.jl

@@ -0,0 +1,101 @@
+"""
+Post-process results to get global forces and aerodynamic coefficients
+
+Parameters
+----------
+Fmag : Lift, Drag and Moment magnitudes
+aero_coeffs : alpha, cl, cd, cm
+ringvec : List of dictionaries containing the vectors that define each ring
+Uinf : Wind speed velocity vector
+controlpoints : List of dictionaries with the variables needed to define each wing section
+Atot : Planform area
+
+Returns
+-------
+F_rel : Lift and drag forces relative to the local angle of attack
+F_gl : Lift and drag forces relative to the wind direction
+Ltot : Total lift
+Dtot : Total drag
+CL : Global CL
+CD : Global CD
+"""
+function output_results(Fmag, aero_coeffs, ringvec, Uinf, controlpoints, Atot)
+    rho = 1.225
+    alpha = aero_coeffs[:, 1]
+    F_rel = []
+    F_gl = []
+    Fmag_gl = []
+    SideF = []
+    Ltot = 0.0
+    Dtot = 0.0
+    SFtot = 0.0
+    
+    for i in eachindex(alpha)
+        r0 = ringvec[i]["r0"]
+        # Relative wind speed direction
+        dir_urel = cos(alpha[i]) * controlpoints[i]["tangential"] +
+                   sin(alpha[i]) * controlpoints[i]["normal"]
+        dir_urel = dir_urel / norm(dir_urel)
+        
+        # Lift direction relative to Urel
+        dir_L = cross(dir_urel, r0)
+        dir_L = dir_L / norm(dir_L)
+        
+        # Drag direction relative to Urel
+        dir_D = cross([0.0, 1.0, 0.0], dir_L)
+        dir_D = dir_D / norm(dir_D)
+        
+        # Lift and drag relative to Urel
+        L_rel = dir_L * Fmag[i, 1]
+        D_rel = dir_D * Fmag[i, 2]
+        push!(F_rel, [L_rel, D_rel])
+        
+        # Lift direction relative to the wind speed
+        dir_L_gl = cross(Uinf, [0.0, 1.0, 0.0])
+        dir_L_gl = dir_L_gl / norm(dir_L_gl)
+        
+        # Lift and drag relative to the windspeed
+        L_gl = vector_projection(L_rel, dir_L_gl) + vector_projection(D_rel, dir_L_gl)
+        D_gl = vector_projection(L_rel, Uinf) + vector_projection(D_rel, Uinf)
+        push!(F_gl, [L_gl, D_gl])
+        
+        push!(Fmag_gl, [
+            dot(L_rel, dir_L_gl) + dot(D_rel, dir_L_gl),
+            dot(L_rel, Uinf / norm(Uinf)) + dot(D_rel, Uinf / norm(Uinf))
+        ])
+        
+        push!(SideF, dot(L_rel, [0.0, 1.0, 0.0]) + dot(D_rel, [0.0, 1.0, 0.0]))
+    end
+
+    # Calculate total aerodynamic forces
+    for i in eachindex(Fmag_gl)
+        Ltot += Fmag_gl[i][1] * norm(ringvec[i]["r0"])
+        Dtot += Fmag_gl[i][2] * norm(ringvec[i]["r0"])
+        SFtot += SideF[i] * norm(ringvec[i]["r0"])
+    end
+
+    Umag = norm(Uinf)
+    CL = Ltot / (0.5 * Umag^2 * Atot * rho)
+    CD = Dtot / (0.5 * Umag^2 * Atot * rho)
+    CS = SFtot / (0.5 * Umag^2 * Atot * rho)
+
+    return F_rel, F_gl, Ltot, Dtot, CL, CD, CS
+end
+
+"""
+Find the projection of a vector into a direction
+
+Parameters
+----------
+v : vector to be projected
+u : direction
+
+Returns
+-------
+proj : projection of the vector v onto u
+"""
+function vector_projection(v, u)
+    unit_u = u / norm(u)
+    proj = dot(v, unit_u) * unit_u
+    return proj
+end