Fix moment and force coefficients (#110)

* correct coefficient calculations and add moment coefficient distribution

* add tests to check solve

* fix com test and throw error on wrong com

* fix tests

Author: 1-Bart-1

src/kite_geometry.jl

@@ -175,7 +175,7 @@ function center_to_com!(vertices, faces)
     end
 
     com = com / area_total
-    !(abs(com[2]) < 0.01) && @warn "Center of mass $com has to lie on the xz-plane."
+    !(abs(com[2]) < 0.01) && throw(ArgumentError("Center of mass $com of .obj file has to lie on the xz-plane."))
     @info "Centering vertices of .obj file to the center of mass: $com"
     for v in vertices
         v .-= com

src/solver.jl

@@ -7,20 +7,24 @@ Struct for storing the solution of the [solve!](@ref) function. Must contain all
 # Attributes
 - aero_force::MVec3: Aerodynamic force vector in KB reference frame [N]
 - aero_moments::MVec3: Aerodynamic moments [Mx, My, Mz] around the reference point [Nm]
-- force_coefficients::MVec3: Aerodynamic force coefficients [CL, CD, CS] [-]
+- force_coefficients::MVec3: Aerodynamic force coefficients [CFx, CFy, CFz] [-]
+- moment_coefficients::MVec3: Aerodynamic moment coefficients [CMx, CMy, CMz] [-]
 - moment_distribution::Vector{Float64}: Pitching moments around the spanwise vector of each panel. [Nm]
+- moment_coefficient_distribution::Vector{Float64}: Pitching moment coefficient around the spanwise vector of each panel. [-]
 - solver_status::SolverStatus: enum, see [SolverStatus](@ref)
 """
 mutable struct VSMSolution    
     aero_force::MVec3          
     aero_moments::MVec3       
     force_coefficients::MVec3  
+    moment_coefficients::MVec3  
     moment_distribution::Vector{Float64}
+    moment_coefficient_distribution::Vector{Float64}
     solver_status::SolverStatus 
 end
 
 function VSMSolution()
-    VSMSolution(zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(3), FAILURE)
+    VSMSolution(zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(MVec3), zeros(3), zeros(3), FAILURE)
 end
 
 """
@@ -112,7 +116,9 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
     cm_array = zeros(n_panels)
     panel_width_array = zeros(n_panels)
     solver.sol.moment_distribution = zeros(n_panels)
+    solver.sol.moment_coefficient_distribution = zeros(n_panels)
     moment_distribution = solver.sol.moment_distribution
+    moment_coefficient_distribution = solver.sol.moment_coefficient_distribution
 
     # Calculate coefficients for each panel
     for (i, panel) in enumerate(panels)                                               # zero bytes
@@ -163,6 +169,9 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
     va = body_aero.va
     va_unit = va / va_mag
     q_inf = 0.5 * density * va_mag^2
+
+    # Calculate wing geometry properties
+    projected_area = sum(wing -> calculate_projected_area(wing), body_aero.wings)
 
     for (i, panel) in enumerate(panels)                                               # 30625 bytes
 
@@ -188,25 +197,25 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
         drag_induced_va = drag[i] * dir_drag_induced_va
         ftotal_induced_va = lift_induced_va + drag_induced_va
 
-        # Calculate forces in prescribed wing frame
-        dir_lift_prescribed_va = cross(va, spanwise_direction)
-        dir_lift_prescribed_va = dir_lift_prescribed_va / norm(dir_lift_prescribed_va)
+        # # Calculate forces in prescribed wing frame
+        # dir_lift_prescribed_va = cross(va, spanwise_direction)
+        # dir_lift_prescribed_va = dir_lift_prescribed_va / norm(dir_lift_prescribed_va)
 
-        # Calculate force components
-        lift_prescribed_va = dot(lift_induced_va, dir_lift_prescribed_va) + 
-                           dot(drag_induced_va, dir_lift_prescribed_va)
-        drag_prescribed_va = dot(lift_induced_va, va_unit) + 
-                           dot(drag_induced_va, va_unit)
-        side_prescribed_va = dot(lift_induced_va, spanwise_direction) + 
-                           dot(drag_induced_va, spanwise_direction)
+        # # Calculate force components
+        # lift_prescribed_va = dot(lift_induced_va, dir_lift_prescribed_va) + 
+        #                    dot(drag_induced_va, dir_lift_prescribed_va)
+        # drag_prescribed_va = dot(lift_induced_va, va_unit) + 
+        #                    dot(drag_induced_va, va_unit)
+        # side_prescribed_va = dot(lift_induced_va, spanwise_direction) + 
+        #                    dot(drag_induced_va, spanwise_direction)
 
         # Body frame forces
         f_body_3D[:,i] .= ftotal_induced_va .* panel.width
 
-        # Update sums
-        lift_wing_3D_sum += lift_prescribed_va * panel.width
-        drag_wing_3D_sum += drag_prescribed_va * panel.width  
-        side_wing_3D_sum += side_prescribed_va * panel.width
+        # # Update sums
+        # lift_wing_3D_sum += lift_prescribed_va * panel.width
+        # drag_wing_3D_sum += drag_prescribed_va * panel.width  
+        # side_wing_3D_sum += side_prescribed_va * panel.width
 
         # Calculate the moments
         # (1) Panel aerodynamic center in body frame:
@@ -226,26 +235,22 @@ function solve!(solver::Solver, body_aero::BodyAerodynamics, gamma_distribution=
         # Calculate the moment distribution (moment on each panel)
         arm = (moment_frac - 0.25) * panel.chord
         moment_distribution[i] = dot(ftotal_induced_va, panel.z_airf) * arm
+        moment_coefficient_distribution[i] = moment_distribution[i] ./ (q_inf * projected_area)
     end
 
-    # Calculate wing geometry properties
-    projected_area = sum(wing -> calculate_projected_area(wing), body_aero.wings)
-
-    Fx = sum(f_body_3D[1,:])
-    Fy = sum(f_body_3D[2,:])
-    Fz = sum(f_body_3D[3,:])
-    Mx = sum(m_body_3D[1,:])
-    My = sum(m_body_3D[2,:])
-    Mz = sum(m_body_3D[3,:])
-
-    CL = lift_wing_3D_sum / (q_inf * projected_area)
-    CD = drag_wing_3D_sum / (q_inf * projected_area)
-    CS = side_wing_3D_sum / (q_inf * projected_area)
-
     # update the result struct
-    solver.sol.aero_force .= MVec3([Fx, Fy, Fz])
-    solver.sol.aero_moments .= MVec3([Mx, My, Mz])
-    solver.sol.force_coefficients .= MVec3([CL, CD, CS])
+    solver.sol.aero_force .= [
+        sum(f_body_3D[1,:]),
+        sum(f_body_3D[2,:]),
+        sum(f_body_3D[3,:])
+    ]
+    solver.sol.aero_moments .= [
+        sum(m_body_3D[1,:]),
+        sum(m_body_3D[2,:]),
+        sum(m_body_3D[3,:])
+    ]
+    solver.sol.force_coefficients .= solver.sol.aero_force ./ (q_inf * projected_area)
+    solver.sol.moment_coefficients .= solver.sol.aero_moments ./ (q_inf * projected_area)
     if converged
         # TODO: Check if the result if feasible if converged
         solver.sol.solver_status = FEASIBLE

test/test_body_aerodynamics.jl

@@ -164,10 +164,13 @@ end
     @test sol.aero_moments.y ≈ 0.0 atol=1e-10
     @test sol.aero_moments.z ≈ -117.97225244011435
 
-    @test sol.force_coefficients[1] ≈ 0.4920964685099385    # CL
-    @test sol.force_coefficients[2] ≈ 0.0038533739066069946 # CD
-    @test sol.force_coefficients[3] ≈ 0.0 atol=1e-10        # CS
+    @test sol.force_coefficients[1] ≈ -0.039050322560956294 # CFx
+    @test sol.force_coefficients[2] ≈ 0.0                   # CFy
+    @test sol.force_coefficients[3] ≈ 0.49055973654418716   # CFz
+    @test sol.force_coefficients[3] / sol.force_coefficients[1] ≈ sol.aero_force[3] / sol.aero_force[1]
     @test sol.moment_distribution[1] ≈ -0.026242454074087297
+    @test sol.moment_coefficient_distribution[1] ≈ -8.686587525353832e-6
+    @test sol.moment_distribution[1] / sol.moment_distribution[2] ≈ sol.moment_coefficient_distribution[1] / sol.moment_coefficient_distribution[2]
 
     @test sol.solver_status == FEASIBLE
 

test/test_kite_geometry.jl

@@ -32,11 +32,11 @@ using Serialization
     end
 
     @testset "Center of Mass Calculation" begin
-        vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]]
+        vertices = [[0.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 0.0, 1.0]]
         faces = [[1, 2, 3]]
 
         com = center_to_com!(vertices, faces)
