Add the turnrate vector omega (#65)

* add omega

* use set_va

* fix tests

* remove second parameter yaw rate comment

* rename Uinf to v_a

Author: 1-Bart-1

README.md

@@ -114,7 +114,7 @@ wa = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, (vel_app, 0.0))  # Second parameter is yaw rate
+set_va!(wa, vel_app)
 ```
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`.
 

docs/src/index.md

@@ -115,7 +115,7 @@ wa = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, (vel_app, 0.0))  # Second parameter is yaw rate
+set_va!(wa, vel_app)
 ```
 It is possible to import the wing geometry using an `.obj` file as shown in the example `ram_air_kite.jl`.
 

examples/bench.jl

@@ -31,7 +31,7 @@ wa = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, (vel_app, 0.0))  # Second parameter is yaw rate
+set_va!(wa, vel_app)
 
 # Step 4: Initialize solvers for both LLT and VSM methods
 llt_solver = Solver(aerodynamic_model_type=:LLT)

examples/ram_air_kite.jl

@@ -36,7 +36,7 @@ vel_app = [
     sin(side_slip),
     sin(aoa_rad)
 ] * v_a
-body_aero.va = vel_app
+set_va!(body_aero, vel_app)
 
 # Plotting geometry
 plot && plot_geometry(

examples/rectangular_wing.jl

@@ -31,7 +31,7 @@ wa = BodyAerodynamics([wing])
 
 # Set inflow conditions
 vel_app = [cos(alpha), 0.0, sin(alpha)] .* v_a
-set_va!(wa, (vel_app, 0.0))  # Second parameter is yaw rate
+set_va!(wa, vel_app, [0, 0, 0.1])
 
 # Step 4: Initialize solvers for both LLT and VSM methods
 llt_solver = Solver(aerodynamic_model_type=:LLT)

examples/stall_model.jl

@@ -63,7 +63,7 @@ vel_app = [
     sin(side_slip),
     sin(aoa_rad)
 ] * v_a
-body_aero_CAD_19ribs.va = vel_app
+set_va!(body_aero_CAD_19ribs, vel_app)
 
 # Plotting geometry
 plot && plot_geometry(

src/body_aerodynamics.jl

@@ -6,8 +6,8 @@ Main structure for calculating aerodynamic properties of bodies.
 # Fields
 - panels::Vector{Panel}: Vector of Panel structs
 - wings::Vector{AbstractWing}: a vector of wings; a body can have multiple wings
-- `_va`::Union{Nothing, Vector{Float64}, Tuple{MVec3, Float64}}: A vector of the apparent wind speed,  
-                                                      or a tuple of the v_a vector and yaw rate (rad/s).
+- `va`::Union{Nothing, MVec3}: A vector of the apparent wind speed
+- `omega`::Union{Nothing, MVec3}: A vector of the turn rate around the kite body axes
 - `gamma_distribution`::Union{Nothing, Vector{Float64}}: unclear, please defined
 - `alpha_uncorrected`::Union{Nothing, Vector{Float64}}: unclear, please define
 - `alpha_corrected`::Union{Nothing, Vector{Float64}}: unclear, please define
@@ -16,7 +16,8 @@ Main structure for calculating aerodynamic properties of bodies.
 mutable struct BodyAerodynamics
     panels::Vector{Panel}
     wings::Vector{AbstractWing} # can be a vector of Wings, or of KiteWings
-    _va::Union{Nothing, Vector{Float64}, Tuple{MVec3, Float64}}
+    _va::Union{Nothing, Vector{Float64}}
+    omega::Union{Nothing, MVec3}
     gamma_distribution::Union{Nothing, Vector{Float64}}
     alpha_uncorrected::Union{Nothing, Vector{Float64}}
     alpha_corrected::Union{Nothing, Vector{Float64}}
@@ -75,6 +76,7 @@ mutable struct BodyAerodynamics
             panels,
             wings,
             nothing,  # va
+            nothing,  # omega
             nothing,  # gamma_distribution
             nothing,  # alpha_uncorrected
             nothing,  # alpha_corrected
@@ -95,7 +97,7 @@ function Base.getproperty(obj::BodyAerodynamics, sym::Symbol)
 end
 
 function Base.setproperty!(obj::BodyAerodynamics, sym::Symbol, val)
-    if sym === :va
+    if sym === :va || sym === :omega
         set_va!(obj, val)
     else
         setfield!(obj, sym, val)
@@ -639,49 +641,43 @@ end
 
 
 """
-    set_va!(body_aero::BodyAerodynamics, va::Union{Vector{Float64}, Tuple{Vector{Float64}, Float64}})
+    set_va!(body_aero::BodyAerodynamics, va, omega=zeros(3))
 
 Set velocity array and update wake filaments.
 
 # Arguments
-- `va`: Either a velocity vector or tuple of (velocity vector, yaw_rate)
+- `va`: Either a velocity vector or a list of velocity vectors
+- `omega`: Turn rate around x y and z axis
 """
-function set_va!(body_aero::BodyAerodynamics, va)
-    # Add length check for va_vec
+function set_va!(body_aero::BodyAerodynamics, va, omega=zeros(MVec3))
+    omega = MVec3(omega)
+
+    # Add length check for va
     if va isa Vector{Float64} && length(va) != 3 && length(va) != length(body_aero.panels)
         throw(ArgumentError("va must be length 3 or match number of panels"))
     end
-    # Handle input types
-    va_vec, yaw_rate = if va isa Tuple && length(va) == 2
-        va
-    else
-        (va, 0.0)
-    end
-    
-    # Validate input
-    va_vec = convert(Vector{Float64}, va_vec)
 
     # Calculate va_distribution based on input type
-    va_distribution = if length(va_vec) == 3 && yaw_rate == 0.0
-        repeat(reshape(va_vec, 1, 3), length(body_aero.panels))
-    elseif length(va_vec) == length(body_aero.panels)
-        va_vec
-    elseif yaw_rate != 0.0 && length(va_vec) == 3
-        va_dist = Vector{Float64}[]
+    va_distribution = if length(va) == 3 && all(omega .== 0.0)
+        repeat(reshape(va, 1, 3), length(body_aero.panels))
+    elseif length(va) == length(body_aero.panels)
+        va
+    elseif !all(omega .== 0.0) && length(va) == 3
+        va_dist = zeros(length(body_aero.panels), 3)
 
         for wing in body_aero.wings
             # Get spanwise positions
-            spanwise_positions = [panel.control_point[2] for panel in body_aero.panels]
+            spanwise_positions = [panel.control_point for panel in body_aero.panels]
 
             # Calculate velocities for each panel
             for i in 1:wing.n_panels
-                yaw_rate_apparent_velocity = [-yaw_rate * spanwise_positions[i], 0.0, 0.0]
-                push!(va_dist, yaw_rate_apparent_velocity + va_vec)
+                omega_va = -omega × spanwise_positions[i]
+                va_dist[i, :] .= omega_va .+ va
             end
         end
-        reduce(vcat, va_dist)
+        va_dist
     else
-        throw(ArgumentError("Invalid va distribution: length(va)=$(length(va_vec)) ≠ n_panels=$(length(body_aero.panels))"))
+        throw(ArgumentError("Invalid va distribution: length(va)=$(length(va)) ≠ n_panels=$(length(body_aero.panels))"))
     end
 
     # Update panel velocities
@@ -691,5 +687,6 @@ function set_va!(body_aero::BodyAerodynamics, va)
 
     # Update wake elements
     body_aero.panels = frozen_wake(va_distribution, body_aero.panels)
-    body_aero._va = va_vec
+    body_aero._va = va
+    return nothing
 end

src/filament.jl

@@ -86,14 +86,14 @@ end
     velocity_3D_trailing_vortex(filament::BoundFilament, 
                               XVP::Vector{Float64}, 
                               gamma::Float64, 
-                              Uinf::Float64)
+                              v_a::Float64)
 
 Calculate induced velocity by a trailing vortex filament.
 
 # Arguments
 - `XVP`: Control point coordinates
 - `gamma`: Vortex strength
-- `Uinf`: Inflow velocity magnitude
+- `v_a`: Inflow velocity magnitude
 
 Reference: Rick Damiani et al. "A vortex step method for nonlinear airfoil polar data 
 as implemented in KiteAeroDyn".
@@ -103,7 +103,7 @@ as implemented in KiteAeroDyn".
     filament::BoundFilament,
     XVP,
     gamma,
-    Uinf,
+    v_a,
     work_vectors
 )
     r0, r1, r2, r_perp, r1Xr2, r1Xr0, r2Xr0, normr1r2 = work_vectors[1:8]
@@ -115,7 +115,7 @@ as implemented in KiteAeroDyn".
     r_perp .= dot(r1, r0) .* r0 ./ (norm(r0)^2)
 
     # Cut-off radius
-    epsilon = sqrt(4 * ALPHA0 * NU * norm(r_perp) / Uinf)
+    epsilon = sqrt(4 * ALPHA0 * NU * norm(r_perp) / v_a)
 
     cross3!(r1Xr2, r1, r2)
     cross3!(r1Xr0, r1, r0)
@@ -161,7 +161,7 @@ end
 """
     velocity_3D_trailing_vortex_semiinfinite(filament::SemiInfiniteFilament, 
                                            Vf::Vector{Float64}, XVP::Vector{Float64},
-                                           GAMMA::Float64, Uinf::Float64)
+                                           GAMMA::Float64, v_a::Float64)
 
 Calculate induced velocity by a semi-infinite trailing vortex filament.
 """
@@ -171,7 +171,7 @@ function velocity_3D_trailing_vortex_semiinfinite!(
     Vf,
     XVP,
     GAMMA,
-    Uinf,
+    v_a,
     work_vectors
 )
     r1, r_perp, r1XVf = work_vectors[1:3]
@@ -180,7 +180,7 @@ function velocity_3D_trailing_vortex_semiinfinite!(
 
     # Calculate core radius
     r_perp .= dot(r1, Vf) .* Vf
-    epsilon = sqrt(4 * ALPHA0 * NU * norm(r_perp) / Uinf)
+    epsilon = sqrt(4 * ALPHA0 * NU * norm(r_perp) / v_a)
 
     cross3!(r1XVf, r1, Vf)
 

src/plotting.jl

@@ -639,7 +639,6 @@ function plot_polars(
             # Read all data first
             data = readdlm(path, ',')
             # Skip the header row by taking data from row 2 onwards
-            @show data[1, :]
             data = data[2:end, :]
             push!(polar_data_list, [data[:, 3], data[:, 1], data[:, 2], zeros(length(data[:, 1]))])
         end

src/solver.jl

@@ -221,7 +221,7 @@ function gamma_loop(
     induced_velocity_all = zeros(n_panels, 3)
     relative_velocity_array = similar(va_array)
     relative_velocity_crossz = similar(relative_velocity_array)
-    Uinfcrossz_array = similar(va_array)
+    v_acrossz_array = similar(va_array)
     cl_array = zeros(n_panels)
     damp = zeros(length(gamma))
     v_normal_array = zeros(n_panels)
@@ -249,7 +249,7 @@ function gamma_loop(
                 view(relative_velocity_array, i, :),
                 view(z_airf_array, i, :)
             )
-            Uinfcrossz_array[i, :] .= cross3(
+            v_acrossz_array[i, :] .= cross3(
                 view(va_array, i, :),
                 view(z_airf_array, i, :)
             )
@@ -263,7 +263,7 @@ function gamma_loop(
 
         for i in 1:n_panels
             @views v_a_array[i] = norm(relative_velocity_crossz[i, :])
-            @views Umagw_array[i] = norm(Uinfcrossz_array[i, :])
+            @views Umagw_array[i] = norm(v_acrossz_array[i, :])
         end
 
         for (i, (panel, alpha)) in enumerate(zip(panels, alpha_array))