Add result tests

Author: 1-Bart-1

examples/ram_air_kite.jl

@@ -49,7 +49,7 @@ if LINEARIZE
         body_aero, 
         wing, 
         [zeros(4); vel_app; zeros(3)]; 
-        alpha_idxs=1:4, 
+        theta_idxs=1:4, 
         va_idxs=5:7, 
         omega_idxs=8:10,
         moment_frac=0.1)
@@ -58,7 +58,7 @@ if LINEARIZE
         body_aero, 
         wing, 
         [zeros(4); vel_app; zeros(3)]; 
-        alpha_idxs=1:4, 
+        theta_idxs=1:4, 
         va_idxs=5:7, 
         omega_idxs=8:10,
         moment_frac=0.1)

src/body_aerodynamics.jl

@@ -48,13 +48,17 @@ Construct a [BodyAerodynamics](@ref) object for aerodynamic calculations.
 
 # Keyword Arguments
 - `kite_body_origin=zeros(MVec3)`: Origin point of kite body reference frame in CAD reference frame
+- `va=[15.0, 0.0, 0.0]`: Apparent wind vector
+- `omega=zeros(3)`: Turn rate in kite body frame x y and z
 
 # Returns
 - [BodyAerodynamics](@ref) object initialized with panels and wings
 """
 function BodyAerodynamics(
     wings::Vector{T};
-    kite_body_origin=zeros(MVec3)
+    kite_body_origin=zeros(MVec3),
+    va=[15.0, 0.0, 0.0],
+    omega=zeros(MVec3)
 ) where T <: AbstractWing
     # Initialize panels
     panels = Panel[]
@@ -78,7 +82,7 @@ function BodyAerodynamics(
     end
 
     body_aero = BodyAerodynamics{length(panels)}(; panels, wings)
-    init!(body_aero)
+    init!(body_aero; va, omega)
     return body_aero
 end
 
@@ -107,11 +111,17 @@ Initialize a BodyAerodynamics struct in-place by setting up panels and coefficie
 
 # Keyword Arguments
 - `init_aero::Bool`: Wether to initialize the aero data or not
+- `va=[15.0, 0.0, 0.0]`: Apparent wind vector
+- `omega=zeros(3)`: Turn rate in kite body frame x y and z
 
 # Returns
 nothing
 """
-function init!(body_aero::BodyAerodynamics; init_aero=true)
+function init!(body_aero::BodyAerodynamics; 
+    init_aero=true,
+    va=[15.0, 0.0, 0.0],
+    omega=zeros(MVec3)
+)
     idx = 1
     vec = zeros(MVec3)
     for wing in body_aero.wings
@@ -151,6 +161,7 @@ function init!(body_aero::BodyAerodynamics; init_aero=true)
     body_aero.alpha_array .= 0.0
     body_aero.v_a_array .= 0.0 
     body_aero.AIC .= 0.0
+    set_va!(body_aero, va, omega)
     return nothing
 end
 

src/kite_geometry.jl

@@ -368,7 +368,8 @@ Represents a curved wing that inherits from Wing with additional geometric prope
 
 # Fields
 - All fields from Wing:
-  - `n_panels::Int64`: Number of panels in aerodynamic mesh
+  - `n_panels::Int16`: Number of panels in aerodynamic mesh
+  - `n_groups::Int16`: Number of panel groups, each panel group has it's own twisting moment
   - `spanwise_distribution`::PanelDistribution: see: [PanelDistribution](@ref)
   - `spanwise_direction::MVec3`: Wing span direction vector
   - `sections::Vector{Section}`: List of wing sections, see: [Section](@ref)