-        expected_com = [1/3, 1/3, 0.0]
+        expected_com = [1/3, 0.0, 1/3]
 
         @test isapprox(com, expected_com, rtol=1e-5)
     end
@@ -75,9 +75,15 @@ using Serialization
     @testset "Interpolation Creation" begin
         vertices = []
         z_center = 2.0
-        for θ in range(-π/4, π/4, length=10)
-            push!(vertices, [0.0, r*sin(θ), z_center + r*cos(θ)])
-            push!(vertices, [1.0, r*sin(θ), z_center + r*cos(θ)])
+        Δθ = π/2 / 1000
+        for θ in range(-π/4, π/4-Δθ, 1000)
+            frac = 1.0 - abs(4θ/π)  # Goes from 0 to 1 to 0
+            x_le = 0.1 - 0.1 * frac      # Goes from 0.1 to 0.0 to 0.1
+            x_te = 0.9 + 0.1 * frac # Goes from 0.9 to 1.0 to 0.9
+        
+            push!(vertices, [x_le, r*sin(θ), z_center + r*cos(θ)])
+            push!(vertices, [x_te, r*sin(θ+0.5Δθ), z_center + r*cos(θ+0.5Δθ)])
+            push!(vertices, [x_le, r*sin(θ+Δθ), z_center + r*cos(θ+Δθ)])
         end
 
         # Create test airfoil data file
@@ -108,7 +114,7 @@ using Serialization
         serialize(polar_path, (alphas, d_trailing_edge_angles, cl_matrix, cd_matrix, cm_matrix))
 
         # Create and serialize obj file
-        faces = [[i, i+1, i+2] for i in 1:2:length(vertices)-2]
+        faces = [[i, i+1, i+2] for i in 1:3:length(vertices)-2]
         open(test_obj_path, "w") do io
             for v in vertices
                 println(io, "v $(v[1]) $(v[2]) $(v[3])")