From 0704473e86285e12d6eaaaa32872ae88c699a8ac Mon Sep 17 00:00:00 2001 From: Ryan Anderson Date: Tue, 12 Mar 2024 12:43:22 -0600 Subject: [PATCH 1/3] fix repeated points check (continue keyword) --- src/wake.jl | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/wake.jl b/src/wake.jl index 153770a..27ff533 100644 --- a/src/wake.jl +++ b/src/wake.jl @@ -147,15 +147,19 @@ end # velocity at the wake shedding locations for is = 1:ns+1 + # skip loop if true + continue_loop = false + # check if this point is a duplicate, skip if it is if (isurf, is) in keys(repeated_points) for (jsurf, js) in repeated_points[(isurf, is)] # NOTE: we assume that a point is not repeated on the same surface if jsurf < isurf wake_velocities[isurf][1, is] = wake_velocities[jsurf][1, js] - continue + continue_loop = true end end + continue_loop && continue end # get vertex location From 57048a2ba3d46361682aa613fd2842ca8d382a40 Mon Sep 17 00:00:00 2001 From: Ryan Anderson Date: Sat, 16 Mar 2024 14:59:54 -0600 Subject: [PATCH 2/3] benchmark copying wake panels --- Project.toml | 1 + src/VortexLattice.jl | 3 +- src/analyses.jl | 5 +- src/fmm.jl | 138 +++ src/induced.jl | 35 + src/panel.jl | 40 +- src/system.jl | 16 +- test/runtests.jl | 2784 +++++++++++++++++++++--------------------- 8 files changed, 1651 insertions(+), 1371 deletions(-) create mode 100644 src/fmm.jl diff --git a/Project.toml b/Project.toml index eb8f62c..3a12923 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ authors = ["Taylor McDonnell "] version = "0.1.3" [deps] +FLOWFMM = "ce07d0d3-2b9f-49ba-89eb-12c800257c85" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/src/VortexLattice.jl b/src/VortexLattice.jl index e703030..18b4076 100644 --- a/src/VortexLattice.jl +++ b/src/VortexLattice.jl @@ -1,8 +1,9 @@ module VortexLattice +using FLOWFMM +using Interpolations using LinearAlgebra using StaticArrays -using Interpolations using WriteVTK # value for dimensionalizing, included just for clarity in the algorithms diff --git a/src/analyses.jl b/src/analyses.jl index 232f267..cf8b490 100644 --- a/src/analyses.jl +++ b/src/analyses.jl @@ -565,6 +565,7 @@ function propagate_system!(system, surfaces, fs, dt; system.nwake .= nwake # update the wake shedding location for this time step + # (based on freestream/kinematic/other velocity only) update_wake_shedding_locations!(wakes, wake_shedding_locations, current_surfaces, ref, fs, dt, additional_velocity, Vte, nwake, eta) @@ -585,7 +586,7 @@ function propagate_system!(system, surfaces, fs, dt; wake_shedding_locations = wake_shedding_locations, trailing_vortices = trailing_vortices) - # calculate RHS + # calculate RHS <-- add FMM here: wake-on-all if derivatives normal_velocity_derivatives!(w, dw, current_surfaces, wakes, ref, fs; additional_velocity, Vcp, symmetric, nwake, @@ -611,6 +612,7 @@ function propagate_system!(system, surfaces, fs, dt; dΓdt ./= dt # divide by corresponding time step # compute transient forces on each panel (if necessary) + # also compute surface-on-all induced velocity if near_field_analysis if derivatives near_field_forces_derivatives!(properties, dproperties, @@ -633,6 +635,7 @@ function propagate_system!(system, surfaces, fs, dt; system.derivatives[] = derivatives # update wake velocities + # ---should already be done if I've done this right--- get_wake_velocities!(wake_velocities, current_surfaces, wakes, ref, fs, Γ, additional_velocity, Vte, symmetric, repeated_points, nwake, surface_id, wake_finite_core, diff --git a/src/fmm.jl b/src/fmm.jl new file mode 100644 index 0000000..5311204 --- /dev/null +++ b/src/fmm.jl @@ -0,0 +1,138 @@ +##### +##### overload fmm compatibility functions for `system::System.fmm_panels` +##### + +Base.getindex(system::System, i, ::FLOWFMM.Position) = system.fmm_panels[i].rcp +Base.getindex(system::System, i, ::FLOWFMM.Radius) = system.fmm_panels[i].radius +# Base.getindex(system::System, i, ::FLOWFMM.VectorPotential) = view(g.potential,2:4,i) +# Base.getindex(system::System, i, ::FLOWFMM.ScalarPotential) = g.potential[1,i] +Base.getindex(system::System, i, ::FLOWFMM.Velocity) = system.fmm_Vcp[i] +# Base.getindex(system::System, i, ::FLOWFMM.VelocityGradient) = reshape(view(g.potential,i_VELOCITY_GRADIENT,i),3,3) +Base.getindex(system::System, i, ::FLOWFMM.ScalarStrength) = system.fmm_panels[i].gamma +Base.getindex(system::System, i, ::FLOWFMM.Body) = system.fmm_panels[i], system.fmm_Vcp[i] +function Base.getindex(system::System, i, ::FLOWFMM.Vertex, i_vertex) + if i_vertex == 1 + return system.fmm_panels[i].rtl + elseif i_vertex == 2 + return system.fmm_panels[i].rbl + elseif i_vertex == 3 + return system.fmm_panels[i].rbr + else # i_vertex == 4 + return system.fmm_panels[i].rtr + end +end +Base.getindex(system::System, i, ::FLOWFMM.Normal) = system.fmm_panels[i].ncp +function Base.setindex!(system::System, val, i, ::FLOWFMM.Body) + panel, velocity = val + system.fmm_panels[i] = panel + system.fmm_Vcp[i] = velocity + return nothing +end +# function Base.setindex!(system::System, val, i, ::FLOWFMM.ScalarPotential) +# g.potential[i_POTENTIAL[1],i] = val +# end +# function Base.setindex!(system::System, val, i, ::FLOWFMM.VectorPotential) +# g.potential[i_POTENTIAL[2:4],i] .= val +# end +function Base.setindex!(system::System, val, i, ::FLOWFMM.Velocity) + system.fmm_Vcp[i] = val +end +# function Base.setindex!(system::System, val, i, ::FLOWFMM.VelocityGradient) +# reshape(g.potential[i_VELOCITY_GRADIENT,i],3,3) .= val +# end +FLOWFMM.get_n_bodies(system::System) = length(system.fmm_panels) +Base.eltype(::System{TF}) where TF = TF + +FLOWFMM.buffer_element(system::System) = deepcopy(system.fmm_panels[1]), deepcopy(system.fmm_Vcp[1]) + +FLOWFMM.B2M!(system::System, branch, bodies_index, harmonics, expansion_order) = + FLOWFMM.B2M!_quadpanel(system, branch, bodies_index, harmonics, system.fmm_p, FLOWFMM.UniformNormalDipolePanel()) + +# no velocity gradient needed if the target is a `::VortexLattice.System` or a FLOWFMM.ProbeSystem without velocity gradient +function FLOWFMM.direct!(target_system::Union{System,FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,Vector{SVector{3,<:Any}},Nothing}}, target_index, source_system::System{TF}, source_index) where TF + # loop over targets + for j_target in target_index + rcp = target_system[j_target,FLOWFMM.POSITION] + V = zero(SVector{3,TF}) + + # loop over source panels + for i_source in source_index + panel = source_system.fmm_panels[i_source] + + # get panel vertices, strength, and core size + r11 = panel.rtl + r12 = panel.rtr + r21 = panel.rbl + r22 = panel.rbr + gamma = panel.gamma + core_size = panel.core_size + + # move origin to control point + r1 = rcp - r11 + r2 = rcp - r12 + r3 = rcp - r22 + r4 = rcp - r21 + + # add induced velocity, ignoring finite core size for now + finite_core = false + this_V = zero(SVector{3,TF}) + this_V += bound_induced_velocity(r1, r2, finite_core, core_size) + this_V += bound_induced_velocity(r2, r3, finite_core, core_size) + this_V += bound_induced_velocity(r3, r4, finite_core, core_size) + this_V += bound_induced_velocity(r4, r1, finite_core, core_size) + V += this_V * gamma + end + + # add to target + target_system[j_target,FLOWFMM.VELOCITY] = target_system[j_target,FLOWFMM.VELOCITY] + V + end +end + +function FLOWFMM.direct!(target_system, target_index, source_system::System{TF}, source_index) where TF + # loop over targets + for j_target in target_index + rcp = target_system[j_target,FLOWFMM.POSITION] + V = zero(SVector{3,TF}) + Vgrad = zero(SMatrix{3,3,TF,9}) + + # loop over source panels + for i_source in source_index + panel = source_system.fmm_panels[i_source] + + # get panel vertices, strength, and core size + r11 = panel.rtl + r12 = panel.rtr + r21 = panel.rbl + r22 = panel.rbr + gamma = panel.gamma + core_size = panel.core_size + + # move origin to control point + r1 = rcp - r11 + r2 = rcp - r12 + r3 = rcp - r22 + r4 = rcp - r21 + + # add induced velocity, ignoring finite core size for now + finite_core = false + this_V = zero(SVector{3,TF}) + this_V += bound_induced_velocity(r1, r2, finite_core, core_size) + this_V += bound_induced_velocity(r2, r3, finite_core, core_size) + this_V += bound_induced_velocity(r3, r4, finite_core, core_size) + this_V += bound_induced_velocity(r4, r1, finite_core, core_size) + V += this_V * gamma + + # velocity gradient + this_Vgrad = zero(SMatrix{3,3,TF,9}) + this_Vgrad += bound_velocity_gradient(r1, r2) + this_Vgrad += bound_velocity_gradient(r2, r3) + this_Vgrad += bound_velocity_gradient(r3, r4) + this_Vgrad += bound_velocity_gradient(r4, r1) + Vgrad += this_Vgrad + end + + # add to target + target_system[j_target,FLOWFMM.VELOCITY] = target_system[j_target,FLOWFMM.VELOCITY] + V + target_system[j_target,FLOWFMM.VELOCITY_GRADIENT] = target_system[j_target,FLOWFMM.VELOCITY_GRADIENT] + Vgrad + end +end \ No newline at end of file diff --git a/src/induced.jl b/src/induced.jl index fd6ae44..f936e41 100644 --- a/src/induced.jl +++ b/src/induced.jl @@ -27,6 +27,41 @@ relative to the end of the bound vortex return Vhat end +function bound_velocity_gradient(r1::SVector{3,TF},r2) where TF + # zeta + r1norm = norm(r1) + r2norm = norm(r2) + dotprod = dot(r1,r2) + t1 = 1/(r1norm*r2norm + dotprod) + t2 = 1/r1norm + 1/r2norm + t3 = 1/4/pi + z = t1*t2*t3 + + # zeta gradient + r1norm3 = r1norm^3 + r2norm3 = r2norm^3 + t4 = SVector{3,TF}(r1[i]/r1norm^3 + r2[i]/r2norm^3 for i in 1:3) + t5 = SVector{3,TF}(r1norm/r2norm*r2[i] + r2norm/r1norm*r1[i] + r1[i] + r2[i] for i in 1:3) + zgrad = t3*(-t1*t4 - t2*t5*t1^2) + + # Omega + o = cross(r1,r2) + + # Omega gradient + ograd = SMatrix{3,3,TF,9}( + 0.0,# 1,1 + r1[3]-r2[3], # 2,1 + r2[2]-r1[2], # 3,1 + r2[3]-r1[3], # 1,2 + 0.0, # 2,2 + r1[1]-r2[1], # 3,2 + r1[2]-r2[2], # 1,3 + r2[1]-r1[1], # 2,3 + 0.0 # 3,3 + ) + return transpose(zgrad * transpose(o)) + z * ograd +end + """ trailing_induced_velocity(r1, r2, xhat, finite_core, core_size) diff --git a/src/panel.jl b/src/panel.jl index f5514ee..34ef0a0 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -390,7 +390,7 @@ end @inline Base.eltype(::Type{WakePanel{TF}}) where TF = TF @inline Base.eltype(::WakePanel{TF}) where TF = TF -WakePanel{TF}(p::WakePanel) where TF = WakePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.core_size, p.chord) +WakePanel{TF}(p::WakePanel) where TF = WakePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.core_size, p.gamma) Base.convert(::Type{WakePanel{TF}}, p::WakePanel) where {TF} = WakePanel{TF}(p) @inline top_left(panel::WakePanel) = panel.rtl @@ -614,3 +614,41 @@ Return induced drag from vortex `j` induced on panel `i` return Di end + +""" + FastMultipolePanel{TF} + +SurfacePanel used for fast multipole acceleration. + +**Fields** + - `rtl`: position of the left side of the top bound vortex + - `rtr`: position of the right side of the top bound vortex + - `rbl`: position of the left side of the bottom bound vortex + - `rbr`: position of the right side of the bottom bound vortex + - `rcp`: position of the panel control point + - `ncp`: normal vector at the panel control point + - `radius`: effective panel radius for determining multipole acceptance + - `core_size`: finite core size (for use when the finite core smoothing model is enabled) + - `gamma`: circulation strength of the panel +""" +struct FastMultipolePanel{TF} + rtl::SVector{3,TF} + rtr::SVector{3,TF} + rbl::SVector{3,TF} + rbr::SVector{3,TF} + ncp::SVector{3,TF} + radius::TF + core_size::TF + gamma::TF +end + +FastMultipolePanel(p::SurfacePanel{TF}, gamma) where TF = FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.ncp, max(norm(p.rtl-p.rbr)/2, norm(p.rbl-p.rtr)/2), p.core_size, gamma) + +function FastMultipolePanel(p::WakePanel{TF}) where TF + v1 = p.rtr-p.rbl + v2 = p.rtl-p.rbr + ncp = cross(v1, v2) + ncp /= norm(ncp) + radius = max(norm(v1)/2, norm(v2)/2) + return FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, ncp, radius, p.core_size, p.gamma) +end \ No newline at end of file diff --git a/src/system.jl b/src/system.jl index ad06a3f..1b1108e 100644 --- a/src/system.jl +++ b/src/system.jl @@ -45,7 +45,7 @@ Contains pre-allocated storage for internal system variables. - `AIC`: Aerodynamic influence coefficient matrix from the surface panels - `w`: Normal velocity at the control points from external sources and wakes - `Γ`: Circulation strength of the surface panels - - `V`: Velocity at the wake vertices for each surface + - `V`: Velocity at the wake vertices for each surface (includes velocity at wake-shed locations) - `surfaces`: Surfaces, represented by matrices of surface panels - `properties`: Surface panel properties for each surface - `wakes`: Wake panel properties for each surface @@ -73,6 +73,13 @@ Contains pre-allocated storage for internal system variables. - `Vv`: Velocity due to surface motion at the vertical bound vortex centers - `Vte`: Velocity due to surface motion at the trailing edge vertices - `dΓdt`: Derivative of the circulation strength with respect to non-dimensional time + - `fmm_panels`: Fast-access copies of all surface panels (including surface-wake transition panels) for fast multipole acceleration. + - `fmm_Vcp`: Velocity at the control points of `fmm_panels` + - `fmm_force_probes`: Locations and induced velocity of the center of the top vortices + - `fmm_wake_probes`: Locations, induced velocity, and induced velocity gradient of the wake particle shed locations + - `fmm_p`: Multipole expansion order + - `fmm_ncrit`: Maximum number of panels in the leaf level of the FMM + - `fmm_theta`: Multipole acceptance criterion for the FMM, between 0 and 1 """ struct System{TF} AIC::Matrix{TF} @@ -103,6 +110,13 @@ struct System{TF} Vv::Vector{Matrix{SVector{3, TF}}} Vte::Vector{Vector{SVector{3, TF}}} dΓdt::Vector{TF} + # fmm_panels::Vector{FastMultipolePanel{TF}} + # fmm_Vcp::Vector{SVector{3,TF}} + # fmm_velocity_probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,Vector{SVector{3,TF}},Nothing} + # fmm_gradient_probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,Vector{SVector{3,TF}},Vector{SMatrix{3,3,TF,9}}} + # fmm_p::Val{Int} + # fmm_ncrit::Int + # fmm_theta::Float64 end @inline Base.eltype(::Type{System{TF}}) where TF = TF diff --git a/test/runtests.jl b/test/runtests.jl index f1b3212..196ca81 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,6 +1,7 @@ -using Test -using VortexLattice +using ForwardDiff using LinearAlgebra +using VortexLattice +using Test ztol = sqrt(eps()) @@ -17,354 +18,1418 @@ function avl_normal_vector(dr, theta) return ncp / norm(ncp) # normal vector used by AVL end -@testset "AVL - Run 1 - Wing with Uniform Spacing" begin - - # Simple Wing with Uniform Spacing - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) - - # vortex rings with symmetry - mirror = false - symmetric = true - - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.24324, atol=1e-3) - @test isapprox(CD, 0.00243, atol=1e-5) - @test isapprox(CDiff, 0.00245, atol=1e-5) - @test isapprox(Cm, -0.02252, atol=1e-4) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) - - # vortex rings with mirrored geometry - mirror = true - symmetric = false - - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.24324, atol=1e-3) - @test isapprox(CD, 0.00243, atol=1e-5) - @test isapprox(CDiff, 0.00245, atol=1e-5) - @test isapprox(Cm, -0.02252, atol=1e-4) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end - -@testset "AVL - Run 2 - Wing with Cosine Spacing" begin - - # Run 2: Simple Wing with Cosine Spacing - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Cosine() - spacing_c = Uniform() - mirror = true - symmetric = false - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) - - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) +# @testset "AVL - Run 1 - Wing with Uniform Spacing" begin - CD, CY, CL = CF - Cl, Cm, Cn = CM +# # Simple Wing with Uniform Spacing - @test isapprox(CL, 0.23744, atol=1e-3) - @test isapprox(CD, 0.00254, atol=1e-5) - @test isapprox(CDiff, 0.00243, atol=1e-5) - @test isapprox(Cm, -0.02165, atol=1e-4) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end - -@testset "AVL - Run 3 - Wing at High Angle of Attack" begin +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() - # Simple Wing at High Angle of Attack - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - symmetric = true +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) - alpha = 8.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) +# # vortex rings with symmetry +# mirror = false +# symmetric = true - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - surfaces = [surface] +# surfaces = [surface] - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - CF, CM = body_forces(system; frame=Stability()) +# CF, CM = body_forces(system; frame=Stability()) - CDiff = far_field_drag(system) +# CDiff = far_field_drag(system) - CD, CY, CL = CF - Cl, Cm, Cn = CM +# CD, CY, CL = CF +# Cl, Cm, Cn = CM - @test isapprox(CL, 0.80348, atol=1e-3) - @test isapprox(CD, 0.02651, atol=1e-4) - @test isapprox(CDiff, 0.02696, atol=1e-5) - @test isapprox(Cm, -0.07399, atol=1e-3) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end +# @test isapprox(CL, 0.24324, atol=1e-3) +# @test isapprox(CD, 0.00243, atol=1e-5) +# @test isapprox(CDiff, 0.00245, atol=1e-5) +# @test isapprox(Cm, -0.02252, atol=1e-4) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) -@testset "AVL - Run 4 - Wing with Dihedral" begin +# # vortex rings with mirrored geometry +# mirror = true +# symmetric = false + +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.24324, atol=1e-3) +# @test isapprox(CD, 0.00243, atol=1e-5) +# @test isapprox(CDiff, 0.00245, atol=1e-5) +# @test isapprox(Cm, -0.02252, atol=1e-4) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end - # Simple Wing with Dihedral +# @testset "AVL - Run 2 - Wing with Cosine Spacing" begin - # NOTE: There is some interaction between twist, dihedral, and chordwise - # position which causes the normal vectors found by AVL to differ from those - # computed by this package. We therefore manually overwrite the normal - # vectors when this occurs in order to get a better comparison. +# # Run 2: Simple Wing with Cosine Spacing - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 3.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - symmetric = true +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Cosine() +# spacing_c = Uniform() +# mirror = true +# symmetric = false - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.23744, atol=1e-3) +# @test isapprox(CD, 0.00254, atol=1e-5) +# @test isapprox(CDiff, 0.00243, atol=1e-5) +# @test isapprox(Cm, -0.02165, atol=1e-4) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 3 - Wing at High Angle of Attack" begin + +# # Simple Wing at High Angle of Attack + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false +# symmetric = true + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 8.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.80348, atol=1e-3) +# @test isapprox(CD, 0.02651, atol=1e-4) +# @test isapprox(CDiff, 0.02696, atol=1e-5) +# @test isapprox(Cm, -0.07399, atol=1e-3) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 4 - Wing with Dihedral" begin + +# # Simple Wing with Dihedral + +# # NOTE: There is some interaction between twist, dihedral, and chordwise +# # position which causes the normal vectors found by AVL to differ from those +# # computed by this package. We therefore manually overwrite the normal +# # vectors when this occurs in order to get a better comparison. + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 3.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false +# symmetric = true + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) + +# # also get normal vector as AVL defines it +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# for (ip, p) in enumerate(surface) +# # check that our normal vector is approximately the same as AVL's +# @test isapprox(p.ncp, ncp, rtol=0.01) +# # replace our normal vector with AVL's normal vector for this test +# surface[ip] = set_normal(p, ncp) +# end + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.24787, atol=1e-3) +# @test isapprox(CD, 0.00246, atol=1e-4) +# @test isapprox(CDiff, 0.00245, atol=1e-5) +# @test isapprox(Cm, -0.02395, atol=1e-4) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 5 - Wing with Dihedral at Very High Angle of Attack" begin + +# # Simple Wing with Dihedral at Very High Angle of Attack + +# # NOTE: this test case is nonphysical, so it just tests the numerics + +# # NOTE: There is some interaction between twist, dihedral, and chordwise +# # position which causes the normal vectors found by AVL to differ from those +# # computed by this package. We therefore manually overwrite the normal +# # vectors when this occurs in order to get a better comparison. + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 3.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false +# symmetric = true + +# # adjust chord length to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 20.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# # vortex rings, untwisted geometry +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# for (ip, p) in enumerate(surface) +# # check that our normal vector is approximately the same as AVL's +# @test isapprox(p.ncp, ncp, rtol=0.01) +# # replace our normal vector with AVL's normal vector for this test +# surface[ip] = set_normal(p, ncp) +# end + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 1.70982, rtol=0.02) +# @test isapprox(CD, 0.12904, rtol=0.02) +# @test isapprox(CDiff, 0.11502, rtol=0.02) +# @test isapprox(Cm, -0.45606, rtol=0.02) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 6 - Wing and Tail without Finite Core Model" begin + +# # Wing and Tail without Finite Core Model + +# # NOTE: AVL's finite-core model is turned off for these tests + +# # NOTE: There is some interaction between twist, dihedral, and chordwise +# # position which causes the normal vectors found by AVL to differ from those +# # computed by this package. We therefore manually overwrite the normal +# # vectors when this occurs in order to get a better comparison. + +# # wing +# xle = [0.0, 0.2] +# yle = [0.0, 5.0] +# zle = [0.0, 1.0] +# chord = [1.0, 0.6] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false + +# # horizontal stabilizer +# xle_h = [0.0, 0.14] +# yle_h = [0.0, 1.25] +# zle_h = [0.0, 0.0] +# chord_h = [0.7, 0.42] +# theta_h = [0.0, 0.0] +# phi_h = [0.0, 0.0] +# ns_h = 6 +# nc_h = 1 +# spacing_s_h = Uniform() +# spacing_c_h = Uniform() +# mirror_h = false + +# # vertical stabilizer +# xle_v = [0.0, 0.14] +# yle_v = [0.0, 0.0] +# zle_v = [0.0, 1.0] +# chord_v = [0.7, 0.42] +# theta_v = [0.0, 0.0] +# phi_v = [0.0, 0.0] +# ns_v = 5 +# nc_v = 1 +# spacing_s_v = Uniform() +# spacing_c_v = Uniform() +# mirror_v = false + +# # adjust chord lengths to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) +# chord_h = @. chord_h/cos(theta_h) +# chord_v = @. chord_v/cos(theta_v) + +# Sref = 9.0 +# cref = 0.9 +# bref = 10.0 +# rref = [0.5, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 5.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = [true, true, false] + +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# # vortex rings - finite core deactivated +# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# for (ip, p) in enumerate(wing) +# # check that our normal vector is approximately the same as AVL's +# @test isapprox(p.ncp, ncp, rtol=0.01) +# # replace our normal vector with AVL's normal vector for this test +# wing[ip] = set_normal(p, ncp) +# end + +# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; +# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) +# translate!(htail, [4.0, 0.0, 0.0]) + +# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; +# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) +# translate!(vtail, [4.0, 0.0, 0.0]) + +# surfaces = [wing, htail, vtail] +# surface_id = [1, 1, 1] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.60408, atol=1e-3) +# @test isapprox(CD, 0.01058, atol=1e-4) +# @test isapprox(CDiff, 0.010378, atol=1e-4) +# @test isapprox(Cm, -0.02778, atol=2e-3) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 7 - Wing and Tail with Finite Core Model" begin + +# # Wing and Tail with Finite Core Model + +# # NOTE: There is some interaction between twist, dihedral, and chordwise +# # position which causes the normal vectors found by AVL to differ from those +# # computed by this package. We therefore manually overwrite the normal +# # vectors when this occurs in order to get a better comparison. + +# # wing +# xle = [0.0, 0.2] +# yle = [0.0, 5.0] +# zle = [0.0, 1.0] +# chord = [1.0, 0.6] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false + +# # horizontal stabilizer +# xle_h = [0.0, 0.14] +# yle_h = [0.0, 1.25] +# zle_h = [0.0, 0.0] +# chord_h = [0.7, 0.42] +# theta_h = [0.0, 0.0] +# phi_h = [0.0, 0.0] +# ns_h = 6 +# nc_h = 1 +# spacing_s_h = Uniform() +# spacing_c_h = Uniform() +# mirror_h = false + +# # vertical stabilizer +# xle_v = [0.0, 0.14] +# yle_v = [0.0, 0.0] +# zle_v = [0.0, 1.0] +# chord_v = [0.7, 0.42] +# theta_v = [0.0, 0.0] +# phi_v = [0.0, 0.0] +# ns_v = 5 +# nc_v = 1 +# spacing_s_v = Uniform() +# spacing_c_v = Uniform() +# mirror_v = false + +# # adjust chord lengths to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) +# chord_h = @. chord_h/cos(theta_h) +# chord_v = @. chord_v/cos(theta_v) + +# Sref = 9.0 +# cref = 0.9 +# bref = 10.0 +# rref = [0.5, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 5.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = [true, true, false] + +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c, fcore = (c,Δs)->max(c/4, Δs/2)) + +# for (ip, p) in enumerate(wing) +# # check that our normal vector is approximately the same as AVL's +# @test isapprox(p.ncp, ncp, rtol=0.01) +# # replace our normal vector with AVL's normal vector for this test +# wing[ip] = set_normal(p, ncp) +# end + +# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; +# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) +# translate!(htail, [4.0, 0.0, 0.0]) + +# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; +# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) +# translate!(vtail, [4.0, 0.0, 0.0]) + +# surfaces = [wing, htail, vtail] +# surface_id = [1, 2, 3] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, derivatives = false) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.60562, atol=1e-3) +# @test isapprox(CD, 0.01058, atol=1e-4) +# @test isapprox(CDiff, 0.0104855, atol=1e-4) +# @test isapprox(Cm, -0.03377, atol=2e-3) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "AVL - Run 8 - Wing with Chordwise Panels" begin + +# # Simple Wing with Chordwise Panels + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 6 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false + +# # adjust chord length to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = true + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.24454, atol=1e-3) +# @test isapprox(CD, 0.00247, atol=1e-5) +# @test isapprox(CDiff, 0.00248, atol=1e-5) +# @test isapprox(Cm, -0.02091, atol=1e-4) +# @test isapprox(CY, 0.0, atol=1e-16) +# @test isapprox(Cl, 0.0, atol=1e-16) +# @test isapprox(Cn, 0.0, atol=1e-16) +# end - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) +# @testset "AVL - Run 9 - Wing with Cosine-Spaced Spanwise and Chordwise Panels" begin - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) +# # Simple Wing with Cosine-Spaced Spanwise and Chordwise Panels - # also get normal vector as AVL defines it - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 6 +# spacing_s = Cosine() +# spacing_c = Cosine() +# mirror = false - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) +# # adjust chord length to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) - for (ip, p) in enumerate(surface) - # check that our normal vector is approximately the same as AVL's - @test isapprox(p.ncp, ncp, rtol=0.01) - # replace our normal vector with AVL's normal vector for this test - surface[ip] = set_normal(p, ncp) - end +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) - surfaces = [surface] +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = true + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.23879, atol=1e-3) +# @test isapprox(CD, 0.00249, atol=1e-5) +# @test isapprox(CDiff, 0.0024626, atol=1e-5) +# @test isapprox(Cm, -0.01995, atol=1e-4) +# @test isapprox(CY, 0.0, atol=1e-16) +# @test isapprox(Cl, 0.0, atol=1e-16) +# @test isapprox(Cn, 0.0, atol=1e-16) +# end - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) +# @testset "AVL - Run 10 - Wing with Sideslip" begin - CF, CM = body_forces(system; frame=Stability()) +# # Simple Wing with Sideslip + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() + +# # adjust chord length to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 15.0*pi/180 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# # vortex rings +# mirror = true +# symmetric = false + +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.22695, atol=1e-3) +# @test isapprox(CD, 0.00227, atol=1e-5) +# @test isapprox(CDiff, 0.0022852, atol=1e-5) +# @test isapprox(Cm, -0.02101, atol=1e-4) +# @test isapprox(CY, 0.0, atol=1e-5) +# @test isapprox(Cl, -0.00644, atol=1e-4) +# @test isapprox(Cn, 0.00012, atol=2e-4) +# end + +# @testset "AVL - Run 11 - Wing Stability Derivatives" begin + +# # Run 11: Simple Wing Stability Derivatives + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = true +# symmetric = false + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# dCF, dCM = stability_derivatives(system) + +# CDa, CYa, CLa = dCF.alpha +# Cla, Cma, Cna = dCM.alpha +# CDb, CYb, CLb = dCF.beta +# Clb, Cmb, Cnb = dCM.beta +# CDp, CYp, CLp = dCF.p +# Clp, Cmp, Cnp = dCM.p +# CDq, CYq, CLq = dCF.q +# Clq, Cmq, Cnq = dCM.q +# CDr, CYr, CLr = dCF.r +# Clr, Cmr, Cnr = dCM.r + +# @test isapprox(CLa, 4.638088, rtol=0.01) +# @test isapprox(CLb, 0.0, atol=ztol) +# @test isapprox(CYa, 0.0, atol=ztol) +# @test isapprox(CYb, -0.000007, atol=1e-4) +# @test isapprox(Cla, 0.0, atol=ztol) +# @test isapprox(Clb, -0.025749, atol=0.001) +# @test isapprox(Cma, -0.429247, rtol=0.01) +# @test isapprox(Cmb, 0.0, atol=ztol) +# @test isapprox(Cna, 0.0, atol=ztol) +# @test isapprox(Cnb, 0.000466, atol=1e-3) +# @test isapprox(Clp, -0.518725, rtol=0.01) +# @test isapprox(Clq, 0.0, atol=ztol) +# @test isapprox(Clr, 0.064243, rtol=0.01) +# @test isapprox(Cmp, 0.0, atol=ztol) +# @test isapprox(Cmq, -0.517094, rtol=0.01) +# @test isapprox(Cmr, 0.0, atol=ztol) +# @test isapprox(Cnp, -0.019846, rtol=0.01) +# @test isapprox(Cnq, 0.0, atol=ztol) +# @test isapprox(Cnr, -0.000898, rtol=0.01) +# end + +# @testset "AVL - Run 12 - Rotational Velocity" begin + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = true +# symmetric = false + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [2*Vinf*0.05/bref; 0.0; 0.0] # nondimensional pbar = 0.05 +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.24323, atol=1e-3) +# @test isapprox(CD, 0.00069, atol=1e-5) +# @test isapprox(Cm, -0.02251, atol=1e-4) +# @test isapprox(CY, 0.00235, atol=2e-4) +# @test isapprox(Cl, -0.02594, atol=1e-4) +# @test isapprox(Cn, -0.00099, atol=2e-4) + +# end + +# @testset "Lifting Line Coefficients" begin +# # Simple Wing with Uniform Spacing + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) + +# # vortex rings with symmetry +# mirror = false +# symmetric = true + +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# grids = [grid] +# surfaces = [surface] + +# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + +# r_ll, c_ll = lifting_line_geometry(grids) + +# cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability()) + +# cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513, +# 0.2404, 0.2233, 0.1952, 0.1434] +# cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024, +# 0.0025, 0.0026, 0.0026, 0.0022] +# cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - CDiff = far_field_drag(system) +# @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf)) +# @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf)) +# @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf)) +# end + +# @testset "Geometry Generation" begin + +# # Tests of the geometry generation functions - CD, CY, CL = CF - Cl, Cm, Cn = CM +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 6 +# spacing_s = Uniform() +# spacing_c = Uniform() - @test isapprox(CL, 0.24787, atol=1e-3) - @test isapprox(CD, 0.00246, atol=1e-4) - @test isapprox(CDiff, 0.00245, atol=1e-5) - @test isapprox(Cm, -0.02395, atol=1e-4) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end +# # Symmetric Geometry -@testset "AVL - Run 5 - Wing with Dihedral at Very High Angle of Attack" begin +# halfgrid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# spacing_s=spacing_s, spacing_c=spacing_c) - # Simple Wing with Dihedral at Very High Angle of Attack +# halfgrid2, surface2 = grid_to_surface_panels(halfgrid1) - # NOTE: this test case is nonphysical, so it just tests the numerics +# halfgrid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; +# spacing_s=spacing_s, spacing_c=spacing_c) - # NOTE: There is some interaction between twist, dihedral, and chordwise - # position which causes the normal vectors found by AVL to differ from those - # computed by this package. We therefore manually overwrite the normal - # vectors when this occurs in order to get a better comparison. +# surface4 = similar(surface1) +# VortexLattice.update_surface_panels!(surface4, halfgrid1) - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 3.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - symmetric = true +# @test isapprox(halfgrid1, halfgrid2) - # adjust chord length to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) +# end +# end - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) +# @test isapprox(halfgrid1, halfgrid3) - alpha = 20.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) +# end +# end - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) +# end +# end - # vortex rings, untwisted geometry - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) +# # Mirrored Geometry - for (ip, p) in enumerate(surface) - # check that our normal vector is approximately the same as AVL's - @test isapprox(p.ncp, ncp, rtol=0.01) - # replace our normal vector with AVL's normal vector for this test - surface[ip] = set_normal(p, ncp) - end +# mirror = true - surfaces = [surface] +# grid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) +# grid2, surface2 = grid_to_surface_panels(halfgrid1; mirror = mirror) - CF, CM = body_forces(system; frame=Stability()) +# grid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; mirror=mirror, +# spacing_s=spacing_s, spacing_c=spacing_c) - CDiff = far_field_drag(system) +# surface4 = similar(surface1) +# VortexLattice.update_surface_panels!(surface4, grid1) - CD, CY, CL = CF - Cl, Cm, Cn = CM +# @test isapprox(grid1, grid2) - @test isapprox(CL, 1.70982, rtol=0.02) - @test isapprox(CD, 0.12904, rtol=0.02) - @test isapprox(CDiff, 0.11502, rtol=0.02) - @test isapprox(Cm, -0.45606, rtol=0.02) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) +# end +# end -@testset "AVL - Run 6 - Wing and Tail without Finite Core Model" begin +# @test isapprox(grid1, grid3) - # Wing and Tail without Finite Core Model +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) +# end +# end - # NOTE: AVL's finite-core model is turned off for these tests +# for I in CartesianIndices(surface1) +# for field in fieldnames(SurfacePanel) +# @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) +# end +# end +# end + +# @testset "Grid Input" begin + +# # Tests for inputting a grid rather than a surface + +# xle = [0.0, 0.4] +# yle = [0.0, 7.5] +# zle = [0.0, 0.0] +# chord = [2.2, 1.8] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 6 +# spacing_s = Uniform() +# spacing_c = Uniform() + +# Sref = 30.0 +# cref = 2.0 +# bref = 15.0 +# rref = [0.50, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 1.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # adjust chord length so x-chord length matches AVL +# chord = @. chord/cos(theta) + +# mirror = false +# symmetric = true + +# halfgrid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# # steady analysis, single grid + +# grids = [halfgrid] + +# system = steady_analysis(grids, ref, fs; symmetric=symmetric) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# @test isapprox(CL, 0.24454, atol=1e-3) +# @test isapprox(CD, 0.00247, atol=1e-5) +# @test isapprox(CDiff, 0.00248, atol=1e-5) +# @test isapprox(Cm, -0.02091, atol=1e-4) +# @test isapprox(CY, 0.0, atol=1e-16) +# @test isapprox(Cl, 0.0, atol=1e-16) +# @test isapprox(Cn, 0.0, atol=1e-16) + +# # steady analysis multiple grids + +# # NOTE: AVL's finite-core model is turned off for these tests + +# # NOTE: There is some interaction between twist, dihedral, and chordwise +# # position which causes the normal vectors found by AVL to differ from those +# # computed by this package. We therefore manually overwrite the normal +# # vectors when this occurs in order to get a better comparison. + +# # wing +# xle = [0.0, 0.2] +# yle = [0.0, 5.0] +# zle = [0.0, 1.0] +# chord = [1.0, 0.6] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false + +# # horizontal stabilizer +# xle_h = [0.0, 0.14] +# yle_h = [0.0, 1.25] +# zle_h = [0.0, 0.0] +# chord_h = [0.7, 0.42] +# theta_h = [0.0, 0.0] +# phi_h = [0.0, 0.0] +# ns_h = 6 +# nc_h = 1 +# spacing_s_h = Uniform() +# spacing_c_h = Uniform() +# mirror_h = false + +# # vertical stabilizer +# xle_v = [0.0, 0.14] +# yle_v = [0.0, 0.0] +# zle_v = [0.0, 1.0] +# chord_v = [0.7, 0.42] +# theta_v = [0.0, 0.0] +# phi_v = [0.0, 0.0] +# ns_v = 5 +# nc_v = 1 +# spacing_s_v = Uniform() +# spacing_c_v = Uniform() +# mirror_v = false + +# # adjust chord lengths to match AVL (which uses chord length in the x-direction) +# chord = @. chord/cos(theta) +# chord_h = @. chord_h/cos(theta_h) +# chord_v = @. chord_v/cos(theta_v) + +# Sref = 9.0 +# cref = 0.9 +# bref = 10.0 +# rref = [0.5, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 5.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = [true, true, false] + +# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + +# # vortex rings - finite core deactivated +# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; +# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) +# translate!(hgrid, [4.0, 0.0, 0.0]) +# translate!(htail, [4.0, 0.0, 0.0]) + +# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; +# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) +# translate!(vgrid, [4.0, 0.0, 0.0]) +# translate!(vtail, [4.0, 0.0, 0.0]) + +# grids = [wgrid, hgrid, vgrid] +# surface_id = [1, 1, 1] + +# system = steady_analysis(grids, ref, fs; symmetric=symmetric, surface_id=surface_id) + +# CF, CM = body_forces(system; frame=Stability()) + +# CDiff = far_field_drag(system) + +# CD, CY, CL = CF +# Cl, Cm, Cn = CM + +# # some differences are expected since we don't set the normal vector in +# # VortexLattice equal to that used by AVL + +# @test isapprox(CL, 0.60408, rtol=0.01) +# @test isapprox(CD, 0.01058, rtol=0.01) +# @test isapprox(CDiff, 0.010378, atol=0.02) +# @test isapprox(Cm, -0.02778, rtol=0.01) +# @test isapprox(CY, 0.0, atol=ztol) +# @test isapprox(Cl, 0.0, atol=ztol) +# @test isapprox(Cn, 0.0, atol=ztol) +# end + +# @testset "Update Trailing Edge Coefficients" begin + +# # This test checks whether the function which updates the trailing edge +# # coefficients `update_trailing_edge_coefficients!` results in the same +# # trailing edge coefficients as `influnece_coefficients!` + +# # wing +# xle = [0.0, 0.2] +# yle = [0.0, 5.0] +# zle = [0.0, 1.0] +# chord = [1.0, 0.6] +# theta = [2.0*pi/180, 2.0*pi/180] +# phi = [0.0, 0.0] +# ns = 12 +# nc = 1 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false + +# # horizontal stabilizer +# xle_h = [0.0, 0.14] +# yle_h = [0.0, 1.25] +# zle_h = [0.0, 0.0] +# chord_h = [0.7, 0.42] +# theta_h = [0.0, 0.0] +# phi_h = [0.0, 0.0] +# ns_h = 6 +# nc_h = 1 +# spacing_s_h = Uniform() +# spacing_c_h = Uniform() +# mirror_h = false + +# # vertical stabilizer +# xle_v = [0.0, 0.14] +# yle_v = [0.0, 0.0] +# zle_v = [0.0, 1.0] +# chord_v = [0.7, 0.42] +# theta_v = [0.0, 0.0] +# phi_v = [0.0, 0.0] +# ns_v = 5 +# nc_v = 1 +# spacing_s_v = Uniform() +# spacing_c_v = Uniform() +# mirror_v = false + +# Sref = 9.0 +# cref = 0.9 +# bref = 10.0 +# rref = [0.5, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# alpha = 5.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# symmetric = [true, true, false] + +# # horseshoe vortices +# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; +# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) +# translate!(htail, [4.0, 0.0, 0.0]) + +# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; +# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) +# translate!(vtail, [4.0, 0.0, 0.0]) + +# surfaces = [wing, htail, vtail] +# surface_id = [1, 2, 2] + +# # number of panels +# N = nc*ns + nc_h*ns_h + nc_v*ns_v +# AIC1 = zeros(N, N) + +# VortexLattice.influence_coefficients!(AIC1, surfaces; +# symmetric = symmetric, +# trailing_vortices = fill(false, length(surfaces)), +# surface_id = surface_id) + +# AIC2 = copy(AIC1) + +# VortexLattice.update_trailing_edge_coefficients!(AIC2, surfaces; +# symmetric = symmetric, +# trailing_vortices = fill(false, length(surfaces)), +# surface_id = surface_id) + +# @test isapprox(AIC1, AIC2) +# end + +# @testset "Wake Induced Velocity" begin + +# # This test constructs two surfaces, one using surface panels, +# # and one using wake panels. It then tests that the wake panel +# # implementations of `induced_velocity` yields identical results to the +# # surface panel implementations + +# # generate the surface/wake geometry +# nc = 5 +# ns = 10 + +# x = range(0, 1, length = nc+1) +# y = range(0, 2, length = ns+1) +# z = range(0, 3, length = ns+1) + +# surface = Matrix{SurfacePanel{Float64}}(undef, nc, ns) +# wake = Matrix{WakePanel{Float64}}(undef, nc, ns) +# Γ = rand(nc, ns) + +# for i = 1:nc, j = 1:ns +# rtl = [x[i], y[j], z[j]] +# rtr = [x[i], y[j+1], z[j+1]] +# rbl = [x[i+1], y[j], z[j]] +# rbr = [x[i+1], y[j+1], z[j+1]] +# rcp = (rtl + rtr + rbl + rbr)/4 +# ncp = cross(rcp - rtr, rcp - rtl) +# core_size = 0.1 +# chord = 0.0 # only used for unsteady simuulations + +# surface[i,j] = SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord) +# wake[i,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, Γ[i, j]) +# end + +# # Test induced velocity calculation at an arbitrary point in space: +# rcp = [10, 11, 12] + +# # no finite core model, no symmetry, no trailing vortices + +# Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; +# finite_core = false, +# symmetric = false, +# trailing_vortices = false, +# xhat = [1, 0, 0]) + +# Vw = VortexLattice.induced_velocity(rcp, wake; +# finite_core = false, +# symmetric = false, +# trailing_vortices = false, +# xhat = [1, 0, 0]) + +# @test isapprox(Vs, Vw) + +# # finite core model, symmetry, and trailing vortices + +# Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; +# finite_core = true, +# symmetric = true, +# trailing_vortices = true, +# xhat = [1, 0, 0]) + +# Vw = VortexLattice.induced_velocity(rcp, wake; +# finite_core = true, +# symmetric = true, +# trailing_vortices = true, +# xhat = [1, 0, 0]) + +# @test isapprox(Vs, Vw) + +# # Test induced velocity calculation at a trailing edge point: +# is = 3 # trailing edge index +# I = CartesianIndex(nc+1, is) # index on surface + +# # no finite core model, no symmetry, no trailing vortices +# Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; +# finite_core = false, +# symmetric = false, +# trailing_vortices = false, +# xhat = [1, 0, 0]) + +# Vw = VortexLattice.induced_velocity(I, wake; +# finite_core = false, +# symmetric = false, +# trailing_vortices = false, +# xhat = [1, 0, 0]) + +# @test isapprox(Vs, Vw) + +# # finite core model, symmetry, and trailing vortices +# Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; +# finite_core = true, +# symmetric = true, +# trailing_vortices = true, +# xhat = [1, 0, 0]) + +# Vw = VortexLattice.induced_velocity(I, wake; +# finite_core = true, +# symmetric = true, +# trailing_vortices = true, +# xhat = [1, 0, 0]) + +# @test isapprox(Vs, Vw) + +# end + +# @testset "Unsteady Vortex Lattice Method - Rectangular Wing" begin + +# Uinf = 1.0 +# AR = 4 + +# # reference parameters +# cref = 1.0 +# bref = AR +# Sref = bref*cref +# rref = [0.0, 0.0, 0.0] +# Vinf = 1.0 +# ref = Reference(Sref, cref, bref, rref, Vinf) + +# # freestream parameters +# alpha = 5.0*pi/180 +# beta = 0.0 +# Omega = [0.0; 0.0; 0.0] +# fs = Freestream(Vinf, alpha, beta, Omega) + +# # geometry +# xle = [0.0, 0.0] +# yle = [-bref/2, bref/2] +# zle = [0.0, 0.0] +# chord = [cref, cref] +# theta = [0.0, 0.0] +# phi = [0.0, 0.0] +# ns = 13 +# nc = 4 +# spacing_s = Uniform() +# spacing_c = Uniform() +# mirror = false +# symmetric = false + +# # non-dimensional time +# t = range(0.0, 10.0, step=0.2) +# dt = [t[i+1]-t[i] for i = 1:length(t)-1] + +# # create vortex rings +# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; +# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + +# surfaces = [surface] + +# # run analysis +# system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; +# symmetric=symmetric) + +# # extract forces at each time step +# CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) +# end + +# @testset "Unsteady Vortex Lattice Method - Wing + Tail" begin - # NOTE: There is some interaction between twist, dihedral, and chordwise - # position which causes the normal vectors found by AVL to differ from those - # computed by this package. We therefore manually overwrite the normal - # vectors when this occurs in order to get a better comparison. + # Unsteady Wing and Tail + ns_multiplier = 1 + nc_multiplier = 1 + time_step = 6.4 # in [6.4, 3.2, 1.6, 0.8, 0.4, 0.2, 0.1, 0.05] # wing xle = [0.0, 0.2] @@ -373,8 +1438,8 @@ end chord = [1.0, 0.6] theta = [2.0*pi/180, 2.0*pi/180] phi = [0.0, 0.0] - ns = 12 - nc = 1 + ns = 12 * ns_multiplier + nc = 1 * nc_multiplier spacing_s = Uniform() spacing_c = Uniform() mirror = false @@ -386,8 +1451,8 @@ end chord_h = [0.7, 0.42] theta_h = [0.0, 0.0] phi_h = [0.0, 0.0] - ns_h = 6 - nc_h = 1 + ns_h = 6 * ns_multiplier + nc_h = 1 * nc_multiplier spacing_s_h = Uniform() spacing_c_h = Uniform() mirror_h = false @@ -399,8 +1464,8 @@ end chord_v = [0.7, 0.42] theta_v = [0.0, 0.0] phi_v = [0.0, 0.0] - ns_v = 5 - nc_v = 1 + ns_v = 5 * ns_multiplier + nc_v = 1 * nc_multiplier spacing_s_v = Uniform() spacing_c_v = Uniform() mirror_v = false @@ -424,127 +1489,10 @@ end symmetric = [true, true, false] - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - - # vortex rings - finite core deactivated + # horseshoe vortices wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - for (ip, p) in enumerate(wing) - # check that our normal vector is approximately the same as AVL's - @test isapprox(p.ncp, ncp, rtol=0.01) - # replace our normal vector with AVL's normal vector for this test - wing[ip] = set_normal(p, ncp) - end - - hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; - mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) - translate!(htail, [4.0, 0.0, 0.0]) - - vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; - mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) - translate!(vtail, [4.0, 0.0, 0.0]) - - surfaces = [wing, htail, vtail] - surface_id = [1, 1, 1] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.60408, atol=1e-3) - @test isapprox(CD, 0.01058, atol=1e-4) - @test isapprox(CDiff, 0.010378, atol=1e-4) - @test isapprox(Cm, -0.02778, atol=2e-3) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end - -@testset "AVL - Run 7 - Wing and Tail with Finite Core Model" begin - - # Wing and Tail with Finite Core Model - - # NOTE: There is some interaction between twist, dihedral, and chordwise - # position which causes the normal vectors found by AVL to differ from those - # computed by this package. We therefore manually overwrite the normal - # vectors when this occurs in order to get a better comparison. - - # wing - xle = [0.0, 0.2] - yle = [0.0, 5.0] - zle = [0.0, 1.0] - chord = [1.0, 0.6] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - - # horizontal stabilizer - xle_h = [0.0, 0.14] - yle_h = [0.0, 1.25] - zle_h = [0.0, 0.0] - chord_h = [0.7, 0.42] - theta_h = [0.0, 0.0] - phi_h = [0.0, 0.0] - ns_h = 6 - nc_h = 1 - spacing_s_h = Uniform() - spacing_c_h = Uniform() - mirror_h = false - - # vertical stabilizer - xle_v = [0.0, 0.14] - yle_v = [0.0, 0.0] - zle_v = [0.0, 1.0] - chord_v = [0.7, 0.42] - theta_v = [0.0, 0.0] - phi_v = [0.0, 0.0] - ns_v = 5 - nc_v = 1 - spacing_s_v = Uniform() - spacing_c_v = Uniform() - mirror_v = false - - # adjust chord lengths to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) - chord_h = @. chord_h/cos(theta_h) - chord_v = @. chord_v/cos(theta_v) - - Sref = 9.0 - cref = 0.9 - bref = 10.0 - rref = [0.5, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 5.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = [true, true, false] - - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - - wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c, fcore = (c,Δs)->max(c/4, Δs/2)) - - for (ip, p) in enumerate(wing) - # check that our normal vector is approximately the same as AVL's - @test isapprox(p.ncp, ncp, rtol=0.01) - # replace our normal vector with AVL's normal vector for this test - wing[ip] = set_normal(p, ncp) - end - hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) translate!(htail, [4.0, 0.0, 0.0]) @@ -556,956 +1504,58 @@ end surfaces = [wing, htail, vtail] surface_id = [1, 2, 3] - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, derivatives = false) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.60562, atol=1e-3) - @test isapprox(CD, 0.01058, atol=1e-4) - @test isapprox(CDiff, 0.0104855, atol=1e-4) - @test isapprox(Cm, -0.03377, atol=2e-3) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end - -@testset "AVL - Run 8 - Wing with Chordwise Panels" begin - - # Simple Wing with Chordwise Panels - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 6 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - - # adjust chord length to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = true - - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.24454, atol=1e-3) - @test isapprox(CD, 0.00247, atol=1e-5) - @test isapprox(CDiff, 0.00248, atol=1e-5) - @test isapprox(Cm, -0.02091, atol=1e-4) - @test isapprox(CY, 0.0, atol=1e-16) - @test isapprox(Cl, 0.0, atol=1e-16) - @test isapprox(Cn, 0.0, atol=1e-16) -end - -@testset "AVL - Run 9 - Wing with Cosine-Spaced Spanwise and Chordwise Panels" begin - - # Simple Wing with Cosine-Spaced Spanwise and Chordwise Panels - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 6 - spacing_s = Cosine() - spacing_c = Cosine() - mirror = false - - # adjust chord length to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = true - - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM + # t + t = range(0.0, 10.0, step=time_step) + dt = t[2:end] - t[1:end-1] - @test isapprox(CL, 0.23879, atol=1e-3) - @test isapprox(CD, 0.00249, atol=1e-5) - @test isapprox(CDiff, 0.0024626, atol=1e-5) - @test isapprox(Cm, -0.01995, atol=1e-4) - @test isapprox(CY, 0.0, atol=1e-16) - @test isapprox(Cl, 0.0, atol=1e-16) - @test isapprox(Cn, 0.0, atol=1e-16) -end + println("Running VortexLattice (nsteps = $(length(t)))") + @time begin + system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; symmetric) -@testset "AVL - Run 10 - Wing with Sideslip" begin + # extract forces at each time step + CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) + end + # copy to fast-access vector + fast_access_vector = Vector{VortexLattice.FastMultipolePanel{Float64}}(undef,sum(length.(system.wakes))) + + function update_fmm_panels!(fmm_panels, wake_panels::Vector{Matrix{VortexLattice.WakePanel{TF}}}) where TF + i_start = 0 + for wake in wake_panels + for (i,panel) in enumerate(wake) + fmm_panels[i_start+i] = VortexLattice.FastMultipolePanel(panel) + end + i_start += length(wake) + end + end - # Simple Wing with Sideslip + using BenchmarkTools + @show sum(length.(system.wakes)) + @btime update_fmm_panels!($fast_access_vector, $system.wakes) - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() +# end - # adjust chord length to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) +# @testset "velocity gradient" begin - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) +# # test velocity gradient of a vortex filament of unit strength +# function test_velocity_gradient(r1,r2) +# function velocity(x) +# return VortexLattice.bound_induced_velocity(r1+x, r2+x, false, 1e-2) +# end +# return ForwardDiff.jacobian(velocity, zeros(3)) +# end - alpha = 1.0*pi/180 - beta = 15.0*pi/180 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) +# r1 = SVector{3}(0.21490591034021655, +# 0.5811234218620989, +# 0.8079954052156225) - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) +# r2 = SVector{3}(0.04054862116722213, +# 0.6134119927849612, +# 0.6270695997855336) - # vortex rings - mirror = true - symmetric = false +# analytic_gradient = VortexLattice.bound_velocity_gradient(r1,r2) +# check = ForwardDiff.jacobian((x)->VortexLattice.bound_induced_velocity(r1+x,r2+x,false,1e-2),zero(SVector{3,Float64})) +# for i in eachindex(check) +# @show isapprox(analytic_gradient[i], check[i]; atol=1e-12) +# end - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.22695, atol=1e-3) - @test isapprox(CD, 0.00227, atol=1e-5) - @test isapprox(CDiff, 0.0022852, atol=1e-5) - @test isapprox(Cm, -0.02101, atol=1e-4) - @test isapprox(CY, 0.0, atol=1e-5) - @test isapprox(Cl, -0.00644, atol=1e-4) - @test isapprox(Cn, 0.00012, atol=2e-4) -end - -@testset "AVL - Run 11 - Wing Stability Derivatives" begin - - # Run 11: Simple Wing Stability Derivatives - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = true - symmetric = false - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - dCF, dCM = stability_derivatives(system) - - CDa, CYa, CLa = dCF.alpha - Cla, Cma, Cna = dCM.alpha - CDb, CYb, CLb = dCF.beta - Clb, Cmb, Cnb = dCM.beta - CDp, CYp, CLp = dCF.p - Clp, Cmp, Cnp = dCM.p - CDq, CYq, CLq = dCF.q - Clq, Cmq, Cnq = dCM.q - CDr, CYr, CLr = dCF.r - Clr, Cmr, Cnr = dCM.r - - @test isapprox(CLa, 4.638088, rtol=0.01) - @test isapprox(CLb, 0.0, atol=ztol) - @test isapprox(CYa, 0.0, atol=ztol) - @test isapprox(CYb, -0.000007, atol=1e-4) - @test isapprox(Cla, 0.0, atol=ztol) - @test isapprox(Clb, -0.025749, atol=0.001) - @test isapprox(Cma, -0.429247, rtol=0.01) - @test isapprox(Cmb, 0.0, atol=ztol) - @test isapprox(Cna, 0.0, atol=ztol) - @test isapprox(Cnb, 0.000466, atol=1e-3) - @test isapprox(Clp, -0.518725, rtol=0.01) - @test isapprox(Clq, 0.0, atol=ztol) - @test isapprox(Clr, 0.064243, rtol=0.01) - @test isapprox(Cmp, 0.0, atol=ztol) - @test isapprox(Cmq, -0.517094, rtol=0.01) - @test isapprox(Cmr, 0.0, atol=ztol) - @test isapprox(Cnp, -0.019846, rtol=0.01) - @test isapprox(Cnq, 0.0, atol=ztol) - @test isapprox(Cnr, -0.000898, rtol=0.01) -end - -@testset "AVL - Run 12 - Rotational Velocity" begin - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = true - symmetric = false - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [2*Vinf*0.05/bref; 0.0; 0.0] # nondimensional pbar = 0.05 - fs = Freestream(Vinf, alpha, beta, Omega) - - # vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.24323, atol=1e-3) - @test isapprox(CD, 0.00069, atol=1e-5) - @test isapprox(Cm, -0.02251, atol=1e-4) - @test isapprox(CY, 0.00235, atol=2e-4) - @test isapprox(Cl, -0.02594, atol=1e-4) - @test isapprox(Cn, -0.00099, atol=2e-4) - -end - -@testset "Lifting Line Coefficients" begin - # Simple Wing with Uniform Spacing - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) - - # vortex rings with symmetry - mirror = false - symmetric = true - - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - grids = [grid] - surfaces = [surface] - - system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - - r_ll, c_ll = lifting_line_geometry(grids) - - cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability()) - - cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513, - 0.2404, 0.2233, 0.1952, 0.1434] - cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024, - 0.0025, 0.0026, 0.0026, 0.0022] - cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] - - @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf)) - @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf)) - @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf)) -end - -@testset "Geometry Generation" begin - - # Tests of the geometry generation functions - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 6 - spacing_s = Uniform() - spacing_c = Uniform() - - # Symmetric Geometry - - halfgrid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - spacing_s=spacing_s, spacing_c=spacing_c) - - halfgrid2, surface2 = grid_to_surface_panels(halfgrid1) - - halfgrid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; - spacing_s=spacing_s, spacing_c=spacing_c) - - surface4 = similar(surface1) - VortexLattice.update_surface_panels!(surface4, halfgrid1) - - @test isapprox(halfgrid1, halfgrid2) - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) - end - end - - @test isapprox(halfgrid1, halfgrid3) - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) - end - end - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) - end - end - - # Mirrored Geometry - - mirror = true - - grid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - grid2, surface2 = grid_to_surface_panels(halfgrid1; mirror = mirror) - - grid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; mirror=mirror, - spacing_s=spacing_s, spacing_c=spacing_c) - - surface4 = similar(surface1) - VortexLattice.update_surface_panels!(surface4, grid1) - - @test isapprox(grid1, grid2) - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) - end - end - - @test isapprox(grid1, grid3) - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) - end - end - - for I in CartesianIndices(surface1) - for field in fieldnames(SurfacePanel) - @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) - end - end -end - -@testset "Grid Input" begin - - # Tests for inputting a grid rather than a surface - - xle = [0.0, 0.4] - yle = [0.0, 7.5] - zle = [0.0, 0.0] - chord = [2.2, 1.8] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 6 - spacing_s = Uniform() - spacing_c = Uniform() - - Sref = 30.0 - cref = 2.0 - bref = 15.0 - rref = [0.50, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 1.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # adjust chord length so x-chord length matches AVL - chord = @. chord/cos(theta) - - mirror = false - symmetric = true - - halfgrid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - # steady analysis, single grid - - grids = [halfgrid] - - system = steady_analysis(grids, ref, fs; symmetric=symmetric) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - @test isapprox(CL, 0.24454, atol=1e-3) - @test isapprox(CD, 0.00247, atol=1e-5) - @test isapprox(CDiff, 0.00248, atol=1e-5) - @test isapprox(Cm, -0.02091, atol=1e-4) - @test isapprox(CY, 0.0, atol=1e-16) - @test isapprox(Cl, 0.0, atol=1e-16) - @test isapprox(Cn, 0.0, atol=1e-16) - - # steady analysis multiple grids - - # NOTE: AVL's finite-core model is turned off for these tests - - # NOTE: There is some interaction between twist, dihedral, and chordwise - # position which causes the normal vectors found by AVL to differ from those - # computed by this package. We therefore manually overwrite the normal - # vectors when this occurs in order to get a better comparison. - - # wing - xle = [0.0, 0.2] - yle = [0.0, 5.0] - zle = [0.0, 1.0] - chord = [1.0, 0.6] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - - # horizontal stabilizer - xle_h = [0.0, 0.14] - yle_h = [0.0, 1.25] - zle_h = [0.0, 0.0] - chord_h = [0.7, 0.42] - theta_h = [0.0, 0.0] - phi_h = [0.0, 0.0] - ns_h = 6 - nc_h = 1 - spacing_s_h = Uniform() - spacing_c_h = Uniform() - mirror_h = false - - # vertical stabilizer - xle_v = [0.0, 0.14] - yle_v = [0.0, 0.0] - zle_v = [0.0, 1.0] - chord_v = [0.7, 0.42] - theta_v = [0.0, 0.0] - phi_v = [0.0, 0.0] - ns_v = 5 - nc_v = 1 - spacing_s_v = Uniform() - spacing_c_v = Uniform() - mirror_v = false - - # adjust chord lengths to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) - chord_h = @. chord_h/cos(theta_h) - chord_v = @. chord_v/cos(theta_v) - - Sref = 9.0 - cref = 0.9 - bref = 10.0 - rref = [0.5, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 5.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = [true, true, false] - - ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - - # vortex rings - finite core deactivated - wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; - mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) - translate!(hgrid, [4.0, 0.0, 0.0]) - translate!(htail, [4.0, 0.0, 0.0]) - - vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; - mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) - translate!(vgrid, [4.0, 0.0, 0.0]) - translate!(vtail, [4.0, 0.0, 0.0]) - - grids = [wgrid, hgrid, vgrid] - surface_id = [1, 1, 1] - - system = steady_analysis(grids, ref, fs; symmetric=symmetric, surface_id=surface_id) - - CF, CM = body_forces(system; frame=Stability()) - - CDiff = far_field_drag(system) - - CD, CY, CL = CF - Cl, Cm, Cn = CM - - # some differences are expected since we don't set the normal vector in - # VortexLattice equal to that used by AVL - - @test isapprox(CL, 0.60408, rtol=0.01) - @test isapprox(CD, 0.01058, rtol=0.01) - @test isapprox(CDiff, 0.010378, atol=0.02) - @test isapprox(Cm, -0.02778, rtol=0.01) - @test isapprox(CY, 0.0, atol=ztol) - @test isapprox(Cl, 0.0, atol=ztol) - @test isapprox(Cn, 0.0, atol=ztol) -end - -@testset "Update Trailing Edge Coefficients" begin - - # This test checks whether the function which updates the trailing edge - # coefficients `update_trailing_edge_coefficients!` results in the same - # trailing edge coefficients as `influnece_coefficients!` - - # wing - xle = [0.0, 0.2] - yle = [0.0, 5.0] - zle = [0.0, 1.0] - chord = [1.0, 0.6] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - - # horizontal stabilizer - xle_h = [0.0, 0.14] - yle_h = [0.0, 1.25] - zle_h = [0.0, 0.0] - chord_h = [0.7, 0.42] - theta_h = [0.0, 0.0] - phi_h = [0.0, 0.0] - ns_h = 6 - nc_h = 1 - spacing_s_h = Uniform() - spacing_c_h = Uniform() - mirror_h = false - - # vertical stabilizer - xle_v = [0.0, 0.14] - yle_v = [0.0, 0.0] - zle_v = [0.0, 1.0] - chord_v = [0.7, 0.42] - theta_v = [0.0, 0.0] - phi_v = [0.0, 0.0] - ns_v = 5 - nc_v = 1 - spacing_s_v = Uniform() - spacing_c_v = Uniform() - mirror_v = false - - Sref = 9.0 - cref = 0.9 - bref = 10.0 - rref = [0.5, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 5.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = [true, true, false] - - # horseshoe vortices - wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; - mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) - translate!(htail, [4.0, 0.0, 0.0]) - - vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; - mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) - translate!(vtail, [4.0, 0.0, 0.0]) - - surfaces = [wing, htail, vtail] - surface_id = [1, 2, 2] - - # number of panels - N = nc*ns + nc_h*ns_h + nc_v*ns_v - AIC1 = zeros(N, N) - - VortexLattice.influence_coefficients!(AIC1, surfaces; - symmetric = symmetric, - trailing_vortices = fill(false, length(surfaces)), - surface_id = surface_id) - - AIC2 = copy(AIC1) - - VortexLattice.update_trailing_edge_coefficients!(AIC2, surfaces; - symmetric = symmetric, - trailing_vortices = fill(false, length(surfaces)), - surface_id = surface_id) - - @test isapprox(AIC1, AIC2) -end - -@testset "Wake Induced Velocity" begin - - # This test constructs two surfaces, one using surface panels, - # and one using wake panels. It then tests that the wake panel - # implementations of `induced_velocity` yields identical results to the - # surface panel implementations - - # generate the surface/wake geometry - nc = 5 - ns = 10 - - x = range(0, 1, length = nc+1) - y = range(0, 2, length = ns+1) - z = range(0, 3, length = ns+1) - - surface = Matrix{SurfacePanel{Float64}}(undef, nc, ns) - wake = Matrix{WakePanel{Float64}}(undef, nc, ns) - Γ = rand(nc, ns) - - for i = 1:nc, j = 1:ns - rtl = [x[i], y[j], z[j]] - rtr = [x[i], y[j+1], z[j+1]] - rbl = [x[i+1], y[j], z[j]] - rbr = [x[i+1], y[j+1], z[j+1]] - rcp = (rtl + rtr + rbl + rbr)/4 - ncp = cross(rcp - rtr, rcp - rtl) - core_size = 0.1 - chord = 0.0 # only used for unsteady simuulations - - surface[i,j] = SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord) - wake[i,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, Γ[i, j]) - end - - # Test induced velocity calculation at an arbitrary point in space: - rcp = [10, 11, 12] - - # no finite core model, no symmetry, no trailing vortices - - Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; - finite_core = false, - symmetric = false, - trailing_vortices = false, - xhat = [1, 0, 0]) - - Vw = VortexLattice.induced_velocity(rcp, wake; - finite_core = false, - symmetric = false, - trailing_vortices = false, - xhat = [1, 0, 0]) - - @test isapprox(Vs, Vw) - - # finite core model, symmetry, and trailing vortices - - Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; - finite_core = true, - symmetric = true, - trailing_vortices = true, - xhat = [1, 0, 0]) - - Vw = VortexLattice.induced_velocity(rcp, wake; - finite_core = true, - symmetric = true, - trailing_vortices = true, - xhat = [1, 0, 0]) - - @test isapprox(Vs, Vw) - - # Test induced velocity calculation at a trailing edge point: - is = 3 # trailing edge index - I = CartesianIndex(nc+1, is) # index on surface - - # no finite core model, no symmetry, no trailing vortices - Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; - finite_core = false, - symmetric = false, - trailing_vortices = false, - xhat = [1, 0, 0]) - - Vw = VortexLattice.induced_velocity(I, wake; - finite_core = false, - symmetric = false, - trailing_vortices = false, - xhat = [1, 0, 0]) - - @test isapprox(Vs, Vw) - - # finite core model, symmetry, and trailing vortices - Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; - finite_core = true, - symmetric = true, - trailing_vortices = true, - xhat = [1, 0, 0]) - - Vw = VortexLattice.induced_velocity(I, wake; - finite_core = true, - symmetric = true, - trailing_vortices = true, - xhat = [1, 0, 0]) - - @test isapprox(Vs, Vw) - -end - -@testset "Unsteady Vortex Lattice Method - Rectangular Wing" begin - - Uinf = 1.0 - AR = 4 - - # reference parameters - cref = 1.0 - bref = AR - Sref = bref*cref - rref = [0.0, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - # freestream parameters - alpha = 5.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - # geometry - xle = [0.0, 0.0] - yle = [-bref/2, bref/2] - zle = [0.0, 0.0] - chord = [cref, cref] - theta = [0.0, 0.0] - phi = [0.0, 0.0] - ns = 13 - nc = 4 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - symmetric = false - - # non-dimensional time - t = range(0.0, 10.0, step=0.2) - dt = [t[i+1]-t[i] for i = 1:length(t)-1] - - # create vortex rings - grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - surfaces = [surface] - - # run analysis - system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; - symmetric=symmetric) - - # extract forces at each time step - CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) -end - -@testset "Unsteady Vortex Lattice Method - Wing + Tail" begin - - # Unsteady Wing and Tail - - # wing - xle = [0.0, 0.2] - yle = [0.0, 5.0] - zle = [0.0, 1.0] - chord = [1.0, 0.6] - theta = [2.0*pi/180, 2.0*pi/180] - phi = [0.0, 0.0] - ns = 12 - nc = 1 - spacing_s = Uniform() - spacing_c = Uniform() - mirror = false - - # horizontal stabilizer - xle_h = [0.0, 0.14] - yle_h = [0.0, 1.25] - zle_h = [0.0, 0.0] - chord_h = [0.7, 0.42] - theta_h = [0.0, 0.0] - phi_h = [0.0, 0.0] - ns_h = 6 - nc_h = 1 - spacing_s_h = Uniform() - spacing_c_h = Uniform() - mirror_h = false - - # vertical stabilizer - xle_v = [0.0, 0.14] - yle_v = [0.0, 0.0] - zle_v = [0.0, 1.0] - chord_v = [0.7, 0.42] - theta_v = [0.0, 0.0] - phi_v = [0.0, 0.0] - ns_v = 5 - nc_v = 1 - spacing_s_v = Uniform() - spacing_c_v = Uniform() - mirror_v = false - - # adjust chord lengths to match AVL (which uses chord length in the x-direction) - chord = @. chord/cos(theta) - chord_h = @. chord_h/cos(theta_h) - chord_v = @. chord_v/cos(theta_v) - - Sref = 9.0 - cref = 0.9 - bref = 10.0 - rref = [0.5, 0.0, 0.0] - Vinf = 1.0 - ref = Reference(Sref, cref, bref, rref, Vinf) - - alpha = 5.0*pi/180 - beta = 0.0 - Omega = [0.0; 0.0; 0.0] - fs = Freestream(Vinf, alpha, beta, Omega) - - symmetric = [true, true, false] - - # horseshoe vortices - wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; - mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - - hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; - mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) - translate!(htail, [4.0, 0.0, 0.0]) - - vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; - mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) - translate!(vtail, [4.0, 0.0, 0.0]) - - surfaces = [wing, htail, vtail] - surface_id = [1, 2, 3] - - # t - t = range(0.0, 10.0, step=0.2) - dt = t[2:end] - t[1:end-1] - - system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; symmetric) - - # extract forces at each time step - CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) -end +# end \ No newline at end of file From 8833ebd5a29bd66177bc28de12c35f14eac83201 Mon Sep 17 00:00:00 2001 From: Ryan Anderson Date: Thu, 21 Mar 2024 14:11:08 -0600 Subject: [PATCH 3/3] add fmm acceleration --- src/VortexLattice.jl | 4 +- src/analyses.jl | 118 +- src/circulation.jl | 45 +- src/fmm.jl | 159 ++- src/induced.jl | 30 + src/nearfield.jl | 311 +++-- src/panel.jl | 92 +- src/system.jl | 31 +- src/wake.jl | 189 +-- test/runtests.jl | 2802 +++++++++++++++++++++--------------------- 10 files changed, 2082 insertions(+), 1699 deletions(-) diff --git a/src/VortexLattice.jl b/src/VortexLattice.jl index 18b4076..4acd003 100644 --- a/src/VortexLattice.jl +++ b/src/VortexLattice.jl @@ -10,7 +10,7 @@ using WriteVTK const RHO = 1.0 include("panel.jl") -export SurfacePanel, WakePanel, TrefftzPanel +export SurfacePanel, WakePanel, TrefftzPanel, FastMultipolePanel export reflect, set_normal include("wake.jl") @@ -37,6 +37,8 @@ include("system.jl") export System export PanelProperties, get_surface_properties +include("fmm.jl") + include("analyses.jl") export steady_analysis, steady_analysis! export unsteady_analysis, unsteady_analysis! diff --git a/src/analyses.jl b/src/analyses.jl index cf8b490..3119fca 100644 --- a/src/analyses.jl +++ b/src/analyses.jl @@ -75,7 +75,8 @@ function steady_analysis!(system, surfaces, ref, fs; fcore = (c, Δs) -> 1e-3, calculate_influence_matrix = true, near_field_analysis = true, - derivatives = true) + derivatives = true, + vertical_segments = false) # number of surfaces nsurf = length(surfaces) @@ -181,7 +182,7 @@ function steady_analysis!(system, surfaces, ref, fs; wake_finite_core = wake_finite_core, wake_shedding_locations = nothing, # shedding location at trailing edge trailing_vortices = trailing_vortices, - xhat = xhat) + xhat = xhat, vertical_segments = vertical_segments) else near_field_forces!(properties, surfaces, wakes, ref, fs, Γ; dΓdt = nothing, # no unsteady forces @@ -194,7 +195,7 @@ function steady_analysis!(system, surfaces, ref, fs; wake_finite_core = wake_finite_core, wake_shedding_locations = nothing, # shedding location at trailing edge trailing_vortices = trailing_vortices, - xhat = xhat) + xhat = xhat, vertical_segments = vertical_segments) end end @@ -241,7 +242,7 @@ step. panels in the system. Defaults to `zeros(N)` where `N` is the total number of surface panels in `surfaces`. - `nwake`: Maximum number of wake panels in the chordwise direction for each - surface. Defaults to `length(dx)` for all surfaces. + surface. Defaults to `length(dt)` for all surfaces. - `surface_id`: Surface ID for each surface. The finite core model is disabled when calculating the influence of surfaces/wakes that share the same ID. - `wake_finite_core`: Flag for each wake indicating whether the finite core @@ -264,10 +265,12 @@ unsteady_analysis # same grids/surfaces at each time step function unsteady_analysis(surfaces::AbstractVector{<:AbstractArray}, ref, fs, dt; - nwake = fill(length(dt), length(surfaces)), kwargs...) + nwake = fill(length(dt), length(surfaces)), + fmm_toggle=false, fmm_direct=false, fmm_p=4, fmm_ncrit=20, fmm_theta=0.4, + kwargs...) # pre-allocate system storage - system = System(surfaces; nw = nwake) + system = System(surfaces; nw = nwake, fmm_toggle, fmm_direct, fmm_p, fmm_ncrit, fmm_theta) return unsteady_analysis!(system, surfaces, ref, fs, dt; kwargs..., nwake, calculate_influence_matrix = true) @@ -275,10 +278,12 @@ end # different grids/surfaces at each time step function unsteady_analysis(surfaces::AbstractVector{<:AbstractVector{<:AbstractArray}}, - ref, fs, dt; nwake = fill(length(dt), length(surfaces[1])), kwargs...) + ref, fs, dt; nwake = fill(length(dt), length(surfaces[1])), + fmm_toggle=false, fmm_direct=false, fmm_p=4, fmm_ncrit=20, fmm_theta=0.4, + kwargs...) # pre-allocate system storage - system = System(surfaces[1]; nw = nwake) + system = System(surfaces[1]; nw = nwake, fmm_toggle, fmm_direct, fmm_p, fmm_ncrit, fmm_theta) return unsteady_analysis!(system, surfaces, ref, fs, dt; kwargs..., nwake, calculate_influence_matrix = true) @@ -301,7 +306,8 @@ function unsteady_analysis!(system, surfaces, ref, fs, dt; save = 1:length(dt), calculate_influence_matrix = true, near_field_analysis = true, - derivatives = true) + derivatives = true, + vertical_segments = false) # float number type TF = eltype(system) @@ -406,7 +412,8 @@ function unsteady_analysis!(system, surfaces, ref, fs, dt; eta = 0.25, calculate_influence_matrix = first_step && calculate_influence_matrix, near_field_analysis = it in save || (last_step && near_field_analysis), - derivatives = last_step && derivatives) + derivatives = last_step && derivatives, vertical_segments = vertical_segments, + fmm_velocity_probes = system.fmm_toggle ? system.fmm_velocity_probes : nothing) else propagate_system!(system, fs[it], dt[it]; additional_velocity, @@ -415,7 +422,8 @@ function unsteady_analysis!(system, surfaces, ref, fs, dt; eta = 0.25, calculate_influence_matrix = first_step && calculate_influence_matrix, near_field_analysis = it in save || (last_step && near_field_analysis), - derivatives = last_step && derivatives) + derivatives = last_step && derivatives, vertical_segments = vertical_segments, + fmm_velocity_probes = system.fmm_toggle ? system.fmm_velocity_probes : nothing) end # increment wake panel counter for each surface @@ -484,7 +492,8 @@ function propagate_system!(system, surfaces, fs, dt; eta, calculate_influence_matrix, near_field_analysis, - derivatives) + derivatives, vertical_segments, + fmm_velocity_probes) # NOTE: Each step models the transition from `t = t[it]` to `t = [it+1]` # (e.g. the first step models from `t = 0` to `t = dt`). Properties are @@ -520,6 +529,11 @@ function propagate_system!(system, surfaces, fs, dt; Vh = system.Vh Vv = system.Vv Vte = system.Vte + fmm_toggle = system.fmm_toggle + fmm_panels = system.fmm_panels + fmm_p = system.fmm_p + fmm_ncrit = system.fmm_ncrit + fmm_theta = system.fmm_theta # check if the surfaces are moving surface_motion = !isnothing(surfaces) @@ -566,6 +580,7 @@ function propagate_system!(system, surfaces, fs, dt; # update the wake shedding location for this time step # (based on freestream/kinematic/other velocity only) + # become the top corners of the next generation of wake update_wake_shedding_locations!(wakes, wake_shedding_locations, current_surfaces, ref, fs, dt, additional_velocity, Vte, nwake, eta) @@ -586,13 +601,59 @@ function propagate_system!(system, surfaces, fs, dt; wake_shedding_locations = wake_shedding_locations, trailing_vortices = trailing_vortices) - # calculate RHS <-- add FMM here: wake-on-all + # wake-on-all + if fmm_toggle # wake on all surfaces + # compile sources + update_fmm_panels!(fmm_panels, wakes, nwake) + + # preallocate probes + n_surface_panels = 0 + n_surface_filaments = 0 + for surface in current_surfaces + nc, ns = size(surface) + n_surface_panels += nc * ns # control points + n_surface_filaments += nc * ns + nc * ns + nc # top (nc*ns) + left (nc*ns) + right (nc) + end + + # n_wake_shedding_locations = 0 + # for wake_shedding_location in wake_shedding_locations + # n_wake_shedding_locations += length(wake_shedding_location) + # end + + n_wake_corners = 0 + for (nc, wake) in zip(nwake, wakes) + if nc > 0 + ns = size(wake, 2) + n_wake_corners += (nc+1)*(ns+1) + end + end + + update_n_probes!(fmm_velocity_probes, n_surface_panels + n_surface_filaments + n_wake_corners) + + # reset to zero + reset!(fmm_velocity_probes) + + # compile targets + update_probes!(fmm_velocity_probes, current_surfaces, 0) # control points and filament centers + update_probes!(fmm_velocity_probes, wakes, nwake, n_surface_panels + n_surface_filaments) # wake corners + + # run FMM + if length(fmm_panels) > 0 + if system.fmm_direct + FLOWFMM.direct!(fmm_velocity_probes, system) + else + FLOWFMM.fmm!(fmm_velocity_probes, system; expansion_order=fmm_p, n_per_branch_source=fmm_ncrit, n_per_branch_target=fmm_ncrit, theta=fmm_theta, ndivisions_source=10, ndivisions_target=10) + end + end + end + + # calculate RHS if derivatives normal_velocity_derivatives!(w, dw, current_surfaces, wakes, - ref, fs; additional_velocity, Vcp, symmetric, nwake, - surface_id, wake_finite_core, trailing_vortices, xhat) - else - normal_velocity!(w, current_surfaces, wakes, ref, fs; + ref, fs, fmm_velocity_probes; additional_velocity, Vcp, symmetric, nwake, + surface_id, wake_finite_core, trailing_vortices, xhat) + else + normal_velocity!(w, current_surfaces, wakes, ref, fs, fmm_velocity_probes; additional_velocity, Vcp, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat) end @@ -612,19 +673,28 @@ function propagate_system!(system, surfaces, fs, dt; dΓdt ./= dt # divide by corresponding time step # compute transient forces on each panel (if necessary) - # also compute surface-on-all induced velocity if near_field_analysis + if fmm_toggle + # update sources + update_fmm_panels!(fmm_panels, current_surfaces, wake_shedding_locations, Γ) + # run FMM + if system.fmm_direct + FLOWFMM.direct!(fmm_velocity_probes, system) + else + FLOWFMM.fmm!(fmm_velocity_probes, system; expansion_order=fmm_p, n_per_branch_source=fmm_ncrit, n_per_branch_target=fmm_ncrit, theta=fmm_theta, ndivisions_source=10, ndivisions_target=10) + end + end if derivatives near_field_forces_derivatives!(properties, dproperties, - current_surfaces, wakes, ref, fs, Γ, dΓ; dΓdt, + current_surfaces, wakes, ref, fs, Γ, dΓ, fmm_velocity_probes; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, surface_id, wake_finite_core, wake_shedding_locations, - trailing_vortices, xhat) + trailing_vortices, xhat, vertical_segments) else near_field_forces!(properties, current_surfaces, wakes, - ref, fs, Γ; dΓdt, additional_velocity, Vh, Vv, + ref, fs, Γ, fmm_velocity_probes; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, surface_id, wake_finite_core, - wake_shedding_locations, trailing_vortices, xhat) + wake_shedding_locations, trailing_vortices, xhat, vertical_segments) end # save flag indicating that a near-field analysis has been performed @@ -635,11 +705,11 @@ function propagate_system!(system, surfaces, fs, dt; system.derivatives[] = derivatives # update wake velocities - # ---should already be done if I've done this right--- get_wake_velocities!(wake_velocities, current_surfaces, wakes, ref, fs, Γ, additional_velocity, Vte, symmetric, repeated_points, nwake, surface_id, wake_finite_core, - wake_shedding_locations, trailing_vortices, xhat) + wake_shedding_locations, trailing_vortices, xhat, + fmm_velocity_probes) # shed additional wake panel (and translate existing wake panels) shed_wake!(wakes, wake_shedding_locations, wake_velocities, diff --git a/src/circulation.jl b/src/circulation.jl index ac36afb..e7257f8 100644 --- a/src/circulation.jl +++ b/src/circulation.jl @@ -48,7 +48,7 @@ induced velocity from the wake panels. This forms the right hand side of the circulation linear system solve. """ -@inline function normal_velocity!(w, surfaces, wakes, ref, fs; additional_velocity, +@inline function normal_velocity!(w, surfaces, wakes, ref, fs, probes=nothing; additional_velocity, Vcp, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat) nsurf = length(surfaces) @@ -88,15 +88,19 @@ This forms the right hand side of the circulation linear system solve. end # velocity due to the wake panels - for jsurf = 1:length(wakes) - if nwake[jsurf] > 0 - V += induced_velocity(rcp, wakes[jsurf]; - finite_core = wake_finite_core[isurf] || (surface_id[isurf] != surface_id[jsurf]), - symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], - xhat = xhat) + if isnothing(probes) + for jsurf = 1:length(wakes) + if nwake[jsurf] > 0 + V += induced_velocity(rcp, wakes[jsurf]; + finite_core = wake_finite_core[isurf] || (surface_id[isurf] != surface_id[jsurf]), + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + end end + else # use probes to obtain FMM acceleration + V += probes.velocity[iw+i] end # get normal component @@ -122,7 +126,7 @@ respect to the freestream parameters. This forms the right hand side of the circulation linear system solve (and its derivatives). """ -@inline function normal_velocity_derivatives!(w, dw, surfaces, wakes, ref, fs; +@inline function normal_velocity_derivatives!(w, dw, surfaces, wakes, ref, fs, probes=nothing; additional_velocity, Vcp, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat) @@ -172,15 +176,20 @@ This forms the right hand side of the circulation linear system solve (and its d end # velocity due to the wake panels - for jsurf = 1:length(wakes) - if nwake[jsurf] > 0 - V += induced_velocity(rcp, wakes[jsurf]; - finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]), - symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], - xhat = xhat) + if isnothing(probes) + for jsurf = 1:length(wakes) + if nwake[jsurf] > 0 + V += induced_velocity(rcp, wakes[jsurf]; + finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]), + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + end end + else # use probes for FMM acceleration + V += probes.velocity[iw+i] + @assert probes.position[iw+i] == rcp end # right hand side vector diff --git a/src/fmm.jl b/src/fmm.jl index 5311204..95a3262 100644 --- a/src/fmm.jl +++ b/src/fmm.jl @@ -4,12 +4,12 @@ Base.getindex(system::System, i, ::FLOWFMM.Position) = system.fmm_panels[i].rcp Base.getindex(system::System, i, ::FLOWFMM.Radius) = system.fmm_panels[i].radius -# Base.getindex(system::System, i, ::FLOWFMM.VectorPotential) = view(g.potential,2:4,i) -# Base.getindex(system::System, i, ::FLOWFMM.ScalarPotential) = g.potential[1,i] -Base.getindex(system::System, i, ::FLOWFMM.Velocity) = system.fmm_Vcp[i] -# Base.getindex(system::System, i, ::FLOWFMM.VelocityGradient) = reshape(view(g.potential,i_VELOCITY_GRADIENT,i),3,3) +# Base.getindex(system::System{TF}, i, ::FLOWFMM.VectorPotential) where TF = zero(SVector{3,TF}) +# Base.getindex(system::System{TF}, i, ::FLOWFMM.ScalarPotential) where TF = zero(TF) +Base.getindex(system::System{TF}, i, ::FLOWFMM.Velocity) where TF = zero(SVector{3,TF}) +# Base.getindex(system::System{TF}, i, ::FLOWFMM.VelocityGradient) where TF = zero(SMatrix{3,3,TF,9}) Base.getindex(system::System, i, ::FLOWFMM.ScalarStrength) = system.fmm_panels[i].gamma -Base.getindex(system::System, i, ::FLOWFMM.Body) = system.fmm_panels[i], system.fmm_Vcp[i] +Base.getindex(system::System, i, ::FLOWFMM.Body) = system.fmm_panels[i]#, system.fmm_Vcp[i] function Base.getindex(system::System, i, ::FLOWFMM.Vertex, i_vertex) if i_vertex == 1 return system.fmm_panels[i].rtl @@ -23,33 +23,33 @@ function Base.getindex(system::System, i, ::FLOWFMM.Vertex, i_vertex) end Base.getindex(system::System, i, ::FLOWFMM.Normal) = system.fmm_panels[i].ncp function Base.setindex!(system::System, val, i, ::FLOWFMM.Body) - panel, velocity = val - system.fmm_panels[i] = panel - system.fmm_Vcp[i] = velocity + # panel, velocity = val + system.fmm_panels[i] = val + # system.fmm_Vcp[i] = velocity return nothing end # function Base.setindex!(system::System, val, i, ::FLOWFMM.ScalarPotential) -# g.potential[i_POTENTIAL[1],i] = val +# return nothing # end # function Base.setindex!(system::System, val, i, ::FLOWFMM.VectorPotential) # g.potential[i_POTENTIAL[2:4],i] .= val # end -function Base.setindex!(system::System, val, i, ::FLOWFMM.Velocity) - system.fmm_Vcp[i] = val -end +# function Base.setindex!(system::System, val, i, ::FLOWFMM.Velocity) +# system.fmm_Vcp[i] = val +# end # function Base.setindex!(system::System, val, i, ::FLOWFMM.VelocityGradient) # reshape(g.potential[i_VELOCITY_GRADIENT,i],3,3) .= val # end FLOWFMM.get_n_bodies(system::System) = length(system.fmm_panels) -Base.eltype(::System{TF}) where TF = TF +# Base.eltype(::System{TF}) where TF = TF -FLOWFMM.buffer_element(system::System) = deepcopy(system.fmm_panels[1]), deepcopy(system.fmm_Vcp[1]) +FLOWFMM.buffer_element(system::System) = deepcopy(system.fmm_panels[1])#, deepcopy(system.fmm_Vcp[1]) FLOWFMM.B2M!(system::System, branch, bodies_index, harmonics, expansion_order) = - FLOWFMM.B2M!_quadpanel(system, branch, bodies_index, harmonics, system.fmm_p, FLOWFMM.UniformNormalDipolePanel()) + FLOWFMM.B2M!_quadpanel(system, branch, bodies_index, harmonics, Val(system.fmm_p), FLOWFMM.UniformNormalDipolePanel()) # no velocity gradient needed if the target is a `::VortexLattice.System` or a FLOWFMM.ProbeSystem without velocity gradient -function FLOWFMM.direct!(target_system::Union{System,FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,Vector{SVector{3,<:Any}},Nothing}}, target_index, source_system::System{TF}, source_index) where TF +function FLOWFMM.direct!(target_system::FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,<:Any,Nothing}, target_index, source_system::System{TF}, source_index) where TF # loop over targets for j_target in target_index rcp = target_system[j_target,FLOWFMM.POSITION] @@ -74,17 +74,17 @@ function FLOWFMM.direct!(target_system::Union{System,FLOWFMM.ProbeSystem{<:Any,< r4 = rcp - r21 # add induced velocity, ignoring finite core size for now - finite_core = false + finite_core = true this_V = zero(SVector{3,TF}) - this_V += bound_induced_velocity(r1, r2, finite_core, core_size) - this_V += bound_induced_velocity(r2, r3, finite_core, core_size) - this_V += bound_induced_velocity(r3, r4, finite_core, core_size) - this_V += bound_induced_velocity(r4, r1, finite_core, core_size) + this_V += bound_induced_velocity_fmm(r1, r2, finite_core, core_size) + this_V += bound_induced_velocity_fmm(r2, r3, finite_core, core_size) + this_V += bound_induced_velocity_fmm(r3, r4, finite_core, core_size) + this_V += bound_induced_velocity_fmm(r4, r1, finite_core, core_size) V += this_V * gamma end # add to target - target_system[j_target,FLOWFMM.VELOCITY] = target_system[j_target,FLOWFMM.VELOCITY] + V + target_system[j_target,FLOWFMM.VELOCITY] += V end end @@ -135,4 +135,117 @@ function FLOWFMM.direct!(target_system, target_index, source_system::System{TF}, target_system[j_target,FLOWFMM.VELOCITY] = target_system[j_target,FLOWFMM.VELOCITY] + V target_system[j_target,FLOWFMM.VELOCITY_GRADIENT] = target_system[j_target,FLOWFMM.VELOCITY_GRADIENT] + Vgrad end -end \ No newline at end of file +end + +##### +##### probes +##### +function update_n_probes!(probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,TV,Nothing}, n) where {TF,TV} + resize!(probes.position, n) + resize!(probes.velocity, n) +end + +# function update_n_probes!(probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,TV,TVG}, n) where {TF,TV,TVG} +# resize!(probes.position, n) +# resize!(probes.velocity, n) +# resize!(probes.velocity_gradient, n) +# end + +function reset!(probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,TV,Nothing}) where {TF,TV} + for i in eachindex(probes.velocity) + probes.velocity[i] = zero(eltype(probes.velocity)) + end +end + +# function reset!(probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,TV,TVG}) where {TF,TV,TVG} +# for i in eachindex(probes.velocity) +# probes.velocity[i] = zero(eltype(probes.velocity)) +# probes.velocity_gradient[i] = zero(eltype(probes.velocity_gradient)) +# end +# end + +"Places probes at the control point of all surface panels, and additional probes at the center of all surface filaments" +function update_probes!(fmm_velocity_probes::FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,<:Any,<:Any}, surfaces::Vector{Matrix{SurfacePanel{TF}}}, i_start) where TF # wake corners + i_probe = 1 + + # place probes at control points + for surface in surfaces + for panel in surface + fmm_velocity_probes.position[i_probe + i_start] = panel.rcp + i_probe += 1 + end + end + + # place probes at vortex filament centers (for force computation) + for surface in surfaces + nc, ns = size(surface) + for i in 1:ns + for j in 1:nc + + # top + fmm_velocity_probes.position[i_probe + i_start] = surface[j,i].rtc + i_probe += 1 + + # left + fmm_velocity_probes.position[i_probe + i_start] = left_center(surface[j,i]) + i_probe += 1 + end + + # NOTE: As the trailing edge panels have the same strength as the wake transition panel, their bottom vortex + # has zero strength and can be neglected when computing the force. Similar behavior to a Kutta condition. + # For now, the bottom vortex is ignored. + # # bottom + # fmm_velocity_probes.position[i_probe + i_start] = panel.rbc + # i_probe += 1 + end + + for j in 1:nc + # right + fmm_velocity_probes.position[i_probe + i_start] = right_center(surface[j,ns]) + i_probe += 1 + end + end +end + +"Places probes at the corners of all wake panels" +function update_probes!(fmm_velocity_probes::FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,<:Any,<:Any}, wakes::Vector{Matrix{WakePanel{TF}}}, nwake, i_start) where TF # wake corners + i_corner = i_start + 1 + for (nc, wake) in zip(nwake, wakes) + if nc > 0 + ns = size(wake,2) # number of spanwise panels + + # get wake_shedding_locations (first row of wake points) + for i in 1:ns + fmm_velocity_probes.position[i_corner] = wake[1,i].rtl + i_corner += 1 + end + fmm_velocity_probes.position[i_corner] = wake[1,ns].rtr + i_corner += 1 + + # bottom left of each panel + for i in 1:ns + for j in 1:nc + fmm_velocity_probes.position[i_corner] = wake[j,i].rbl + i_corner += 1 + end + end + + # bottom right of the last column + for j in 1:nc + fmm_velocity_probes.position[i_corner] = wake[j,ns].rbr + i_corner += 1 + end + end + end +end + +# "Places probes at all wake shedding locations" +# function update_probes!(fmm_velocity_probes::FLOWFMM.ProbeSystem{<:Any,<:Any,<:Any,<:Any,<:Any}, wake_shedding_locations::Vector{Vector{SVector{3,TF}}}, i_start) where TF # wake transition points +# i_shed = 1 +# for wsl in wake_shedding_locations +# for location in wsl +# fmm_velocity_probes.position[i_shed + i_start] = location +# i_shed += 1 +# end +# end +# end diff --git a/src/induced.jl b/src/induced.jl index f936e41..e7bd14a 100644 --- a/src/induced.jl +++ b/src/induced.jl @@ -27,6 +27,36 @@ relative to the end of the bound vortex return Vhat end +@inline function bound_induced_velocity_fmm(r1, r2, finite_core, core_size; epsilon=sqrt(eps())) + + # parameters + nr1 = norm(r1) + nr2 = norm(r2) + nr1nr2 = nr1*nr2 + rcross = cross(r1, r2) + rdot = dot(r1, r2) + + # check if evaluation point is colinear with the bound vortex + if norm(rcross) < epsilon # colinear if true + if isapprox(rdot, -nr1nr2; atol=epsilon) # at the midpoint, so return zero + return zero(typeof(r1)) + elseif rdot <= 0.0 # coincident with the filament so use the finite core model + r1s, r2s, εs = nr1^2, nr2^2, core_size^2 + f1 = rcross/(r1s*r2s - rdot^2 + εs*(r1s + r2s - 2*nr1nr2)) + f2 = (r1s - rdot)/sqrt(r1s + εs) + (r2s - rdot)/sqrt(r2s + εs) + Vhat = (f1*f2)/(4*pi) + return Vhat + end + end + + f1 = rcross/(nr1nr2 + rdot) + f2 = (1/nr1 + 1/nr2) + + Vhat = (f1*f2)/(4*pi) + + return Vhat +end + function bound_velocity_gradient(r1::SVector{3,TF},r2) where TF # zeta r1norm = norm(r1) diff --git a/src/nearfield.jl b/src/nearfield.jl index ffcae90..4eac4ae 100644 --- a/src/nearfield.jl +++ b/src/nearfield.jl @@ -5,15 +5,21 @@ Calculate local panel forces in the body frame. """ -@inline function near_field_forces!(props, surfaces, wakes, ref, fs, Γ; +@inline function near_field_forces!(props, surfaces, wakes, ref, fs, Γ, probes=nothing; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, surface_id, - wake_finite_core, wake_shedding_locations, trailing_vortices, xhat) + wake_finite_core, wake_shedding_locations, trailing_vortices, xhat, vertical_segments=false) # number of surfaces nsurf = length(surfaces) # loop through receiving surfaces iΓ = 0 # index for accessing Γ + i_probe = 1 # index for accessing probes + for surface in surfaces + nc, ns = size(surface) + i_probe += nc*ns + end + for isurf = 1:nsurf receiving = surfaces[isurf] @@ -48,57 +54,65 @@ Calculate local panel forces in the body frame. Vi += Vh[isurf][i] end - # induced velocity from surfaces and wakes - jΓ = 0 # index for accessing Γ - for jsurf = 1:nsurf - - # number of panels on sending surface - sending = surfaces[jsurf] - Ns = length(sending) - - # see if wake panels are being used - wake_panels = nwake[jsurf] > 0 - - # check if we need to shift shedding locations - if isnothing(wake_shedding_locations) - shedding_locations = nothing - else - shedding_locations = wake_shedding_locations[jsurf] - end - - # extract circulation values corresonding to the sending surface - vΓ = view(Γ, jΓ+1:jΓ+Ns) - - # induced velocity from this surface - if isurf == jsurf - # induced velocity on self - Vi += induced_velocity(I, surfaces[jsurf], vΓ; - finite_core = surface_id[isurf] != surface_id[jsurf], - wake_shedding_locations = shedding_locations, - symmetric = symmetric[jsurf], - trailing_vortices = trailing_vortices[jsurf] && !wake_panels, - xhat = xhat) - else - # induced velocity on another surface - Vi += induced_velocity(rc, surfaces[jsurf], vΓ; - finite_core = surface_id[isurf] != surface_id[jsurf], - wake_shedding_locations = shedding_locations, - symmetric = symmetric[jsurf], - trailing_vortices = trailing_vortices[jsurf] && !wake_panels, - xhat = xhat) - end + this_Vind = zero(SVector{3,Float64}) - # induced velocity from corresponding wake - if wake_panels - Vi += induced_velocity(rc, wakes[jsurf]; - finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]), - symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], - xhat = xhat) + # induced velocity from surfaces and wakes + if isnothing(probes) # no fast multipole acceleration + jΓ = 0 # index for accessing Γ + for jsurf = 1:nsurf + + # number of panels on sending surface + sending = surfaces[jsurf] + Ns = length(sending) + + # see if wake panels are being used + wake_panels = nwake[jsurf] > 0 + + # check if we need to shift shedding locations + if isnothing(wake_shedding_locations) + shedding_locations = nothing + else + shedding_locations = wake_shedding_locations[jsurf] + end + + # extract circulation values corresonding to the sending surface + vΓ = view(Γ, jΓ+1:jΓ+Ns) + + # induced velocity from this surface + if isurf == jsurf + # induced velocity on self + Vi += induced_velocity(I, surfaces[jsurf], vΓ; + finite_core = surface_id[isurf] != surface_id[jsurf], + wake_shedding_locations = shedding_locations, + symmetric = symmetric[jsurf], + trailing_vortices = trailing_vortices[jsurf] && !wake_panels, + xhat = xhat) + + else + # induced velocity on another surface + Vi += induced_velocity(rc, surfaces[jsurf], vΓ; + finite_core = surface_id[isurf] != surface_id[jsurf], + wake_shedding_locations = shedding_locations, + symmetric = symmetric[jsurf], + trailing_vortices = trailing_vortices[jsurf] && !wake_panels, + xhat = xhat) + + end + + # induced velocity from corresponding wake + if wake_panels # no fast multipole acceleration + Vi += induced_velocity(rc, wakes[jsurf]; + finite_core = wake_finite_core[jsurf] || (surface_id[isurf] != surface_id[jsurf]), + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + end + + jΓ += Ns # increment Γ index for sending panels end - - jΓ += Ns # increment Γ index for sending panels + else # fast multipole acceleration + Vi += probes.velocity[i_probe] end # steady part of Kutta-Joukowski theorem @@ -144,6 +158,14 @@ Calculate local panel forces in the body frame. # for the vertical segments because its influence is likely negligible # once we take the cross product with the bound vortex vector. This # is also assumed in AVL. This could change in the future. + # BUT.... if we did, and we wanted to use fast multipole acceleration, + # this is what it would look like: + # ACTUALLY: We DO include induced velocity in the effective velocity IF THE FMM IS USED + # as its cost is small at this point, and it provides greater generality. + + if !isnothing(probes) && vertical_segments + Veff += probes.velocity[i_probe + 1] # left bound vortex is always 1 index ahead of the bound vortex + end # steady part of Kutta-Joukowski theorem Γli = Γ[iΓ+i] @@ -174,6 +196,17 @@ Calculate local panel forces in the body frame. # for the vertical segments because its influence is likely negligible # once we take the cross product with the bound vortex vector. This # is also assumed in AVL. This could change in the future. + # BUT.... if we did, and we wanted to use fast multipole acceleration, + # what it would look like: + # ACTUALLY... we DO include the induced velocity IF THE FMM IS USED + # as its cost is negligible. + + if !isnothing(probes) && vertical_segments + i_probe_right = i_probe + nr1<<1 + 1 + I[2] == ns && (i_probe_right -= I[1]) + Veff += probes.velocity[i_probe_right] # right bound vortex is always 2*nc+1 indices ahead of the bound vortex + # unless we are in the final column, in which case it is I[1] less + end # steady part of Kutta-Joukowski theorem Γri = Γ[iΓ+i] @@ -185,10 +218,17 @@ Calculate local panel forces in the body frame. props[isurf][i] = PanelProperties(Γ[iΓ+i]/ref.V, Vi/ref.V, Fbi/(q*ref.S), Fbli/(q*ref.S), Fbri/(q*ref.S)) + + + # increment probe index + i_probe += 2 # by 2 since each column has top and left bound vortices end # increment Γ index for receiving panels iΓ += nr + + # increment probe index + i_probe += nr1 # by nc to account for the column of right bound vortices end return props @@ -205,8 +245,8 @@ the local panel forces with respect to the freestream variables. near_field_forces_derivatives! @inline function near_field_forces_derivatives!(props, dprops, surfaces, wakes, - ref, fs, Γ, dΓ; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, - surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat) + ref, fs, Γ, dΓ, probes=nothing; dΓdt, additional_velocity, Vh, Vv, symmetric, nwake, + surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat, vertical_segments=false) # unpack derivatives props_a, props_b, props_p, props_q, props_r = dprops @@ -217,6 +257,11 @@ near_field_forces_derivatives! # loop through receiving surfaces iΓ = 0 # index for accessing Γ + i_probe = 1 # index for accessing probes + for surface in surfaces # filament probes begin after control points of each surface panel + i_probe += length(surface) + end + for isurf = 1:nsurf receiving = surfaces[isurf] @@ -252,74 +297,78 @@ near_field_forces_derivatives! end # induced velocity from surfaces and wakes - jΓ = 0 # index for accessing Γ - for jsurf = 1:nsurf - - # number of panels on sending surface - sending = surfaces[jsurf] - Ns = length(sending) - - # see if wake panels are being used - wake_panels = nwake[jsurf] > 0 - - # check if we need to shift shedding locations - if isnothing(wake_shedding_locations) - shedding_locations = nothing - else - shedding_locations = wake_shedding_locations[jsurf] - end - - # extract circulation values corresonding to the sending surface - vΓ = view(Γ, jΓ+1:jΓ+Ns) - - vΓ_a = view(Γ_a, jΓ+1:jΓ+Ns) - vΓ_b = view(Γ_b, jΓ+1:jΓ+Ns) - vΓ_p = view(Γ_p, jΓ+1:jΓ+Ns) - vΓ_q = view(Γ_q, jΓ+1:jΓ+Ns) - vΓ_r = view(Γ_r, jΓ+1:jΓ+Ns) - - vdΓ = (vΓ_a, vΓ_b, vΓ_p, vΓ_q, vΓ_r) - - # induced velocity from this surface - if isurf == jsurf - # induced velocity on self - Vind, dVind = induced_velocity_derivatives(I, surfaces[jsurf], vΓ, vdΓ; - finite_core = surface_id[isurf] != surface_id[jsurf], - wake_shedding_locations = shedding_locations, - symmetric = symmetric[jsurf], - trailing_vortices = trailing_vortices[jsurf] && !wake_panels, - xhat = xhat) - else - # induced velocity on another surface - Vind, dVind = induced_velocity_derivatives(rc, surfaces[jsurf], vΓ, vdΓ; - finite_core = surface_id[isurf] != surface_id[jsurf], - wake_shedding_locations = shedding_locations, - symmetric = symmetric[jsurf], - trailing_vortices = trailing_vortices[jsurf] && !wake_panels, - xhat = xhat) - end - - Vind_a, Vind_b, Vind_p, Vind_q, Vind_r = dVind - - Vi += Vind - - Vi_a += Vind_a - Vi_b += Vind_b - Vi_p += Vind_p - Vi_q += Vind_q - Vi_r += Vind_r - - # induced velocity from corresponding wake - if wake_panels - Vi += induced_velocity(rc, wakes[jsurf]; - finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], - symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], - xhat = xhat) + if isnothing(probes) + jΓ = 0 # index for accessing Γ + for jsurf = 1:nsurf + + # number of panels on sending surface + sending = surfaces[jsurf] + Ns = length(sending) + + # see if wake panels are being used + wake_panels = nwake[jsurf] > 0 + + # check if we need to shift shedding locations + if isnothing(wake_shedding_locations) + shedding_locations = nothing + else + shedding_locations = wake_shedding_locations[jsurf] + end + + # extract circulation values corresonding to the sending surface + vΓ = view(Γ, jΓ+1:jΓ+Ns) + + vΓ_a = view(Γ_a, jΓ+1:jΓ+Ns) + vΓ_b = view(Γ_b, jΓ+1:jΓ+Ns) + vΓ_p = view(Γ_p, jΓ+1:jΓ+Ns) + vΓ_q = view(Γ_q, jΓ+1:jΓ+Ns) + vΓ_r = view(Γ_r, jΓ+1:jΓ+Ns) + + vdΓ = (vΓ_a, vΓ_b, vΓ_p, vΓ_q, vΓ_r) + + # induced velocity from this surface + if isurf == jsurf + # induced velocity on self + Vind, dVind = induced_velocity_derivatives(I, surfaces[jsurf], vΓ, vdΓ; + finite_core = surface_id[isurf] != surface_id[jsurf], + wake_shedding_locations = shedding_locations, + symmetric = symmetric[jsurf], + trailing_vortices = trailing_vortices[jsurf] && !wake_panels, + xhat = xhat) + else + # induced velocity on another surface + Vind, dVind = induced_velocity_derivatives(rc, surfaces[jsurf], vΓ, vdΓ; + finite_core = surface_id[isurf] != surface_id[jsurf], + wake_shedding_locations = shedding_locations, + symmetric = symmetric[jsurf], + trailing_vortices = trailing_vortices[jsurf] && !wake_panels, + xhat = xhat) + + end + + Vind_a, Vind_b, Vind_p, Vind_q, Vind_r = dVind + Vi += Vind + Vi_a += Vind_a + Vi_b += Vind_b + Vi_p += Vind_p + Vi_q += Vind_q + Vi_r += Vind_r + + # induced velocity from corresponding wake + if wake_panels + Vi += induced_velocity(rc, wakes[jsurf]; + finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + end + + jΓ += Ns # increment Γ index for sending panels end - - jΓ += Ns # increment Γ index for sending panels + else # fast multipole acceleration + Vi += probes.velocity[i_probe] + @assert probes.position[i_probe] == rc end # steady part of Kutta-Joukowski theorem @@ -394,6 +443,12 @@ near_field_forces_derivatives! # for the vertical segments because its influence is likely negligible # once we take the cross product with the bound vortex vector. This # is also assumed in AVL. This could change in the future. + # BUT.... if we did, and we wanted to use fast multipole acceleration, + # this is what it would look like: + + if !isnothing(probes) && vertical_segments + Veff += probes.velocity[i_probe + 1] # left bound vortex is always 1 index ahead of the bound vortex + end # steady part of Kutta-Joukowski theorem Γli = Γ[iΓ+i] @@ -444,6 +499,15 @@ near_field_forces_derivatives! # for the vertical segments because its influence is likely negligible # once we take the cross product with the bound vortex vector. This # is also assumed in AVL. This could change in the future. + # BUT.... if we did, and we wanted to use fast multipole acceleration, + # this is what it would look like: + + if !isnothing(probes) && vertical_segments + i_probe_right = i_probe + nr1<<1 + 1 + I[2] == nr2 && (i_probe_right -= I[1]) + Veff += probes.velocity[i_probe_right] # right bound vortex is always 2*nc+1 indices ahead of the bound vortex + # unless we are in the final column, in which case it is I[1] less + end # steady part of Kutta-Joukowski theorem Γri = Γ[iΓ+i] @@ -482,10 +546,17 @@ near_field_forces_derivatives! Fbli_q/(q*ref.S), Fbri_q/(q*ref.S)) props_r[isurf][I] = PanelProperties(Γ_r[iΓ+i]/ref.V, Vi_r/ref.V, Fbi_r/(q*ref.S), Fbli_r/(q*ref.S), Fbri_r/(q*ref.S)) + + # increment probe index + i_probe += 2 # by 2 since each column has top and left bound vortices + # (note there is a final column containing the right bound vortices) end # increment Γ index for receiving panels iΓ += nr + + # increment probe index + i_probe += nr1 # by nc to account for the column of right bound vortices end return props, dprops diff --git a/src/panel.jl b/src/panel.jl index 34ef0a0..664a642 100644 --- a/src/panel.jl +++ b/src/panel.jl @@ -636,19 +636,105 @@ struct FastMultipolePanel{TF} rtr::SVector{3,TF} rbl::SVector{3,TF} rbr::SVector{3,TF} + rcp::SVector{3,TF} ncp::SVector{3,TF} radius::TF core_size::TF gamma::TF end -FastMultipolePanel(p::SurfacePanel{TF}, gamma) where TF = FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.ncp, max(norm(p.rtl-p.rbr)/2, norm(p.rbl-p.rtr)/2), p.core_size, gamma) +FastMultipolePanel(p::SurfacePanel{TF}, gamma) where TF = FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.rcp, p.ncp, max(norm(p.rtl-p.rbr)/2, norm(p.rbl-p.rtr)/2), p.core_size, gamma) function FastMultipolePanel(p::WakePanel{TF}) where TF v1 = p.rtr-p.rbl v2 = p.rtl-p.rbr + rcp_t = linearinterp(0.5,p.rtr,p.rtl) + rcp_b = linearinterp(0.5,p.rbr,p.rbl) + rcp = linearinterp(0.5,rcp_t,rcp_b) ncp = cross(v1, v2) ncp /= norm(ncp) - radius = max(norm(v1)/2, norm(v2)/2) - return FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, ncp, radius, p.core_size, p.gamma) + radius = max(norm(v1), norm(v2))/2 + return FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, rcp, ncp, radius, p.core_size, p.gamma) +end + +function FastMultipolePanel(p::SurfacePanel{TF}, gamma) where TF + v1 = p.rtr-p.rbl + v2 = p.rtl-p.rbr + radius = max(norm(v1), norm(v2))/2 + return FastMultipolePanel{TF}(p.rtl, p.rtr, p.rbl, p.rbr, p.rcp, p.ncp, radius, p.core_size, gamma) +end + +function update_fmm_panels!(fmm_panels, wake_panels, nwake) + # clear space for fmm panels + n_panels = 0 + for (nc, wake) in zip(nwake, wake_panels) + nc > 0 && (n_panels += size(wake,2) * nc) + end + resize!(fmm_panels, n_panels) + + # copy wake panels to fmm panels + i_fmm = 1 + for (nc, wake) in zip(nwake, wake_panels) + ns = size(wake, 2) + for i in 1:ns + for j in 1:nc + fmm_panels[i_fmm] = FastMultipolePanel(wake[j,i]) + i_fmm += 1 + end + end + end + + return fmm_panels +end + +function update_fmm_panels!(fmm_panels, surface_panels, wake_shedding_locations, gamma) + # clear space for fmm panels + n_panels = 0 + for surface in surface_panels + nc, ns = size(surface) + n_panels += (nc+1)*ns + end + resize!(fmm_panels, n_panels) + + # add surface panels + i_panel = 1 + for surface in surface_panels + for panel in surface + fmm_panels[i_panel] = VortexLattice.FastMultipolePanel(panel, gamma[i_panel]) + i_panel += 1 + end + end + + # add wake transition panels + i_start = 0 + for (wake_shedding_location,surface) in zip(wake_shedding_locations,surface_panels) + + # number of chord/spanwise panels + nc, ns = size(surface) + + # loop over transition panels + i_gamma = i_start + nc + for j in 1:ns + rtl = surface[nc,j].rbl + rtr = surface[nc,j].rbr + core_size = surface[nc,j].core_size + rbl = wake_shedding_location[j] + rbr = wake_shedding_location[j+1] + v1 = rtr-rbl + v2 = rtl-rbr + rcp_t = linearinterp(0.5,rtr,rtl) + rcp_b = linearinterp(0.5,rbr,rbl) + rcp = linearinterp(0.5,rcp_t,rcp_b) + ncp = cross(v1, v2) + ncp /= norm(ncp) + radius = max(norm(v1), norm(v2))/2 + fmm_panels[i_panel] = FastMultipolePanel(rtl,rtr,rbl,rbr,rcp,ncp,radius,core_size,gamma[i_gamma]) + i_gamma += nc + i_panel += 1 + end + + i_start += length(surface) + end + + return fmm_panels end \ No newline at end of file diff --git a/src/system.jl b/src/system.jl index 1b1108e..c4a17e1 100644 --- a/src/system.jl +++ b/src/system.jl @@ -73,10 +73,10 @@ Contains pre-allocated storage for internal system variables. - `Vv`: Velocity due to surface motion at the vertical bound vortex centers - `Vte`: Velocity due to surface motion at the trailing edge vertices - `dΓdt`: Derivative of the circulation strength with respect to non-dimensional time - - `fmm_panels`: Fast-access copies of all surface panels (including surface-wake transition panels) for fast multipole acceleration. - - `fmm_Vcp`: Velocity at the control points of `fmm_panels` - - `fmm_force_probes`: Locations and induced velocity of the center of the top vortices - - `fmm_wake_probes`: Locations, induced velocity, and induced velocity gradient of the wake particle shed locations + - `fmm_toggle`: Switch for activating fast multipole acceleration + - `fmm_direct`: Switch for activating the `FLOWFMM.direct!` function (for debugging purposes) + - `fmm_panels`: Fast-access copies of all surface panels (including surface-wake transition panels) for fast multipole acceleration + - `fmm_velocity_probes`: Probes for obtaining the induced velocity at surface control points, wake corners, and wake shedding locations - `fmm_p`: Multipole expansion order - `fmm_ncrit`: Maximum number of panels in the leaf level of the FMM - `fmm_theta`: Multipole acceptance criterion for the FMM, between 0 and 1 @@ -110,13 +110,13 @@ struct System{TF} Vv::Vector{Matrix{SVector{3, TF}}} Vte::Vector{Vector{SVector{3, TF}}} dΓdt::Vector{TF} - # fmm_panels::Vector{FastMultipolePanel{TF}} - # fmm_Vcp::Vector{SVector{3,TF}} - # fmm_velocity_probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,Vector{SVector{3,TF}},Nothing} - # fmm_gradient_probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,Vector{SVector{3,TF}},Vector{SMatrix{3,3,TF,9}}} - # fmm_p::Val{Int} - # fmm_ncrit::Int - # fmm_theta::Float64 + fmm_toggle::Bool + fmm_direct::Bool + fmm_panels::Vector{FastMultipolePanel{TF}} + fmm_velocity_probes::FLOWFMM.ProbeSystem{TF,Nothing,Nothing,Vector{SVector{3,TF}},Nothing} + fmm_p::Int + fmm_ncrit::Int + fmm_theta::Float64 end @inline Base.eltype(::Type{System{TF}}) where TF = TF @@ -192,7 +192,8 @@ variables - `nw`: Number of chordwise wake panels for each surface. Defaults to zero wake panels on each surface """ -function System(TF::Type, nc, ns; nw = zero(nc)) +function System(TF::Type, nc, ns; nw = zero(nc), + fmm_toggle=false, fmm_direct=false, fmm_p=4, fmm_ncrit=20, fmm_theta=0.4) @assert length(nc) == length(ns) == length(nw) @@ -230,12 +231,16 @@ function System(TF::Type, nc, ns; nw = zero(nc)) Vv = [fill((@SVector zeros(TF, 3)), nc[i], ns[i]+1) for i = 1:nsurf] Vte = [fill((@SVector zeros(TF, 3)), ns[i]+1) for i = 1:nsurf] dΓdt = zeros(TF, N) + fmm_panels = Vector{FastMultipolePanel{TF}}(undef,0) + fmm_toggle && sizehint!(fmm_panels, max(sum((nc .+1).*ns), sum(length.(wakes)))) + n_velocity_probes = fmm_toggle ? N + sum((nw .+1) .* (ns .+1)) + sum(ns .+1) : 0 # surface control points + wake corners + wake shedding locations + fmm_velocity_probes = FLOWFMM.ProbeSystem(n_velocity_probes, TF; velocity=true) return System{TF}(AIC, w, Γ, V, surfaces, properties, wakes, trefftz, reference, freestream, symmetric, nwake, surface_id, wake_finite_core, trailing_vortices, xhat, near_field_analysis, derivatives, dw, dΓ, dproperties, wake_shedding_locations, previous_surfaces, Vcp, Vh, - Vv, Vte, dΓdt) + Vv, Vte, dΓdt, fmm_toggle, fmm_direct, fmm_panels, fmm_velocity_probes, fmm_p, fmm_ncrit, fmm_theta) end """ diff --git a/src/wake.jl b/src/wake.jl index 27ff533..6d2e662 100644 --- a/src/wake.jl +++ b/src/wake.jl @@ -132,10 +132,18 @@ end """ @inline function get_wake_velocities!(wake_velocities, surfaces, wakes, ref, fs, Γ, additional_velocity, Vte, symmetric, repeated_points, nwake, - surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat) + surface_id, wake_finite_core, wake_shedding_locations, trailing_vortices, xhat, + probes=nothing) # number of surfaces nsurf = length(surfaces) + + # first probe + i_probe = 1 + for surface in surfaces # surface control points and filament centers + nc, ns = size(surface) + i_probe += nc*ns + 2*nc*ns + nc + end # loop through all surfaces for isurf = 1:nsurf @@ -147,20 +155,21 @@ end # velocity at the wake shedding locations for is = 1:ns+1 - # skip loop if true - continue_loop = false + # # skip loop if true + # continue_loop = false # check if this point is a duplicate, skip if it is - if (isurf, is) in keys(repeated_points) - for (jsurf, js) in repeated_points[(isurf, is)] - # NOTE: we assume that a point is not repeated on the same surface - if jsurf < isurf - wake_velocities[isurf][1, is] = wake_velocities[jsurf][1, js] - continue_loop = true - end - end - continue_loop && continue - end + # if (isurf, is) in keys(repeated_points) + # for (jsurf, js) in repeated_points[(isurf, is)] + # # NOTE: we assume that a point is not repeated on the same surface + # if jsurf < isurf + # wake_velocities[isurf][1, is] = wake_velocities[jsurf][1, js] + # # continue_loop = true + # continue + # end + # end + # # continue_loop && continue + # end # get vertex location rc = wake_shedding_locations[isurf][is] @@ -178,6 +187,9 @@ end # velocity at the trailing edge wake_velocities[isurf][1,is] = V + + # we're ignoring the wake- and surface-induced velocity here + i_probe += 1 end # velocity at all other wake vertices @@ -187,15 +199,15 @@ end for I in cr # check if this point is a duplicate, skip if it is - if (isurf, I[2]) in keys(repeated_points) - for (jsurf, js) in repeated_points[(isurf, I[2])] - # NOTE: we assume that a point is not repeated on the same surface - if jsurf < isurf - wake_velocities[isurf][I] = wake_velocities[jsurf][I[1], js] - continue - end - end - end + # if (isurf, I[2]) in keys(repeated_points) + # for (jsurf, js) in repeated_points[(isurf, I[2])] + # # NOTE: we assume that a point is not repeated on the same surface + # if jsurf < isurf + # wake_velocities[isurf][I] = wake_velocities[jsurf][I[1], js] + # continue + # end + # end + # end # get vertex location if I[1] <= nw && I[2] <= ns @@ -220,77 +232,82 @@ end end # induced velocity from each surface and wake - jΓ = 0 # index for accessing Γ - for jsurf = 1:nsurf - - # number of panels on sending surface - Ns = length(surfaces[jsurf]) - - # check if receiving point is repeated on the sending surface - if isurf == jsurf - # the surfaces are the same - same_surface = true - # vertex spanwise coordinate - js = I[2] - else - # the surfaces are different, but the point could still be a duplicate - if (isurf, I[2]) in keys(repeated_points) - # the vertex is duplicated - - # check if the vertex is on the sending surface - idx = findfirst(x -> x[1] == jsurf, repeated_points[(isurf, I[2])]) - - if isnothing(idx) - # the vertex is not duplicated on the sending surface - same_surface = false - else - # the vertex is duplicated on the sending surface - same_surface = true - # vertex spanwise coordinate on sending surface - js = repeated_points[(isurf, I[2])][idx][2] - end + if isnothing(probes) # no fmm acceleration + jΓ = 0 # index for accessing Γ + for jsurf = 1:nsurf + + # number of panels on sending surface + Ns = length(surfaces[jsurf]) + + # check if receiving point is repeated on the sending surface + if isurf == jsurf + # the surfaces are the same + same_surface = true + # vertex spanwise coordinate + js = I[2] else - # the vertex is not duplicated + # the surfaces are different, but the point could still be a duplicate + if (isurf, I[2]) in keys(repeated_points) + # the vertex is duplicated + + # check if the vertex is on the sending surface + idx = findfirst(x -> x[1] == jsurf, repeated_points[(isurf, I[2])]) + + if isnothing(idx) + # the vertex is not duplicated on the sending surface + same_surface = false + else + # the vertex is duplicated on the sending surface + same_surface = true + # vertex spanwise coordinate on sending surface + js = repeated_points[(isurf, I[2])][idx][2] + end + else + # the vertex is not duplicated - # the vertex is not on the sending surface - same_surface = false + # the vertex is not on the sending surface + same_surface = false + end end - end - # extract circulation values corresponding to the sending surface - vΓ = view(Γ, jΓ+1:jΓ+Ns) - - # induced velocity from this surface - wake_velocities[isurf][I] += induced_velocity(rc, surfaces[jsurf], vΓ; - finite_core = surface_id[isurf] != surface_id[jsurf], - wake_shedding_locations = wake_shedding_locations[jsurf], - symmetric = symmetric[jsurf], - trailing_vortices = false, - xhat = xhat) - - # add induced velocity from the wake - if same_surface - # vertex location on wake - J = CartesianIndex(I[1], js) - - # induced velocity from wake on its own vertex - wake_velocities[isurf][I] += induced_velocity(J, wakes[jsurf]; - finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], - symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], - xhat = xhat) - else - # induced velocity from wake on another wake's vertex - wake_velocities[isurf][I] += induced_velocity(rc, wakes[jsurf]; - finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], + # extract circulation values corresponding to the sending surface + vΓ = view(Γ, jΓ+1:jΓ+Ns) + + # induced velocity from this surface + wake_velocities[isurf][I] += induced_velocity(rc, surfaces[jsurf], vΓ; + finite_core = surface_id[isurf] != surface_id[jsurf], + wake_shedding_locations = wake_shedding_locations[jsurf], symmetric = symmetric[jsurf], - nc = nwake[jsurf], - trailing_vortices = trailing_vortices[jsurf], + trailing_vortices = false, xhat = xhat) - end - jΓ += Ns # increment Γ index + # add induced velocity from the wake + if same_surface + # vertex location on wake + J = CartesianIndex(I[1], js) + + # induced velocity from wake on its own vertex + wake_velocities[isurf][I] += induced_velocity(J, wakes[jsurf]; + finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + else + # induced velocity from wake on another wake's vertex + wake_velocities[isurf][I] += induced_velocity(rc, wakes[jsurf]; + finite_core = wake_finite_core[jsurf] || surface_id[isurf] != surface_id[jsurf], + symmetric = symmetric[jsurf], + nc = nwake[jsurf], + trailing_vortices = trailing_vortices[jsurf], + xhat = xhat) + end + + jΓ += Ns # increment Γ index + end + else # fmm acceleration + wake_velocities[isurf][I] += probes.velocity[i_probe] + i_probe += 1 end end end diff --git a/test/runtests.jl b/test/runtests.jl index 196ca81..ec48cb3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,1418 +18,1419 @@ function avl_normal_vector(dr, theta) return ncp / norm(ncp) # normal vector used by AVL end -# @testset "AVL - Run 1 - Wing with Uniform Spacing" begin +@testset "AVL - Run 1 - Wing with Uniform Spacing" begin -# # Simple Wing with Uniform Spacing - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() + # Simple Wing with Uniform Spacing -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + # vortex rings with symmetry + mirror = false + symmetric = true + + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.24324, atol=1e-3) + @test isapprox(CD, 0.00243, atol=1e-5) + @test isapprox(CDiff, 0.00245, atol=1e-5) + @test isapprox(Cm, -0.02252, atol=1e-4) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) + + # vortex rings with mirrored geometry + mirror = true + symmetric = false + + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.24324, atol=1e-3) + @test isapprox(CD, 0.00243, atol=1e-5) + @test isapprox(CDiff, 0.00245, atol=1e-5) + @test isapprox(Cm, -0.02252, atol=1e-4) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 2 - Wing with Cosine Spacing" begin + + # Run 2: Simple Wing with Cosine Spacing + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Cosine() + spacing_c = Uniform() + mirror = true + symmetric = false + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.23744, atol=1e-3) + @test isapprox(CD, 0.00254, atol=1e-5) + @test isapprox(CDiff, 0.00243, atol=1e-5) + @test isapprox(Cm, -0.02165, atol=1e-4) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 3 - Wing at High Angle of Attack" begin + + # Simple Wing at High Angle of Attack + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + symmetric = true + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 8.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.80348, atol=1e-3) + @test isapprox(CD, 0.02651, atol=1e-4) + @test isapprox(CDiff, 0.02696, atol=1e-5) + @test isapprox(Cm, -0.07399, atol=1e-3) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 4 - Wing with Dihedral" begin + + # Simple Wing with Dihedral + + # NOTE: There is some interaction between twist, dihedral, and chordwise + # position which causes the normal vectors found by AVL to differ from those + # computed by this package. We therefore manually overwrite the normal + # vectors when this occurs in order to get a better comparison. + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 3.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + symmetric = true + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + # also get normal vector as AVL defines it + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + for (ip, p) in enumerate(surface) + # check that our normal vector is approximately the same as AVL's + @test isapprox(p.ncp, ncp, rtol=0.01) + # replace our normal vector with AVL's normal vector for this test + surface[ip] = set_normal(p, ncp) + end + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.24787, atol=1e-3) + @test isapprox(CD, 0.00246, atol=1e-4) + @test isapprox(CDiff, 0.00245, atol=1e-5) + @test isapprox(Cm, -0.02395, atol=1e-4) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 5 - Wing with Dihedral at Very High Angle of Attack" begin + + # Simple Wing with Dihedral at Very High Angle of Attack + + # NOTE: this test case is nonphysical, so it just tests the numerics + + # NOTE: There is some interaction between twist, dihedral, and chordwise + # position which causes the normal vectors found by AVL to differ from those + # computed by this package. We therefore manually overwrite the normal + # vectors when this occurs in order to get a better comparison. + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 3.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + symmetric = true + + # adjust chord length to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 20.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + # vortex rings, untwisted geometry + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + for (ip, p) in enumerate(surface) + # check that our normal vector is approximately the same as AVL's + @test isapprox(p.ncp, ncp, rtol=0.01) + # replace our normal vector with AVL's normal vector for this test + surface[ip] = set_normal(p, ncp) + end + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 1.70982, rtol=0.02) + @test isapprox(CD, 0.12904, rtol=0.02) + @test isapprox(CDiff, 0.11502, rtol=0.02) + @test isapprox(Cm, -0.45606, rtol=0.02) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 6 - Wing and Tail without Finite Core Model" begin + + # Wing and Tail without Finite Core Model + + # NOTE: AVL's finite-core model is turned off for these tests + + # NOTE: There is some interaction between twist, dihedral, and chordwise + # position which causes the normal vectors found by AVL to differ from those + # computed by this package. We therefore manually overwrite the normal + # vectors when this occurs in order to get a better comparison. + + # wing + xle = [0.0, 0.2] + yle = [0.0, 5.0] + zle = [0.0, 1.0] + chord = [1.0, 0.6] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + + # horizontal stabilizer + xle_h = [0.0, 0.14] + yle_h = [0.0, 1.25] + zle_h = [0.0, 0.0] + chord_h = [0.7, 0.42] + theta_h = [0.0, 0.0] + phi_h = [0.0, 0.0] + ns_h = 6 + nc_h = 1 + spacing_s_h = Uniform() + spacing_c_h = Uniform() + mirror_h = false + + # vertical stabilizer + xle_v = [0.0, 0.14] + yle_v = [0.0, 0.0] + zle_v = [0.0, 1.0] + chord_v = [0.7, 0.42] + theta_v = [0.0, 0.0] + phi_v = [0.0, 0.0] + ns_v = 5 + nc_v = 1 + spacing_s_v = Uniform() + spacing_c_v = Uniform() + mirror_v = false + + # adjust chord lengths to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + chord_h = @. chord_h/cos(theta_h) + chord_v = @. chord_v/cos(theta_v) + + Sref = 9.0 + cref = 0.9 + bref = 10.0 + rref = [0.5, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 5.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + symmetric = [true, true, false] + + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + # vortex rings - finite core deactivated + wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + for (ip, p) in enumerate(wing) + # check that our normal vector is approximately the same as AVL's + @test isapprox(p.ncp, ncp, rtol=0.01) + # replace our normal vector with AVL's normal vector for this test + wing[ip] = set_normal(p, ncp) + end + + hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; + mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) + translate!(htail, [4.0, 0.0, 0.0]) + + vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; + mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) + translate!(vtail, [4.0, 0.0, 0.0]) + + surfaces = [wing, htail, vtail] + surface_id = [1, 1, 1] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.60408, atol=1e-3) + @test isapprox(CD, 0.01058, atol=1e-4) + @test isapprox(CDiff, 0.010378, atol=1e-4) + @test isapprox(Cm, -0.02778, atol=2e-3) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 7 - Wing and Tail with Finite Core Model" begin + + # Wing and Tail with Finite Core Model + + # NOTE: There is some interaction between twist, dihedral, and chordwise + # position which causes the normal vectors found by AVL to differ from those + # computed by this package. We therefore manually overwrite the normal + # vectors when this occurs in order to get a better comparison. + + # wing + xle = [0.0, 0.2] + yle = [0.0, 5.0] + zle = [0.0, 1.0] + chord = [1.0, 0.6] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + + # horizontal stabilizer + xle_h = [0.0, 0.14] + yle_h = [0.0, 1.25] + zle_h = [0.0, 0.0] + chord_h = [0.7, 0.42] + theta_h = [0.0, 0.0] + phi_h = [0.0, 0.0] + ns_h = 6 + nc_h = 1 + spacing_s_h = Uniform() + spacing_c_h = Uniform() + mirror_h = false + + # vertical stabilizer + xle_v = [0.0, 0.14] + yle_v = [0.0, 0.0] + zle_v = [0.0, 1.0] + chord_v = [0.7, 0.42] + theta_v = [0.0, 0.0] + phi_v = [0.0, 0.0] + ns_v = 5 + nc_v = 1 + spacing_s_v = Uniform() + spacing_c_v = Uniform() + mirror_v = false + + # adjust chord lengths to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + chord_h = @. chord_h/cos(theta_h) + chord_v = @. chord_v/cos(theta_v) + + Sref = 9.0 + cref = 0.9 + bref = 10.0 + rref = [0.5, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 5.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + symmetric = [true, true, false] + + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c, fcore = (c,Δs)->max(c/4, Δs/2)) + + for (ip, p) in enumerate(wing) + # check that our normal vector is approximately the same as AVL's + @test isapprox(p.ncp, ncp, rtol=0.01) + # replace our normal vector with AVL's normal vector for this test + wing[ip] = set_normal(p, ncp) + end + + hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; + mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) + translate!(htail, [4.0, 0.0, 0.0]) + + vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; + mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) + translate!(vtail, [4.0, 0.0, 0.0]) + + surfaces = [wing, htail, vtail] + surface_id = [1, 2, 3] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, derivatives = false) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.60562, atol=1e-3) + @test isapprox(CD, 0.01058, atol=1e-4) + @test isapprox(CDiff, 0.0104855, atol=1e-4) + @test isapprox(Cm, -0.03377, atol=2e-3) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "AVL - Run 8 - Wing with Chordwise Panels" begin + + # Simple Wing with Chordwise Panels + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 6 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + + # adjust chord length to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + symmetric = true + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.24454, atol=1e-3) + @test isapprox(CD, 0.00247, atol=1e-5) + @test isapprox(CDiff, 0.00248, atol=1e-5) + @test isapprox(Cm, -0.02091, atol=1e-4) + @test isapprox(CY, 0.0, atol=1e-16) + @test isapprox(Cl, 0.0, atol=1e-16) + @test isapprox(Cn, 0.0, atol=1e-16) +end + +@testset "AVL - Run 9 - Wing with Cosine-Spaced Spanwise and Chordwise Panels" begin + + # Simple Wing with Cosine-Spaced Spanwise and Chordwise Panels + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 6 + spacing_s = Cosine() + spacing_c = Cosine() + mirror = false + + # adjust chord length to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + symmetric = true + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.23879, atol=1e-3) + @test isapprox(CD, 0.00249, atol=1e-5) + @test isapprox(CDiff, 0.0024626, atol=1e-5) + @test isapprox(Cm, -0.01995, atol=1e-4) + @test isapprox(CY, 0.0, atol=1e-16) + @test isapprox(Cl, 0.0, atol=1e-16) + @test isapprox(Cn, 0.0, atol=1e-16) +end + +@testset "AVL - Run 10 - Wing with Sideslip" begin + + # Simple Wing with Sideslip + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + + # adjust chord length to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 15.0*pi/180 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + # vortex rings + mirror = true + symmetric = false + + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.22695, atol=1e-3) + @test isapprox(CD, 0.00227, atol=1e-5) + @test isapprox(CDiff, 0.0022852, atol=1e-5) + @test isapprox(Cm, -0.02101, atol=1e-4) + @test isapprox(CY, 0.0, atol=1e-5) + @test isapprox(Cl, -0.00644, atol=1e-4) + @test isapprox(Cn, 0.00012, atol=2e-4) +end + +@testset "AVL - Run 11 - Wing Stability Derivatives" begin + + # Run 11: Simple Wing Stability Derivatives + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = true + symmetric = false + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + dCF, dCM = stability_derivatives(system) + + CDa, CYa, CLa = dCF.alpha + Cla, Cma, Cna = dCM.alpha + CDb, CYb, CLb = dCF.beta + Clb, Cmb, Cnb = dCM.beta + CDp, CYp, CLp = dCF.p + Clp, Cmp, Cnp = dCM.p + CDq, CYq, CLq = dCF.q + Clq, Cmq, Cnq = dCM.q + CDr, CYr, CLr = dCF.r + Clr, Cmr, Cnr = dCM.r + + @test isapprox(CLa, 4.638088, rtol=0.01) + @test isapprox(CLb, 0.0, atol=ztol) + @test isapprox(CYa, 0.0, atol=ztol) + @test isapprox(CYb, -0.000007, atol=1e-4) + @test isapprox(Cla, 0.0, atol=ztol) + @test isapprox(Clb, -0.025749, atol=0.001) + @test isapprox(Cma, -0.429247, rtol=0.01) + @test isapprox(Cmb, 0.0, atol=ztol) + @test isapprox(Cna, 0.0, atol=ztol) + @test isapprox(Cnb, 0.000466, atol=1e-3) + @test isapprox(Clp, -0.518725, rtol=0.01) + @test isapprox(Clq, 0.0, atol=ztol) + @test isapprox(Clr, 0.064243, rtol=0.01) + @test isapprox(Cmp, 0.0, atol=ztol) + @test isapprox(Cmq, -0.517094, rtol=0.01) + @test isapprox(Cmr, 0.0, atol=ztol) + @test isapprox(Cnp, -0.019846, rtol=0.01) + @test isapprox(Cnq, 0.0, atol=ztol) + @test isapprox(Cnr, -0.000898, rtol=0.01) +end + +@testset "AVL - Run 12 - Rotational Velocity" begin + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = true + symmetric = false + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [2*Vinf*0.05/bref; 0.0; 0.0] # nondimensional pbar = 0.05 + fs = Freestream(Vinf, alpha, beta, Omega) + + # vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) + CF, CM = body_forces(system; frame=Stability()) -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) + CD, CY, CL = CF + Cl, Cm, Cn = CM -# # vortex rings with symmetry -# mirror = false -# symmetric = true + @test isapprox(CL, 0.24323, atol=1e-3) + @test isapprox(CD, 0.00069, atol=1e-5) + @test isapprox(Cm, -0.02251, atol=1e-4) + @test isapprox(CY, 0.00235, atol=2e-4) + @test isapprox(Cl, -0.02594, atol=1e-4) + @test isapprox(Cn, -0.00099, atol=2e-4) -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) +end + +@testset "Lifting Line Coefficients" begin + # Simple Wing with Uniform Spacing + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + # vortex rings with symmetry + mirror = false + symmetric = true + + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + grids = [grid] + surfaces = [surface] + + system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + + r_ll, c_ll = lifting_line_geometry(grids) + + cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability()) + + cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513, + 0.2404, 0.2233, 0.1952, 0.1434] + cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024, + 0.0025, 0.0026, 0.0026, 0.0022] + cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + + @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf)) + @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf)) + @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf)) +end + +@testset "Geometry Generation" begin + + # Tests of the geometry generation functions + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 6 + spacing_s = Uniform() + spacing_c = Uniform() + + # Symmetric Geometry + + halfgrid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + spacing_s=spacing_s, spacing_c=spacing_c) + + halfgrid2, surface2 = grid_to_surface_panels(halfgrid1) + + halfgrid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; + spacing_s=spacing_s, spacing_c=spacing_c) + + surface4 = similar(surface1) + VortexLattice.update_surface_panels!(surface4, halfgrid1) + + @test isapprox(halfgrid1, halfgrid2) + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) + end + end + + @test isapprox(halfgrid1, halfgrid3) + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) + end + end + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) + end + end + + # Mirrored Geometry + + mirror = true + + grid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + grid2, surface2 = grid_to_surface_panels(halfgrid1; mirror = mirror) + + grid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; mirror=mirror, + spacing_s=spacing_s, spacing_c=spacing_c) + + surface4 = similar(surface1) + VortexLattice.update_surface_panels!(surface4, grid1) + + @test isapprox(grid1, grid2) + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) + end + end + + @test isapprox(grid1, grid3) + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) + end + end + + for I in CartesianIndices(surface1) + for field in fieldnames(SurfacePanel) + @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) + end + end +end + +@testset "Grid Input" begin + + # Tests for inputting a grid rather than a surface + + xle = [0.0, 0.4] + yle = [0.0, 7.5] + zle = [0.0, 0.0] + chord = [2.2, 1.8] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 6 + spacing_s = Uniform() + spacing_c = Uniform() + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + # adjust chord length so x-chord length matches AVL + chord = @. chord/cos(theta) + + mirror = false + symmetric = true + + halfgrid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + # steady analysis, single grid + + grids = [halfgrid] + + system = steady_analysis(grids, ref, fs; symmetric=symmetric) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + @test isapprox(CL, 0.24454, atol=1e-3) + @test isapprox(CD, 0.00247, atol=1e-5) + @test isapprox(CDiff, 0.00248, atol=1e-5) + @test isapprox(Cm, -0.02091, atol=1e-4) + @test isapprox(CY, 0.0, atol=1e-16) + @test isapprox(Cl, 0.0, atol=1e-16) + @test isapprox(Cn, 0.0, atol=1e-16) + + # steady analysis multiple grids + + # NOTE: AVL's finite-core model is turned off for these tests + + # NOTE: There is some interaction between twist, dihedral, and chordwise + # position which causes the normal vectors found by AVL to differ from those + # computed by this package. We therefore manually overwrite the normal + # vectors when this occurs in order to get a better comparison. + + # wing + xle = [0.0, 0.2] + yle = [0.0, 5.0] + zle = [0.0, 1.0] + chord = [1.0, 0.6] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + + # horizontal stabilizer + xle_h = [0.0, 0.14] + yle_h = [0.0, 1.25] + zle_h = [0.0, 0.0] + chord_h = [0.7, 0.42] + theta_h = [0.0, 0.0] + phi_h = [0.0, 0.0] + ns_h = 6 + nc_h = 1 + spacing_s_h = Uniform() + spacing_c_h = Uniform() + mirror_h = false + + # vertical stabilizer + xle_v = [0.0, 0.14] + yle_v = [0.0, 0.0] + zle_v = [0.0, 1.0] + chord_v = [0.7, 0.42] + theta_v = [0.0, 0.0] + phi_v = [0.0, 0.0] + ns_v = 5 + nc_v = 1 + spacing_s_v = Uniform() + spacing_c_v = Uniform() + mirror_v = false + + # adjust chord lengths to match AVL (which uses chord length in the x-direction) + chord = @. chord/cos(theta) + chord_h = @. chord_h/cos(theta_h) + chord_v = @. chord_v/cos(theta_v) + + Sref = 9.0 + cref = 0.9 + bref = 10.0 + rref = [0.5, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) + + alpha = 5.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) + + symmetric = [true, true, false] + + ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) + + # vortex rings - finite core deactivated + wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + + hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; + mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) + translate!(hgrid, [4.0, 0.0, 0.0]) + translate!(htail, [4.0, 0.0, 0.0]) + + vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; + mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) + translate!(vgrid, [4.0, 0.0, 0.0]) + translate!(vtail, [4.0, 0.0, 0.0]) + + grids = [wgrid, hgrid, vgrid] + surface_id = [1, 1, 1] + + system = steady_analysis(grids, ref, fs; symmetric=symmetric, surface_id=surface_id) + + CF, CM = body_forces(system; frame=Stability()) + + CDiff = far_field_drag(system) + + CD, CY, CL = CF + Cl, Cm, Cn = CM + + # some differences are expected since we don't set the normal vector in + # VortexLattice equal to that used by AVL + + @test isapprox(CL, 0.60408, rtol=0.01) + @test isapprox(CD, 0.01058, rtol=0.01) + @test isapprox(CDiff, 0.010378, atol=0.02) + @test isapprox(Cm, -0.02778, rtol=0.01) + @test isapprox(CY, 0.0, atol=ztol) + @test isapprox(Cl, 0.0, atol=ztol) + @test isapprox(Cn, 0.0, atol=ztol) +end + +@testset "Update Trailing Edge Coefficients" begin + + # This test checks whether the function which updates the trailing edge + # coefficients `update_trailing_edge_coefficients!` results in the same + # trailing edge coefficients as `influnece_coefficients!` + + # wing + xle = [0.0, 0.2] + yle = [0.0, 5.0] + zle = [0.0, 1.0] + chord = [1.0, 0.6] + theta = [2.0*pi/180, 2.0*pi/180] + phi = [0.0, 0.0] + ns = 12 + nc = 1 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + + # horizontal stabilizer + xle_h = [0.0, 0.14] + yle_h = [0.0, 1.25] + zle_h = [0.0, 0.0] + chord_h = [0.7, 0.42] + theta_h = [0.0, 0.0] + phi_h = [0.0, 0.0] + ns_h = 6 + nc_h = 1 + spacing_s_h = Uniform() + spacing_c_h = Uniform() + mirror_h = false + + # vertical stabilizer + xle_v = [0.0, 0.14] + yle_v = [0.0, 0.0] + zle_v = [0.0, 1.0] + chord_v = [0.7, 0.42] + theta_v = [0.0, 0.0] + phi_v = [0.0, 0.0] + ns_v = 5 + nc_v = 1 + spacing_s_v = Uniform() + spacing_c_v = Uniform() + mirror_v = false -# surfaces = [surface] + Sref = 9.0 + cref = 0.9 + bref = 10.0 + rref = [0.5, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) + alpha = 5.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) -# CF, CM = body_forces(system; frame=Stability()) + symmetric = [true, true, false] -# CDiff = far_field_drag(system) + # horseshoe vortices + wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) -# CD, CY, CL = CF -# Cl, Cm, Cn = CM + hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; + mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) + translate!(htail, [4.0, 0.0, 0.0]) -# @test isapprox(CL, 0.24324, atol=1e-3) -# @test isapprox(CD, 0.00243, atol=1e-5) -# @test isapprox(CDiff, 0.00245, atol=1e-5) -# @test isapprox(Cm, -0.02252, atol=1e-4) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) - -# # vortex rings with mirrored geometry -# mirror = true -# symmetric = false - -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.24324, atol=1e-3) -# @test isapprox(CD, 0.00243, atol=1e-5) -# @test isapprox(CDiff, 0.00245, atol=1e-5) -# @test isapprox(Cm, -0.02252, atol=1e-4) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end + vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; + mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) + translate!(vtail, [4.0, 0.0, 0.0]) -# @testset "AVL - Run 2 - Wing with Cosine Spacing" begin + surfaces = [wing, htail, vtail] + surface_id = [1, 2, 2] -# # Run 2: Simple Wing with Cosine Spacing + # number of panels + N = nc*ns + nc_h*ns_h + nc_v*ns_v + AIC1 = zeros(N, N) -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Cosine() -# spacing_c = Uniform() -# mirror = true -# symmetric = false + VortexLattice.influence_coefficients!(AIC1, surfaces; + symmetric = symmetric, + trailing_vortices = fill(false, length(surfaces)), + surface_id = surface_id) -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.23744, atol=1e-3) -# @test isapprox(CD, 0.00254, atol=1e-5) -# @test isapprox(CDiff, 0.00243, atol=1e-5) -# @test isapprox(Cm, -0.02165, atol=1e-4) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 3 - Wing at High Angle of Attack" begin - -# # Simple Wing at High Angle of Attack - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false -# symmetric = true - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 8.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.80348, atol=1e-3) -# @test isapprox(CD, 0.02651, atol=1e-4) -# @test isapprox(CDiff, 0.02696, atol=1e-5) -# @test isapprox(Cm, -0.07399, atol=1e-3) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 4 - Wing with Dihedral" begin - -# # Simple Wing with Dihedral - -# # NOTE: There is some interaction between twist, dihedral, and chordwise -# # position which causes the normal vectors found by AVL to differ from those -# # computed by this package. We therefore manually overwrite the normal -# # vectors when this occurs in order to get a better comparison. - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 3.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false -# symmetric = true - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) - -# # also get normal vector as AVL defines it -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# for (ip, p) in enumerate(surface) -# # check that our normal vector is approximately the same as AVL's -# @test isapprox(p.ncp, ncp, rtol=0.01) -# # replace our normal vector with AVL's normal vector for this test -# surface[ip] = set_normal(p, ncp) -# end - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.24787, atol=1e-3) -# @test isapprox(CD, 0.00246, atol=1e-4) -# @test isapprox(CDiff, 0.00245, atol=1e-5) -# @test isapprox(Cm, -0.02395, atol=1e-4) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 5 - Wing with Dihedral at Very High Angle of Attack" begin - -# # Simple Wing with Dihedral at Very High Angle of Attack - -# # NOTE: this test case is nonphysical, so it just tests the numerics - -# # NOTE: There is some interaction between twist, dihedral, and chordwise -# # position which causes the normal vectors found by AVL to differ from those -# # computed by this package. We therefore manually overwrite the normal -# # vectors when this occurs in order to get a better comparison. - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 3.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false -# symmetric = true - -# # adjust chord length to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 20.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# # vortex rings, untwisted geometry -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# for (ip, p) in enumerate(surface) -# # check that our normal vector is approximately the same as AVL's -# @test isapprox(p.ncp, ncp, rtol=0.01) -# # replace our normal vector with AVL's normal vector for this test -# surface[ip] = set_normal(p, ncp) -# end - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 1.70982, rtol=0.02) -# @test isapprox(CD, 0.12904, rtol=0.02) -# @test isapprox(CDiff, 0.11502, rtol=0.02) -# @test isapprox(Cm, -0.45606, rtol=0.02) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 6 - Wing and Tail without Finite Core Model" begin - -# # Wing and Tail without Finite Core Model - -# # NOTE: AVL's finite-core model is turned off for these tests - -# # NOTE: There is some interaction between twist, dihedral, and chordwise -# # position which causes the normal vectors found by AVL to differ from those -# # computed by this package. We therefore manually overwrite the normal -# # vectors when this occurs in order to get a better comparison. - -# # wing -# xle = [0.0, 0.2] -# yle = [0.0, 5.0] -# zle = [0.0, 1.0] -# chord = [1.0, 0.6] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false - -# # horizontal stabilizer -# xle_h = [0.0, 0.14] -# yle_h = [0.0, 1.25] -# zle_h = [0.0, 0.0] -# chord_h = [0.7, 0.42] -# theta_h = [0.0, 0.0] -# phi_h = [0.0, 0.0] -# ns_h = 6 -# nc_h = 1 -# spacing_s_h = Uniform() -# spacing_c_h = Uniform() -# mirror_h = false - -# # vertical stabilizer -# xle_v = [0.0, 0.14] -# yle_v = [0.0, 0.0] -# zle_v = [0.0, 1.0] -# chord_v = [0.7, 0.42] -# theta_v = [0.0, 0.0] -# phi_v = [0.0, 0.0] -# ns_v = 5 -# nc_v = 1 -# spacing_s_v = Uniform() -# spacing_c_v = Uniform() -# mirror_v = false - -# # adjust chord lengths to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) -# chord_h = @. chord_h/cos(theta_h) -# chord_v = @. chord_v/cos(theta_v) - -# Sref = 9.0 -# cref = 0.9 -# bref = 10.0 -# rref = [0.5, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 5.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = [true, true, false] - -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# # vortex rings - finite core deactivated -# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# for (ip, p) in enumerate(wing) -# # check that our normal vector is approximately the same as AVL's -# @test isapprox(p.ncp, ncp, rtol=0.01) -# # replace our normal vector with AVL's normal vector for this test -# wing[ip] = set_normal(p, ncp) -# end - -# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; -# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) -# translate!(htail, [4.0, 0.0, 0.0]) - -# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; -# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) -# translate!(vtail, [4.0, 0.0, 0.0]) - -# surfaces = [wing, htail, vtail] -# surface_id = [1, 1, 1] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.60408, atol=1e-3) -# @test isapprox(CD, 0.01058, atol=1e-4) -# @test isapprox(CDiff, 0.010378, atol=1e-4) -# @test isapprox(Cm, -0.02778, atol=2e-3) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 7 - Wing and Tail with Finite Core Model" begin - -# # Wing and Tail with Finite Core Model - -# # NOTE: There is some interaction between twist, dihedral, and chordwise -# # position which causes the normal vectors found by AVL to differ from those -# # computed by this package. We therefore manually overwrite the normal -# # vectors when this occurs in order to get a better comparison. - -# # wing -# xle = [0.0, 0.2] -# yle = [0.0, 5.0] -# zle = [0.0, 1.0] -# chord = [1.0, 0.6] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false - -# # horizontal stabilizer -# xle_h = [0.0, 0.14] -# yle_h = [0.0, 1.25] -# zle_h = [0.0, 0.0] -# chord_h = [0.7, 0.42] -# theta_h = [0.0, 0.0] -# phi_h = [0.0, 0.0] -# ns_h = 6 -# nc_h = 1 -# spacing_s_h = Uniform() -# spacing_c_h = Uniform() -# mirror_h = false - -# # vertical stabilizer -# xle_v = [0.0, 0.14] -# yle_v = [0.0, 0.0] -# zle_v = [0.0, 1.0] -# chord_v = [0.7, 0.42] -# theta_v = [0.0, 0.0] -# phi_v = [0.0, 0.0] -# ns_v = 5 -# nc_v = 1 -# spacing_s_v = Uniform() -# spacing_c_v = Uniform() -# mirror_v = false - -# # adjust chord lengths to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) -# chord_h = @. chord_h/cos(theta_h) -# chord_v = @. chord_v/cos(theta_v) - -# Sref = 9.0 -# cref = 0.9 -# bref = 10.0 -# rref = [0.5, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 5.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = [true, true, false] - -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c, fcore = (c,Δs)->max(c/4, Δs/2)) - -# for (ip, p) in enumerate(wing) -# # check that our normal vector is approximately the same as AVL's -# @test isapprox(p.ncp, ncp, rtol=0.01) -# # replace our normal vector with AVL's normal vector for this test -# wing[ip] = set_normal(p, ncp) -# end - -# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; -# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) -# translate!(htail, [4.0, 0.0, 0.0]) - -# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; -# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) -# translate!(vtail, [4.0, 0.0, 0.0]) - -# surfaces = [wing, htail, vtail] -# surface_id = [1, 2, 3] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, derivatives = false) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.60562, atol=1e-3) -# @test isapprox(CD, 0.01058, atol=1e-4) -# @test isapprox(CDiff, 0.0104855, atol=1e-4) -# @test isapprox(Cm, -0.03377, atol=2e-3) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "AVL - Run 8 - Wing with Chordwise Panels" begin - -# # Simple Wing with Chordwise Panels - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 6 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false - -# # adjust chord length to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = true - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.24454, atol=1e-3) -# @test isapprox(CD, 0.00247, atol=1e-5) -# @test isapprox(CDiff, 0.00248, atol=1e-5) -# @test isapprox(Cm, -0.02091, atol=1e-4) -# @test isapprox(CY, 0.0, atol=1e-16) -# @test isapprox(Cl, 0.0, atol=1e-16) -# @test isapprox(Cn, 0.0, atol=1e-16) -# end + AIC2 = copy(AIC1) -# @testset "AVL - Run 9 - Wing with Cosine-Spaced Spanwise and Chordwise Panels" begin + VortexLattice.update_trailing_edge_coefficients!(AIC2, surfaces; + symmetric = symmetric, + trailing_vortices = fill(false, length(surfaces)), + surface_id = surface_id) -# # Simple Wing with Cosine-Spaced Spanwise and Chordwise Panels + @test isapprox(AIC1, AIC2) +end -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 6 -# spacing_s = Cosine() -# spacing_c = Cosine() -# mirror = false +@testset "Wake Induced Velocity" begin + + # This test constructs two surfaces, one using surface panels, + # and one using wake panels. It then tests that the wake panel + # implementations of `induced_velocity` yields identical results to the + # surface panel implementations + + # generate the surface/wake geometry + nc = 5 + ns = 10 + + x = range(0, 1, length = nc+1) + y = range(0, 2, length = ns+1) + z = range(0, 3, length = ns+1) + + surface = Matrix{SurfacePanel{Float64}}(undef, nc, ns) + wake = Matrix{WakePanel{Float64}}(undef, nc, ns) + Γ = rand(nc, ns) + + for i = 1:nc, j = 1:ns + rtl = [x[i], y[j], z[j]] + rtr = [x[i], y[j+1], z[j+1]] + rbl = [x[i+1], y[j], z[j]] + rbr = [x[i+1], y[j+1], z[j+1]] + rcp = (rtl + rtr + rbl + rbr)/4 + ncp = cross(rcp - rtr, rcp - rtl) + core_size = 0.1 + chord = 0.0 # only used for unsteady simuulations + + surface[i,j] = SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord) + wake[i,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, Γ[i, j]) + end -# # adjust chord length to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) + # Test induced velocity calculation at an arbitrary point in space: + rcp = [10, 11, 12] -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) + # no finite core model, no symmetry, no trailing vortices -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = true - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.23879, atol=1e-3) -# @test isapprox(CD, 0.00249, atol=1e-5) -# @test isapprox(CDiff, 0.0024626, atol=1e-5) -# @test isapprox(Cm, -0.01995, atol=1e-4) -# @test isapprox(CY, 0.0, atol=1e-16) -# @test isapprox(Cl, 0.0, atol=1e-16) -# @test isapprox(Cn, 0.0, atol=1e-16) -# end + Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; + finite_core = false, + symmetric = false, + trailing_vortices = false, + xhat = [1, 0, 0]) -# @testset "AVL - Run 10 - Wing with Sideslip" begin + Vw = VortexLattice.induced_velocity(rcp, wake; + finite_core = false, + symmetric = false, + trailing_vortices = false, + xhat = [1, 0, 0]) -# # Simple Wing with Sideslip - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() - -# # adjust chord length to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 15.0*pi/180 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# # vortex rings -# mirror = true -# symmetric = false - -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.22695, atol=1e-3) -# @test isapprox(CD, 0.00227, atol=1e-5) -# @test isapprox(CDiff, 0.0022852, atol=1e-5) -# @test isapprox(Cm, -0.02101, atol=1e-4) -# @test isapprox(CY, 0.0, atol=1e-5) -# @test isapprox(Cl, -0.00644, atol=1e-4) -# @test isapprox(Cn, 0.00012, atol=2e-4) -# end - -# @testset "AVL - Run 11 - Wing Stability Derivatives" begin - -# # Run 11: Simple Wing Stability Derivatives - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = true -# symmetric = false - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# dCF, dCM = stability_derivatives(system) - -# CDa, CYa, CLa = dCF.alpha -# Cla, Cma, Cna = dCM.alpha -# CDb, CYb, CLb = dCF.beta -# Clb, Cmb, Cnb = dCM.beta -# CDp, CYp, CLp = dCF.p -# Clp, Cmp, Cnp = dCM.p -# CDq, CYq, CLq = dCF.q -# Clq, Cmq, Cnq = dCM.q -# CDr, CYr, CLr = dCF.r -# Clr, Cmr, Cnr = dCM.r - -# @test isapprox(CLa, 4.638088, rtol=0.01) -# @test isapprox(CLb, 0.0, atol=ztol) -# @test isapprox(CYa, 0.0, atol=ztol) -# @test isapprox(CYb, -0.000007, atol=1e-4) -# @test isapprox(Cla, 0.0, atol=ztol) -# @test isapprox(Clb, -0.025749, atol=0.001) -# @test isapprox(Cma, -0.429247, rtol=0.01) -# @test isapprox(Cmb, 0.0, atol=ztol) -# @test isapprox(Cna, 0.0, atol=ztol) -# @test isapprox(Cnb, 0.000466, atol=1e-3) -# @test isapprox(Clp, -0.518725, rtol=0.01) -# @test isapprox(Clq, 0.0, atol=ztol) -# @test isapprox(Clr, 0.064243, rtol=0.01) -# @test isapprox(Cmp, 0.0, atol=ztol) -# @test isapprox(Cmq, -0.517094, rtol=0.01) -# @test isapprox(Cmr, 0.0, atol=ztol) -# @test isapprox(Cnp, -0.019846, rtol=0.01) -# @test isapprox(Cnq, 0.0, atol=ztol) -# @test isapprox(Cnr, -0.000898, rtol=0.01) -# end - -# @testset "AVL - Run 12 - Rotational Velocity" begin - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = true -# symmetric = false - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [2*Vinf*0.05/bref; 0.0; 0.0] # nondimensional pbar = 0.05 -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.24323, atol=1e-3) -# @test isapprox(CD, 0.00069, atol=1e-5) -# @test isapprox(Cm, -0.02251, atol=1e-4) -# @test isapprox(CY, 0.00235, atol=2e-4) -# @test isapprox(Cl, -0.02594, atol=1e-4) -# @test isapprox(Cn, -0.00099, atol=2e-4) - -# end - -# @testset "Lifting Line Coefficients" begin -# # Simple Wing with Uniform Spacing - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) - -# # vortex rings with symmetry -# mirror = false -# symmetric = true - -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# grids = [grid] -# surfaces = [surface] - -# system = steady_analysis(surfaces, ref, fs; symmetric=symmetric) - -# r_ll, c_ll = lifting_line_geometry(grids) - -# cf, cm = lifting_line_coefficients(system, r_ll, c_ll; frame=Stability()) - -# cl_avl = [0.2618, 0.2646, 0.2661, 0.2664, 0.2654, 0.2628, 0.2584, 0.2513, -# 0.2404, 0.2233, 0.1952, 0.1434] -# cd_avl = [0.0029, 0.0024, 0.0023, 0.0023, 0.0023, 0.0023, 0.0024, 0.0024, -# 0.0025, 0.0026, 0.0026, 0.0022] -# cm_avl = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] + @test isapprox(Vs, Vw) -# @test isapprox(cf[1][3,:], cl_avl, atol=1e-3, norm=(x)->norm(x, Inf)) -# @test isapprox(cf[1][1,:], cd_avl, atol=1e-4, norm=(x)->norm(x, Inf)) -# @test isapprox(cm[1][2,:], cm_avl, atol=1e-4, norm=(x)->norm(x, Inf)) -# end - -# @testset "Geometry Generation" begin - -# # Tests of the geometry generation functions + # finite core model, symmetry, and trailing vortices -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 6 -# spacing_s = Uniform() -# spacing_c = Uniform() + Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; + finite_core = true, + symmetric = true, + trailing_vortices = true, + xhat = [1, 0, 0]) -# # Symmetric Geometry + Vw = VortexLattice.induced_velocity(rcp, wake; + finite_core = true, + symmetric = true, + trailing_vortices = true, + xhat = [1, 0, 0]) -# halfgrid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# spacing_s=spacing_s, spacing_c=spacing_c) + @test isapprox(Vs, Vw) -# halfgrid2, surface2 = grid_to_surface_panels(halfgrid1) + # Test induced velocity calculation at a trailing edge point: + is = 3 # trailing edge index + I = CartesianIndex(nc+1, is) # index on surface -# halfgrid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; -# spacing_s=spacing_s, spacing_c=spacing_c) + # no finite core model, no symmetry, no trailing vortices + Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; + finite_core = false, + symmetric = false, + trailing_vortices = false, + xhat = [1, 0, 0]) -# surface4 = similar(surface1) -# VortexLattice.update_surface_panels!(surface4, halfgrid1) + Vw = VortexLattice.induced_velocity(I, wake; + finite_core = false, + symmetric = false, + trailing_vortices = false, + xhat = [1, 0, 0]) -# @test isapprox(halfgrid1, halfgrid2) + @test isapprox(Vs, Vw) -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) -# end -# end + # finite core model, symmetry, and trailing vortices + Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; + finite_core = true, + symmetric = true, + trailing_vortices = true, + xhat = [1, 0, 0]) -# @test isapprox(halfgrid1, halfgrid3) + Vw = VortexLattice.induced_velocity(I, wake; + finite_core = true, + symmetric = true, + trailing_vortices = true, + xhat = [1, 0, 0]) -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) -# end -# end + @test isapprox(Vs, Vw) -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) -# end -# end +end -# # Mirrored Geometry +@testset "Unsteady Vortex Lattice Method - Rectangular Wing" begin -# mirror = true + Uinf = 1.0 + AR = 4 -# grid1, surface1 = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) + # reference parameters + cref = 1.0 + bref = AR + Sref = bref*cref + rref = [0.0, 0.0, 0.0] + Vinf = 1.0 + ref = Reference(Sref, cref, bref, rref, Vinf) -# grid2, surface2 = grid_to_surface_panels(halfgrid1; mirror = mirror) + # freestream parameters + alpha = 5.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = Freestream(Vinf, alpha, beta, Omega) -# grid3, surface3 = grid_to_surface_panels(halfgrid1, ns, nc; mirror=mirror, -# spacing_s=spacing_s, spacing_c=spacing_c) + # geometry + xle = [0.0, 0.0] + yle = [-bref/2, bref/2] + zle = [0.0, 0.0] + chord = [cref, cref] + theta = [0.0, 0.0] + phi = [0.0, 0.0] + ns = 13 + nc = 4 + spacing_s = Uniform() + spacing_c = Uniform() + mirror = false + symmetric = false -# surface4 = similar(surface1) -# VortexLattice.update_surface_panels!(surface4, grid1) + # non-dimensional time + t = range(0.0, 10.0, step=0.2) + dt = [t[i+1]-t[i] for i = 1:length(t)-1] -# @test isapprox(grid1, grid2) + # create vortex rings + grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; + mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface2[I], field)) -# end -# end + surfaces = [surface] -# @test isapprox(grid1, grid3) + # run analysis + system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; + symmetric=symmetric) -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface3[I], field)) -# end -# end + # extract forces at each time step + CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) +end -# for I in CartesianIndices(surface1) -# for field in fieldnames(SurfacePanel) -# @test isapprox(getproperty(surface1[I], field), getproperty(surface4[I], field)) -# end -# end -# end - -# @testset "Grid Input" begin - -# # Tests for inputting a grid rather than a surface - -# xle = [0.0, 0.4] -# yle = [0.0, 7.5] -# zle = [0.0, 0.0] -# chord = [2.2, 1.8] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 6 -# spacing_s = Uniform() -# spacing_c = Uniform() - -# Sref = 30.0 -# cref = 2.0 -# bref = 15.0 -# rref = [0.50, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 1.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # adjust chord length so x-chord length matches AVL -# chord = @. chord/cos(theta) - -# mirror = false -# symmetric = true - -# halfgrid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# # steady analysis, single grid - -# grids = [halfgrid] - -# system = steady_analysis(grids, ref, fs; symmetric=symmetric) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# @test isapprox(CL, 0.24454, atol=1e-3) -# @test isapprox(CD, 0.00247, atol=1e-5) -# @test isapprox(CDiff, 0.00248, atol=1e-5) -# @test isapprox(Cm, -0.02091, atol=1e-4) -# @test isapprox(CY, 0.0, atol=1e-16) -# @test isapprox(Cl, 0.0, atol=1e-16) -# @test isapprox(Cn, 0.0, atol=1e-16) - -# # steady analysis multiple grids - -# # NOTE: AVL's finite-core model is turned off for these tests - -# # NOTE: There is some interaction between twist, dihedral, and chordwise -# # position which causes the normal vectors found by AVL to differ from those -# # computed by this package. We therefore manually overwrite the normal -# # vectors when this occurs in order to get a better comparison. - -# # wing -# xle = [0.0, 0.2] -# yle = [0.0, 5.0] -# zle = [0.0, 1.0] -# chord = [1.0, 0.6] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false - -# # horizontal stabilizer -# xle_h = [0.0, 0.14] -# yle_h = [0.0, 1.25] -# zle_h = [0.0, 0.0] -# chord_h = [0.7, 0.42] -# theta_h = [0.0, 0.0] -# phi_h = [0.0, 0.0] -# ns_h = 6 -# nc_h = 1 -# spacing_s_h = Uniform() -# spacing_c_h = Uniform() -# mirror_h = false - -# # vertical stabilizer -# xle_v = [0.0, 0.14] -# yle_v = [0.0, 0.0] -# zle_v = [0.0, 1.0] -# chord_v = [0.7, 0.42] -# theta_v = [0.0, 0.0] -# phi_v = [0.0, 0.0] -# ns_v = 5 -# nc_v = 1 -# spacing_s_v = Uniform() -# spacing_c_v = Uniform() -# mirror_v = false - -# # adjust chord lengths to match AVL (which uses chord length in the x-direction) -# chord = @. chord/cos(theta) -# chord_h = @. chord_h/cos(theta_h) -# chord_v = @. chord_v/cos(theta_v) - -# Sref = 9.0 -# cref = 0.9 -# bref = 10.0 -# rref = [0.5, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 5.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = [true, true, false] - -# ncp = avl_normal_vector([xle[2]-xle[1], yle[2]-yle[1], zle[2]-zle[1]], 2.0*pi/180) - -# # vortex rings - finite core deactivated -# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; -# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) -# translate!(hgrid, [4.0, 0.0, 0.0]) -# translate!(htail, [4.0, 0.0, 0.0]) - -# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; -# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) -# translate!(vgrid, [4.0, 0.0, 0.0]) -# translate!(vtail, [4.0, 0.0, 0.0]) - -# grids = [wgrid, hgrid, vgrid] -# surface_id = [1, 1, 1] - -# system = steady_analysis(grids, ref, fs; symmetric=symmetric, surface_id=surface_id) - -# CF, CM = body_forces(system; frame=Stability()) - -# CDiff = far_field_drag(system) - -# CD, CY, CL = CF -# Cl, Cm, Cn = CM - -# # some differences are expected since we don't set the normal vector in -# # VortexLattice equal to that used by AVL - -# @test isapprox(CL, 0.60408, rtol=0.01) -# @test isapprox(CD, 0.01058, rtol=0.01) -# @test isapprox(CDiff, 0.010378, atol=0.02) -# @test isapprox(Cm, -0.02778, rtol=0.01) -# @test isapprox(CY, 0.0, atol=ztol) -# @test isapprox(Cl, 0.0, atol=ztol) -# @test isapprox(Cn, 0.0, atol=ztol) -# end - -# @testset "Update Trailing Edge Coefficients" begin - -# # This test checks whether the function which updates the trailing edge -# # coefficients `update_trailing_edge_coefficients!` results in the same -# # trailing edge coefficients as `influnece_coefficients!` - -# # wing -# xle = [0.0, 0.2] -# yle = [0.0, 5.0] -# zle = [0.0, 1.0] -# chord = [1.0, 0.6] -# theta = [2.0*pi/180, 2.0*pi/180] -# phi = [0.0, 0.0] -# ns = 12 -# nc = 1 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false - -# # horizontal stabilizer -# xle_h = [0.0, 0.14] -# yle_h = [0.0, 1.25] -# zle_h = [0.0, 0.0] -# chord_h = [0.7, 0.42] -# theta_h = [0.0, 0.0] -# phi_h = [0.0, 0.0] -# ns_h = 6 -# nc_h = 1 -# spacing_s_h = Uniform() -# spacing_c_h = Uniform() -# mirror_h = false - -# # vertical stabilizer -# xle_v = [0.0, 0.14] -# yle_v = [0.0, 0.0] -# zle_v = [0.0, 1.0] -# chord_v = [0.7, 0.42] -# theta_v = [0.0, 0.0] -# phi_v = [0.0, 0.0] -# ns_v = 5 -# nc_v = 1 -# spacing_s_v = Uniform() -# spacing_c_v = Uniform() -# mirror_v = false - -# Sref = 9.0 -# cref = 0.9 -# bref = 10.0 -# rref = [0.5, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# alpha = 5.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# symmetric = [true, true, false] - -# # horseshoe vortices -# wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h; -# mirror=mirror_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h) -# translate!(htail, [4.0, 0.0, 0.0]) - -# vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v; -# mirror=mirror_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v) -# translate!(vtail, [4.0, 0.0, 0.0]) - -# surfaces = [wing, htail, vtail] -# surface_id = [1, 2, 2] - -# # number of panels -# N = nc*ns + nc_h*ns_h + nc_v*ns_v -# AIC1 = zeros(N, N) - -# VortexLattice.influence_coefficients!(AIC1, surfaces; -# symmetric = symmetric, -# trailing_vortices = fill(false, length(surfaces)), -# surface_id = surface_id) - -# AIC2 = copy(AIC1) - -# VortexLattice.update_trailing_edge_coefficients!(AIC2, surfaces; -# symmetric = symmetric, -# trailing_vortices = fill(false, length(surfaces)), -# surface_id = surface_id) - -# @test isapprox(AIC1, AIC2) -# end - -# @testset "Wake Induced Velocity" begin - -# # This test constructs two surfaces, one using surface panels, -# # and one using wake panels. It then tests that the wake panel -# # implementations of `induced_velocity` yields identical results to the -# # surface panel implementations - -# # generate the surface/wake geometry -# nc = 5 -# ns = 10 - -# x = range(0, 1, length = nc+1) -# y = range(0, 2, length = ns+1) -# z = range(0, 3, length = ns+1) - -# surface = Matrix{SurfacePanel{Float64}}(undef, nc, ns) -# wake = Matrix{WakePanel{Float64}}(undef, nc, ns) -# Γ = rand(nc, ns) - -# for i = 1:nc, j = 1:ns -# rtl = [x[i], y[j], z[j]] -# rtr = [x[i], y[j+1], z[j+1]] -# rbl = [x[i+1], y[j], z[j]] -# rbr = [x[i+1], y[j+1], z[j+1]] -# rcp = (rtl + rtr + rbl + rbr)/4 -# ncp = cross(rcp - rtr, rcp - rtl) -# core_size = 0.1 -# chord = 0.0 # only used for unsteady simuulations - -# surface[i,j] = SurfacePanel(rtl, rtr, rbl, rbr, rcp, ncp, core_size, chord) -# wake[i,j] = WakePanel(rtl, rtr, rbl, rbr, core_size, Γ[i, j]) -# end - -# # Test induced velocity calculation at an arbitrary point in space: -# rcp = [10, 11, 12] - -# # no finite core model, no symmetry, no trailing vortices - -# Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; -# finite_core = false, -# symmetric = false, -# trailing_vortices = false, -# xhat = [1, 0, 0]) - -# Vw = VortexLattice.induced_velocity(rcp, wake; -# finite_core = false, -# symmetric = false, -# trailing_vortices = false, -# xhat = [1, 0, 0]) - -# @test isapprox(Vs, Vw) - -# # finite core model, symmetry, and trailing vortices - -# Vs = VortexLattice.induced_velocity(rcp, surface, Γ[:]; -# finite_core = true, -# symmetric = true, -# trailing_vortices = true, -# xhat = [1, 0, 0]) - -# Vw = VortexLattice.induced_velocity(rcp, wake; -# finite_core = true, -# symmetric = true, -# trailing_vortices = true, -# xhat = [1, 0, 0]) - -# @test isapprox(Vs, Vw) - -# # Test induced velocity calculation at a trailing edge point: -# is = 3 # trailing edge index -# I = CartesianIndex(nc+1, is) # index on surface - -# # no finite core model, no symmetry, no trailing vortices -# Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; -# finite_core = false, -# symmetric = false, -# trailing_vortices = false, -# xhat = [1, 0, 0]) - -# Vw = VortexLattice.induced_velocity(I, wake; -# finite_core = false, -# symmetric = false, -# trailing_vortices = false, -# xhat = [1, 0, 0]) - -# @test isapprox(Vs, Vw) - -# # finite core model, symmetry, and trailing vortices -# Vs = VortexLattice.induced_velocity(is, surface, Γ[:]; -# finite_core = true, -# symmetric = true, -# trailing_vortices = true, -# xhat = [1, 0, 0]) - -# Vw = VortexLattice.induced_velocity(I, wake; -# finite_core = true, -# symmetric = true, -# trailing_vortices = true, -# xhat = [1, 0, 0]) - -# @test isapprox(Vs, Vw) - -# end - -# @testset "Unsteady Vortex Lattice Method - Rectangular Wing" begin - -# Uinf = 1.0 -# AR = 4 - -# # reference parameters -# cref = 1.0 -# bref = AR -# Sref = bref*cref -# rref = [0.0, 0.0, 0.0] -# Vinf = 1.0 -# ref = Reference(Sref, cref, bref, rref, Vinf) - -# # freestream parameters -# alpha = 5.0*pi/180 -# beta = 0.0 -# Omega = [0.0; 0.0; 0.0] -# fs = Freestream(Vinf, alpha, beta, Omega) - -# # geometry -# xle = [0.0, 0.0] -# yle = [-bref/2, bref/2] -# zle = [0.0, 0.0] -# chord = [cref, cref] -# theta = [0.0, 0.0] -# phi = [0.0, 0.0] -# ns = 13 -# nc = 4 -# spacing_s = Uniform() -# spacing_c = Uniform() -# mirror = false -# symmetric = false - -# # non-dimensional time -# t = range(0.0, 10.0, step=0.2) -# dt = [t[i+1]-t[i] for i = 1:length(t)-1] - -# # create vortex rings -# grid, surface = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; -# mirror=mirror, spacing_s=spacing_s, spacing_c=spacing_c) - -# surfaces = [surface] - -# # run analysis -# system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; -# symmetric=symmetric) - -# # extract forces at each time step -# CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) -# end - -# @testset "Unsteady Vortex Lattice Method - Wing + Tail" begin +@testset "FMM acceleration" begin +function test_fmm(; fmm_toggle=true, fmm_direct=false, fmm_p=5, fmm_ncrit=20, fmm_theta=0.4) # Unsteady Wing and Tail ns_multiplier = 1 nc_multiplier = 1 - time_step = 6.4 # in [6.4, 3.2, 1.6, 0.8, 0.4, 0.2, 0.1, 0.05] + time_step = 0.05 # in [6.4, 3.2, 1.6, 0.8, 0.4, 0.2, 0.1, 0.05] # wing xle = [0.0, 0.2] @@ -1438,7 +1439,7 @@ end chord = [1.0, 0.6] theta = [2.0*pi/180, 2.0*pi/180] phi = [0.0, 0.0] - ns = 12 * ns_multiplier + ns = Int(12 * ns_multiplier) nc = 1 * nc_multiplier spacing_s = Uniform() spacing_c = Uniform() @@ -1451,11 +1452,11 @@ end chord_h = [0.7, 0.42] theta_h = [0.0, 0.0] phi_h = [0.0, 0.0] - ns_h = 6 * ns_multiplier + ns_h = Int(6 * ns_multiplier) nc_h = 1 * nc_multiplier spacing_s_h = Uniform() spacing_c_h = Uniform() - mirror_h = false + mirror_h = true # vertical stabilizer xle_v = [0.0, 0.14] @@ -1487,7 +1488,7 @@ end Omega = [0.0; 0.0; 0.0] fs = Freestream(Vinf, alpha, beta, Omega) - symmetric = [true, true, false] + symmetric = [false, false, false] # horseshoe vortices wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc; @@ -1505,57 +1506,36 @@ end surface_id = [1, 2, 3] # t - t = range(0.0, 10.0, step=time_step) + t = range(0.0, length=5, step=time_step) dt = t[2:end] - t[1:end-1] - println("Running VortexLattice (nsteps = $(length(t)))") - @time begin - system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; symmetric) - - # extract forces at each time step - CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) - end - # copy to fast-access vector - fast_access_vector = Vector{VortexLattice.FastMultipolePanel{Float64}}(undef,sum(length.(system.wakes))) - - function update_fmm_panels!(fmm_panels, wake_panels::Vector{Matrix{VortexLattice.WakePanel{TF}}}) where TF - i_start = 0 - for wake in wake_panels - for (i,panel) in enumerate(wake) - fmm_panels[i_start+i] = VortexLattice.FastMultipolePanel(panel) - end - i_start += length(wake) - end - end - - using BenchmarkTools - @show sum(length.(system.wakes)) - @btime update_fmm_panels!($fast_access_vector, $system.wakes) - -# end + system, surface_history, property_history, wake_history = unsteady_analysis(surfaces, ref, fs, dt; + symmetric, fmm_toggle, fmm_direct, fmm_p, fmm_ncrit, fmm_theta, fcore = (c, Δs) -> 1e-3, vertical_segments=false, + wake_finite_core=fill(false, length(surfaces))) -# @testset "velocity gradient" begin + # extract forces at each time step + CF, CM = body_forces_history(system, surface_history, property_history; frame=Wind()) -# # test velocity gradient of a vortex filament of unit strength -# function test_velocity_gradient(r1,r2) -# function velocity(x) -# return VortexLattice.bound_induced_velocity(r1+x, r2+x, false, 1e-2) -# end -# return ForwardDiff.jacobian(velocity, zeros(3)) -# end + return CF, CM, system, surface_history, property_history, wake_history +end -# r1 = SVector{3}(0.21490591034021655, -# 0.5811234218620989, -# 0.8079954052156225) +# test with FLOWFMM.direct! +CF_fmm, CM_fmm, system_fmm, surface_history_fmm, property_history_fmm, wake_history_fmm = test_fmm(; fmm_toggle=true, fmm_direct=true) +CF_direct, CM_direct, system_direct, surface_history_direct, property_history_direct, wake_history_direct = test_fmm(; fmm_toggle=false) -# r2 = SVector{3}(0.04054862116722213, -# 0.6134119927849612, -# 0.6270695997855336) +for i in eachindex(CF_fmm) + for j in 1:3 + @test isapprox(CF_fmm[i][j], CF_direct[i][j]; atol=1e-10) + end +end -# analytic_gradient = VortexLattice.bound_velocity_gradient(r1,r2) -# check = ForwardDiff.jacobian((x)->VortexLattice.bound_induced_velocity(r1+x,r2+x,false,1e-2),zero(SVector{3,Float64})) -# for i in eachindex(check) -# @show isapprox(analytic_gradient[i], check[i]; atol=1e-12) -# end +# test with FLOWFMM.fmm! +CF_fmm, CM_fmm, system_fmm, surface_history_fmm, property_history_fmm, wake_history_fmm = test_fmm(; fmm_toggle=true, fmm_direct=false) +CF_direct, CM_direct, system_direct, surface_history_direct, property_history_direct, wake_history_direct = test_fmm(; fmm_toggle=false) +for i in eachindex(CF_fmm) + for j in 1:3 + @test isapprox(CF_fmm[i][j], CF_direct[i][j]; atol=1e-3) + end +end -# end \ No newline at end of file +end