@@ -387,7 +388,8 @@ Represents a curved wing that inherits from Wing with additional geometric prope
 
 """
 mutable struct RamAirWing <: AbstractWing
-    n_panels::Int64
+    n_panels::Int16
+    n_groups::Int16
     spanwise_distribution::PanelDistribution
     panel_props::PanelProperties
     spanwise_direction::MVec3
@@ -422,19 +424,19 @@ Constructor for a [RamAirWing](@ref) that allows to use an `.obj` and a `.dat` f
 - dat_path: Path to the `.dat` file, a standard format for 2d foil geometry
 
 # Keyword Parameters
-- alpha=0.0: Angle of attack of each segment relative to the x axis [rad]
-- crease_frac=0.75: The x coordinate around which the trailing edge rotates on a normalized 2d foil, 
+- crease_frac=0.9: The x coordinate around which the trailing edge rotates on a normalized 2d foil, 
                     used in the xfoil polar generation
 - wind_vel=10.0: Apparent wind speed in m/s, used in the xfoil polar generation
 - mass=1.0: Mass of the wing in kg, used for the inertia calculations 
 - `n_panels`=56: Number of panels.
-- `n_sections`=n_panels+1: Number of sections (there is a section on each side of each panel.)
+- `n_sections=n_panels+1`: Number of sections (there is a section on each side of each panel.)
+- `n_groups=n_panels ÷ 4`: Number of panel groups
 - `spanwise_distribution`=UNCHANGED: see: [PanelDistribution](@ref)
 - `spanwise_direction`=[0.0, 1.0, 0.0]
 - `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not
 """
-function RamAirWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10., mass=1.0, 
-                  n_panels=56, n_sections=n_panels+1, spanwise_distribution=UNCHANGED, 
+function RamAirWing(obj_path, dat_path; crease_frac=0.9, wind_vel=10., mass=1.0, 
+                  n_panels=56, n_sections=n_panels+1, n_groups=n_panels÷4, spanwise_distribution=UNCHANGED, 
                   spanwise_direction=[0.0, 1.0, 0.0], remove_nan=true, align_to_principal=false,
                   alpha_range=deg2rad.(-5:1:20), delta_range=deg2rad.(-5:1:20), interp_steps=n_sections
                   )
@@ -508,7 +510,7 @@ function RamAirWing(obj_path, dat_path; alpha=0.0, crease_frac=0.75, wind_vel=10
         panel_props = PanelProperties{n_panels}()
         cache = [LazyBufferCache()]
 
-        RamAirWing(n_panels, spanwise_distribution, panel_props, spanwise_direction, sections, 
+        RamAirWing(n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, sections, 
             refined_sections, remove_nan, non_deformed_sections,
             mass, gamma_tip, inertia_tensor, center_of_mass, radius,
             le_interp, te_interp, area_interp, zeros(n_panels), zeros(n_panels), cache)

src/wing_geometry.jl

@@ -194,7 +194,8 @@ end
 Represents a wing composed of multiple sections with aerodynamic properties.
 
 # Fields
-- `n_panels::Int64`: Number of panels in aerodynamic mesh
+- `n_panels::Int16`: Number of panels in aerodynamic mesh
+- `n_groups::Int16`: Number of panel groups
 - `spanwise_distribution`::PanelDistribution: [PanelDistribution](@ref)
 - `spanwise_direction::MVec3`: Wing span direction vector
 - `sections::Vector{Section}`: Vector of wing sections, see: [Section](@ref)
@@ -203,7 +204,8 @@ Represents a wing composed of multiple sections with aerodynamic properties.
 
 """
 mutable struct Wing <: AbstractWing
-    n_panels::Int64
+    n_panels::Int16
+    n_groups::Int16
     spanwise_distribution::PanelDistribution
     panel_props::PanelProperties
     spanwise_direction::MVec3
@@ -214,6 +216,7 @@ end
 
 """
     Wing(n_panels::Int;
+         n_groups=n_panels,
          spanwise_distribution::PanelDistribution=LINEAR,
          spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]),
          remove_nan::Bool=true)
@@ -222,17 +225,19 @@ Constructor for a [Wing](@ref) struct with default values that initializes the s
 and refined sections as empty arrays.
 
 # Parameters
-- `n_panels::Int64`: Number of panels in aerodynamic mesh
+- `n_panels::Int`: Number of panels in aerodynamic mesh
+- `n_groups::Int`: Number of panel groups in aerodynamic mesh
 - `spanwise_distribution`::PanelDistribution = LINEAR: [PanelDistribution](@ref)
 - `spanwise_direction::MVec3` = MVec3([0.0, 1.0, 0.0]): Wing span direction vector
 - `remove_nan::Bool`: Wether to remove the NaNs from interpolations or not
 """
 function Wing(n_panels::Int;
+        n_groups = n_panels,
         spanwise_distribution::PanelDistribution=LINEAR,
         spanwise_direction::PosVector=MVec3([0.0, 1.0, 0.0]),
         remove_nan=true)
     panel_props = PanelProperties{n_panels}()
-    Wing(n_panels, spanwise_distribution, panel_props, spanwise_direction, Section[], Section[], remove_nan)
+    Wing(n_panels, n_groups, spanwise_distribution, panel_props, spanwise_direction, Section[], Section[], remove_nan)
 end
 
 function init!(wing::AbstractWing)

test/runtests.jl

@@ -19,13 +19,19 @@ end::Bool
 cd("..")
 println("Running tests...")
 @testset verbose = true "Testing VortexStepMethod..." begin
+
+    cp("data/ram_air_kite_body.obj", "/tmp/ram_air_kite_body.obj"; force=true)
+    cp("data/ram_air_kite_foil.dat", "/tmp/ram_air_kite_foil.dat"; force=true)
+    ram_wing = RamAirWing("/tmp/ram_air_kite_body.obj", "/tmp/ram_air_kite_foil.dat"; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1))
+
     if build_is_production_build
         include("bench.jl")
     end
     include("test_bound_filament.jl")
     include("test_panel.jl")
     include("test_semi_infinite_filament.jl")
     include("test_body_aerodynamics.jl")
+    include("test_results.jl")
     include("test_kite_geometry.jl")
     include("test_wing_geometry.jl")
     include("test_plotting.jl")

test/test_body_aerodynamics.jl

@@ -113,140 +113,6 @@ include("utils.jl")
 end
 
 
-@testset "Calculate results against output results" begin
-    # Setup
-    density = 1.225
-    N = 40
-    max_chord = 1.0
-    span = 15.709  # AR = 20
-    v_a = 20.0
-    AR = span^2 / (π * span * max_chord / 4)
-    aoa = deg2rad(5)
-    v_a = [cos(aoa), 0.0, sin(aoa)] .* v_a
-    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_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
-    
-    body_aero = BodyAerodynamics([wing])
-    set_va!(body_aero, v_a)
-
-    # Run analysis
-    P = length(body_aero.panels)
-    loop_solver = Solver{P}(
-        aerodynamic_model_type=model, 
-        core_radius_fraction=core_radius_fraction,
-        solver_type=LOOP,
-        atol=1e-8,
-        rtol=1e-8
-    )
-    nonlin_solver = Solver{P}(
-        aerodynamic_model_type=model, 
-        core_radius_fraction=core_radius_fraction,
-        solver_type=NONLIN,
-        atol=1e-8,
-        rtol=1e-8
-    )
-    results_NEW = solve(loop_solver, body_aero; reference_point=[0,1,0])
-    # println(results_NEW)
-
-    @test results_NEW isa Dict
-
-    @testset "Loop and nonlin solve!" begin
-        loop_sol = solve!(loop_solver, body_aero; reference_point=[0,1,0])
-        nonlin_sol = solve!(nonlin_solver, body_aero; reference_point=[0,1,0])
-
-        @test all(isapprox.(nonlin_sol.gamma_distribution, loop_sol.gamma_distribution; atol=1e-4))
-
-        @test loop_sol.force.x ≈ -117.96518414816444 atol=1e-4
-        @test loop_sol.force.y ≈ 0.0 atol=1e-10
-        @test loop_sol.force.z ≈ 1481.996390329679 atol=1e-4 rtol= 1e-4
-
-        @test loop_sol.moment.x ≈ -1481.996390329678 atol=1e-4 rtol= 1e-4
-        @test loop_sol.moment.y ≈ 0.0 atol=1e-10
-        @test loop_sol.moment.z ≈ -117.9651841481644 atol=1e-4
-
-        @test loop_sol.force_coefficients[1] ≈ -0.039050322560956294 atol=1e-4 # CFx
-        @test loop_sol.force_coefficients[2] ≈ 0.0                   atol=1e-4 # CFy
-        @test loop_sol.force_coefficients[3] ≈ 0.49055973654418716   atol=1e-4 # CFz
-        @test loop_sol.force_coefficients[3] / loop_sol.force_coefficients[1] ≈ loop_sol.force[3] / loop_sol.force[1]
-        @test loop_sol.moment_distribution[1] ≈ -0.0006683569356186426 atol=1e-8
-        @test loop_sol.moment_coeff_dist[1] ≈ -2.212405554436003e-7 atol=1e-10
-        @test loop_sol.moment_distribution[1] / loop_sol.moment_distribution[2] ≈ loop_sol.moment_coeff_dist[1] / loop_sol.moment_coeff_dist[2]
-
-        @test loop_sol.solver_status == FEASIBLE
-    end
-
-    # Calculate forces using uncorrected alpha
-    alpha = results_NEW["alpha_uncorrected"]
-    dyn_visc = 0.5 * density * norm(v_a)^2
-    n_panels = length(body_aero.panels)
-    lift = zeros(n_panels)
-    drag = zeros(n_panels)
-    moment = zeros(n_panels)
-    
-    for (i, panel) in enumerate(body_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(body_aero.panels)],
-        [calculate_cl(panel, alpha[i]) for (i, panel) in enumerate(body_aero.panels)],
-        [calculate_cd_cm(panel, alpha[i])[1] for (i, panel) in enumerate(body_aero.panels)],
-        [calculate_cd_cm(panel, alpha[i])[2] for (i, panel) in enumerate(body_aero.panels)]
-    )
-    
-    ringvec = [Dict("r0" => panel.width * panel.y_airf) for panel in body_aero.panels]
-    controlpoints = [Dict("tangential" => panel.x_airf, "normal" => panel.z_airf) 
-                    for panel in body_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, v_a, 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)
-    @test isapprox(results_NEW["Fx"], results_NEW["Mz"], rtol=1e-4) # 1 meter arm
-    @test isapprox(results_NEW["My"], 0.0, atol=1e-3)
-    @test isapprox(results_NEW["Fz"], -results_NEW["Mx"], rtol=1e-4) # 1 meter arm
-
-    # Check array shapes
-    @test length(results_NEW["cl_distribution"]) == length(body_aero.panels)
-    @test length(results_NEW["cd_distribution"]) == length(body_aero.panels)
-end
-
 @testset "Wing Geometry Creation" begin
     @testset "Origin Translation" begin
         # Create minimal wing with three sections (2 panels)

test/test_plotting.jl

@@ -103,10 +103,7 @@ plt.ioff()
         rm("/tmp/Rectangular_Wing_Polars.pdf")
 
         # Step 9: Test polar data plotting
-        cp("data/ram_air_kite_body.obj", "/tmp/ram_air_kite_body.obj"; force=true)
-        cp("data/ram_air_kite_foil.dat", "/tmp/ram_air_kite_foil.dat"; force=true)
-        wing = RamAirWing("/tmp/ram_air_kite_body.obj", "/tmp/ram_air_kite_foil.dat"; alpha_range=deg2rad.(-1:1), delta_range=deg2rad.(-1:1))
-        body_aero = BodyAerodynamics([wing])
+        body_aero = BodyAerodynamics([ram_wing])
         fig = plot_polar_data(body_aero; is_show=false)
         @test fig isa plt.PyPlot.Figure
     end