From a9997c1734e6f770f4fd00fad5e7bebbc14743ad Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 02:03:10 +0800 Subject: [PATCH 01/13] Merging vortex ring functionality --- .../3d_vortex_lattice_method/vlm_benchmark.jl | 110 +++++++++++ .../3d_vortex_lattice_method/vlm_implicit.jl | 120 ++++++++++++ .../3d_vortex_lattice_method/vlm_sand.jl | 123 ++++++++++++ .../vlm_vortex_rings.jl | 55 ++++++ src/AeroFuse.jl | 53 ++++-- .../VortexLattice/VortexLattice.jl | 14 +- src/Aerodynamics/VortexLattice/farfield.jl | 33 ++-- src/Aerodynamics/VortexLattice/horseshoes.jl | 151 --------------- .../VortexLattice/prandtl_glauert.jl | 9 + src/Aerodynamics/VortexLattice/printing.jl | 33 +++- .../VortexLattice/reference_frames.jl | 22 +-- src/Aerodynamics/VortexLattice/system.jl | 23 +-- .../VortexLattice/vortex_rings.jl | 60 ------ src/Aerodynamics/VortexLattice/vortices.jl | 179 ++++++++++++++++++ src/Aerodynamics/vlm_interface.jl | 58 ++++++ 15 files changed, 761 insertions(+), 282 deletions(-) create mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl create mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl create mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl create mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl delete mode 100644 src/Aerodynamics/VortexLattice/horseshoes.jl delete mode 100644 src/Aerodynamics/VortexLattice/vortex_rings.jl create mode 100644 src/Aerodynamics/VortexLattice/vortices.jl create mode 100644 src/Aerodynamics/vlm_interface.jl diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl new file mode 100644 index 00000000..891f2086 --- /dev/null +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl @@ -0,0 +1,110 @@ +## +using BenchmarkTools + +## BYU Flow +using VortexLattice + +function vlm_byu() + # 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 = 6 + spacing_s = VortexLattice.Uniform() + spacing_c = VortexLattice.Uniform() + + Sref = 30.0 + cref = 2.0 + bref = 15.0 + rref = [0.50, 0.0, 0.0] + Vinf = 1.0 + ref = VortexLattice.Reference(Sref, cref, bref, rref, Vinf) + + alpha = 1.0*pi/180 + beta = 0.0 + Omega = [0.0; 0.0; 0.0] + fs = VortexLattice.Freestream(Vinf, alpha, beta, Omega) + + # vortex rings with mirrored geometry + mirror = true + symmetric = false + + _, 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=VortexLattice.Stability()) + CDiff = far_field_drag(system) + + return CF, CM, CDiff, system +end + +## +@time CF, CM, CDiff, byu_sys = vlm_byu(); + +## +# @benchmark vlm_byu() + +## AeroFuse.jl +using AeroFuse + +## Wing section setup +function vlm_aerofuse() + # Wing with same settings + wing = Wing( + foils = fill(naca4((0,0,1,2)), 2), + chords = [2.2, 1.8], + twists = [2.0, 2.0], + spans = [7.5], + dihedrals = [0.], + sweeps = [3.0528], + w_sweep = 0., # Quarter-chord sweep + symmetry = true, + ) + + wing_mesh = WingMesh(wing, [24], 6, span_spacing = AeroFuse.Uniform()); + + # Freestream conditions + fs = Freestream( + alpha = 1.0, # deg + beta = 0.0, # deg + omega = [0.,0.,0.] + ) + + # Reference values + ref = References( + speed = 1., # m/s + density = 1.225, # kg/m³ + viscosity = 1.5e-5, # ??? + area = Sref, # m² + span = bref, # m + chord = cref, # m + location = rref # m + ) + + ## Horseshoes + ac_hs = ComponentVector(wing = make_horseshoes(wing_mesh)) + system = VortexLatticeSystem(ac_hs, fs, ref) + + ## Vortex rings + # ac_vs = ComponentVector(wing = make_vortex_rings(wing_mesh)) + # system = VortexLatticeSystem(ac_vs, fs, ref) + + return nearfield(system), farfield(system), system +end + +## +@time nfs, ffs, sys = vlm_aerofuse(); + +## +print_coefficients(nfs,ffs) + +## +@benchmark vlm_aerofuse() diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl new file mode 100644 index 00000000..63a779fa --- /dev/null +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl @@ -0,0 +1,120 @@ +## VLM implicit analysis and derivatives +using AeroFuse +using LinearAlgebra +using ImplicitAD +using NLsolve +using StructArrays + +## Surfaces +#================================================# + +function make_surfaces(x) + ## Wing + wing = Wing( + foils = fill(naca4(2,4,1,2), 2), + chords = [1.0, 0.6], + twists = [2.0, 0.0], + spans = [4.0], + dihedrals = [5.], + sweeps = [5.], + symmetry = true + ) + + ## Horizontal tail + htail = Wing( + foils = fill(naca4(0,0,1,2), 2), + chords = [0.7, 0.42], + twists = [0.0, 0.0], + spans = [1.25], + dihedrals = [0.], + sweeps = [6.39], + position = [4., 0, -0.1], + angle = -2., + axis = [0., 1., 0.], + symmetry = true + ) + + ## Vertical tail + vtail = Wing( + foils = fill(naca4(0,0,0,9), 2), + chords = [0.7, 0.42], + twists = [0.0, 0.0], + spans = [1.0], + dihedrals = [0.], + sweeps = [7.97], + position = [4., 0, 0], + angle = 90., + axis = [1., 0., 0.] + ) + + ## Meshing + #================================================# + + wing_mesh = WingMesh(wing, [24], 12, span_spacing = Cosine()) + htail_mesh = WingMesh(htail, [12], 6, span_spacing = Cosine()) + vtail_mesh = WingMesh(vtail, [12], 6, span_spacing = Cosine()) + + # Assemble + aircraft = ComponentArray( + wing = make_horseshoes(wing_mesh), + htail = make_horseshoes(htail_mesh), + vtail = make_horseshoes(vtail_mesh) + ) + + return aircraft +end + +ac = make_surfaces([1.0, 0.6]) + +## VLM analysis +#================================================# + +fs = Freestream( + alpha = 0.0, + beta = 0.0, + omega = [0., 0., 0.] +) + +ref = References( + speed = 150.0, + density = 1.225, + viscosity = 1.5e-5, +) + +function solve(system) + rwrap(R, p) = AeroFuse.VortexLattice.residual!(R, p, system) + res = nlsolve( + rwrap, + zero(system.circulations), + method = :newton, + autodiff = :forward, + show_trace = true + ) + return res.zero +end + +## +function program(x) + ac = make_surfaces(x) + fs = Freestream( + alpha = convert(eltype(x), 5.0), + beta = 0.0, + omega = [0., 0., 0.] + ) + + ref = References( + speed = 150.0, + density = 1.225, + viscosity = 1.5e-5, + ) + + sys = VortexLatticeSystem(ac, fs, ref) + + solve(sys) + # [ velocity(fs) sum(ac.rc) ] +end + +## +using ForwardDiff + +ForwardDiff.jacobian(program, [1.1, 0.6]) diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl new file mode 100644 index 00000000..1e2662cb --- /dev/null +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl @@ -0,0 +1,123 @@ +## +using AeroFuse +using JuMP +using Ipopt +using NLsolve + +## Wing section setup +wing = Wing( + foils = fill(naca4((2,4,1,2)), 3), + chords = [2.0, 1.6, 0.2], + twists = [0.0, 0.0, 0.0], + spans = [5., 0.6], + dihedrals = [5., 5.], + sweeps = [20.,20.], + w_sweep = 0.25, # Quarter-chord sweep + symmetry = true, + # flip = true +) + +## +x_w, y_w, z_w = wing_mac = mean_aerodynamic_center(wing) +print_info(wing, "Wing") + +## Meshing and assembly +wing_mesh = WingMesh(wing, [24,12], 6, + span_spacing = Cosine() + ); +aircraft = ComponentVector(wing = make_horseshoes(wing_mesh)) + +# Freestream conditions +fs = Freestream( + alpha = 2.0, # deg + beta = 2.0, # deg + omega = [0.,0.,0.] +) + +# Reference values +ref = References( + speed = 150, # m/s + density = 1.225, # kg/m³ + viscosity = 1.5e-5, # ??? + area = projected_area(wing), # m² + span = span(wing), # m + chord = mean_aerodynamic_chord(wing), # m + location = mean_aerodynamic_center(wing) # m +) + +## Solve system +system = VortexLatticeSystem(aircraft, fs, ref) + +## NLsolve +R = similar(system.circulations) + +## +AeroFuse.solve_nonlinear!(R, x) = solve_nonlinear!(R, system.vortices, x ./ ref.speed, velocity(fs) / ref.speed, fs.omega) + +## +xz = nlsolve( + solve_nonlinear!, + system.circulations, + autodiff = :forward, + show_trace = true, +) + +n = length(system.circulations) + +## THIS IS THE IMPLICIT FUNCTION THEOREM APPLIED, R(x) = R(x, u(x)) +obtain_root(x) = nlsolve(u -> solve_nonlinear!(R, x, u), system.circulations).zero # Rootfinding method (doesn't matter which) + + +## ChainRules adjoint +## Evaluate +x0 = [1.2, -2.0] +obtain_root(x0) +loss_functional(obtain_root(x0)) + +## Define adjoint rule +using ChainRulesCore + +function ChainRulesCore.rrule(::typeof(obtain_root), x) + # Solve residual equation + # R(x,u) =(by implicit function theorem)= u -> R(u(x)) = 0 + u_class = obtain_root(x) + + # Define pullback + function obtain_root_pullback(∂J_∂u) + # Evaluate gradients ∂R/∂x, ∂R/∂u + ∂R_∂x, ∂R_∂u = ForwardDiff.jacobian!(residual, x, u_class) + + # Adjoint variable: + # λᵀ = ∂J/∂R = (∂R/∂u)^(-1) * ∂J/∂u + λᵀ = ∂R_∂u' \ ∂J_∂u + + # Loss sensitivity to input: + # ∂J/∂x = -λ^T * ∂R/∂x + dJ_dx = -λᵀ' * ∂R_∂x + + # Loss sensitivity to residual equation: + # ∂J/∂R = 0 + ∂J_∂R = NoTangent() + + return ∂J_∂R, dJ_dx' + end + + return u_class, obtain_root_pullback +end + + + + +## Create JuMP model, using Ipopt as the solver +model = Model(optimizer_with_attributes(Ipopt.Optimizer)) + +@variables(model, begin + Γs[j = 1:n] == system.circulations[j] +end) + +register(model, :R_A, n, solve_nonlinear; autodiff = true) + +# @NLconstraint(model, [j = 1:n], solve_nonlinear(Γs)[i] == 0) + +## +optimize!(model) \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl new file mode 100644 index 00000000..0fe17ce4 --- /dev/null +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl @@ -0,0 +1,55 @@ +## Wing analysis case +using AeroFuse +using Plots +using StructArrays + +## Wing section setup +wing = Wing( + foils = fill(naca4((0,0,1,2)), 2), + chords = [2.2, 1.8], + twists = [0.0, 0.0], + spans = [7.5], + dihedrals = [0.], + sweeps = [3.0528], + w_sweep = 0., # Quarter-chord sweep + symmetry = true, +) + +wing_mesh = WingMesh(wing, [24], 6, span_spacing = Uniform()); + +# Freestream conditions +fs = Freestream( + alpha = 1.0, # deg + beta = 0.0, # deg + omega = [0.,0.,0.] +) + +# Reference values +ref = References( + speed = 15., # m/s + density = 1.225, # kg/m³ + viscosity = 1.5e-5, # ??? + area = projected_area(wing), # m² + span = span(wing), # m + chord = mean_aerodynamic_chord(wing), # m + location = [0.5,0.,0.] # m +) + +## Horseshoes +ac_hs = ComponentVector(wing = make_horseshoes(wing_mesh)) +system = VortexLatticeSystem(ac_hs, fs, ref) +print_coefficients(nearfield(system), farfield(system)) + +## Vortex rings +ac_vs = ComponentVector(wing = make_vortex_rings(wing_mesh)) +sys = VortexLatticeSystem(ac_vs, fs, ref) +print_coefficients(nearfield(sys), farfield(sys)) + +## Manual +A = influence_matrix(ac_vs) +b = boundary_condition(ac_vs, -velocity(fs), fs.omega) +gam = A \ b + +## +sys = VortexLatticeSystem(ac_vs, gam * ref.speed, A, b, fs, ref, Wind()) +print_coefficients(nearfield(sys), farfield(sys)) diff --git a/src/AeroFuse.jl b/src/AeroFuse.jl index e4c601d8..5e0ed7bb 100644 --- a/src/AeroFuse.jl +++ b/src/AeroFuse.jl @@ -20,11 +20,10 @@ export combinedimsview, combinedims, splitdimsview, splitdims using ComponentArrays export ComponentVector, ComponentArray -using Accessors -export @set +using LabelledArrays using Accessors -using LabelledArrays +export @set ## Methods to be extended in submodules #==========================================================================================# @@ -100,9 +99,9 @@ import .AircraftGeometry: AbstractAircraft, AbstractWing, AbstractFoil export AbstractAircraft, AbstractWing, AbstractFoil # Foil -import .AircraftGeometry: Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_interpolation, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, interpolate, rotate, affine, scale, reflect, camber, camber_line, control_surface +import .AircraftGeometry: Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_interpolation, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, interpolate, rotate, affine, scale, reflect, camber, camber_line, thickness_line, control_surface -export Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_interpolation, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, interpolate, rotate, affine, scale, reflect, camber, camber_line, control_surface +export Foil, arc_length, kulfan_CST, naca4, camber_CST, make_panels, read_foil, leading_edge_index, upper_surface, lower_surface, split_surface, coordinates_to_camber_thickness, camber_thickness_to_coordinates, camber_thickness, camber_thickness_to_coordinates, cosine_interpolation, camber_thickness_to_CST, coordinates_to_CST, maximum_thickness_to_chord, translate, interpolate, rotate, affine, scale, reflect, camber, camber_line, thickness_line, control_surface # Fuselage import .AircraftGeometry: Fuselage, projected_area, length, cosine_interpolation, volume @@ -129,13 +128,42 @@ export Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, a make_horseshoes(wing :: WingMesh) = map((cho,cam) -> Horseshoe(cho, normal_vector(cam)), chord_panels(wing), camber_panels(wing)) function make_vortex_rings(wing :: WingMesh) - cams = camber_panels(wing) + cam_coo = camber_coordinates(wing) + cams = combinedimsview(cam_coo, (1,2)) + + # Generate vortex ring mesh + vor_cams = similar(cams) + vor_cams = 0.75 * cams[1:end-1,:,:] + 0.25 * cams[2:end,:,:] + vor_cams[end,:,:] = cams[end,:,:] - # Vortex rings on surface - rings = VortexRing.(cams) + @views rings = VortexRing.(make_panels(splitdimsview(vor_cams, (1,2)))) + + # NOTE: JUST ADD A TRAILING TYPE TO THE VORTEX RING AND ADD THE HORSESHOE VELOCITY METHOD # Trailing edge semi-infinite vortex lines - horsies = Horseshoe.(cams[end,:], normal_vector.(cams[end,:])) + cam_pan = camber_panels(wing) + @views horsies = Horseshoe.(cam_pan[end,:], normal_vector.(cam_pan[end,:])) + + # Generate vortex ring mesh + # vor_cams = similar(cams) + # vor_cams[1:end-1,:,:] = 0.75 * cams[1:end-1,:,:] + 0.25 * cams[2:end,:,:] + # vor_cams[end,:,:] = cams[end,:,:] + + # # display(vor_cams) + + # vor_coo = splitdimsview(vor_cams, (1,2)) + # vor_pan = make_panels(vor_coo) + # @views rings = @. VortexRing(vor_pan) + + # # Trailing edge semi-infinite vortex lines + # hor_cams = similar(cams[end-1:end,:,:]) + # hor_cams[end-1,:,:] = vor_cams[end,:,:] + # hor_cams[end,:,:] = (cams[end,:,:] - cams[end-1,:,:]) * 0.25 + cams[end,:,:] + + # hor_coo = splitdimsview(hor_cams, (1,2)) + # hor_pan = make_panels(hor_coo) + # @views horsies = @. Horseshoe(hor_pan[end,:], normal_vector(hor_pan[end,:])) + # Assemble [ rings; permutedims(horsies) ] @@ -163,9 +191,12 @@ export total_velocity, source_velocity, vortex_velocity, vortex_influence_matrix ## Vortex lattice include("Aerodynamics/VortexLattice/VortexLattice.jl") -import .VortexLattice: Horseshoe, AbstractVortexLatticeSystem, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, transform, bound_leg, bound_leg_center, bound_leg_vector, r1, r2, control_point, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes, rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, velocity, kinematic_viscosity, mach_number, stream_velocity, center_of_pressure, freestream_derivatives!, freestream_derivatives, print_coefficients, print_derivatives +import .VortexLattice: Horseshoe, AbstractVortexLatticeSystem, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, bound_leg_center, bound_leg_vector, control_point, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes, rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, velocity, kinematic_viscosity, mach_number, stream_velocity, center_of_pressure, freestream_derivatives!, freestream_derivatives, print_coefficients, print_derivatives + +export Horseshoe, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, bound_leg_center, bound_leg_vector, control_point, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes, rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, velocity, kinematic_viscosity, mach_number, stream_velocity, center_of_pressure, freestream_derivatives!, freestream_derivatives, print_coefficients, print_derivatives -export Horseshoe, VortexLatticeSystem, References, AbstractAxisSystem, Stability, Wind, Body, Geometry, streamlines, influence_coefficient, influence_matrix, boundary_condition, solve_system, transform, bound_leg, bound_leg_center, bound_leg_vector, r1, r2, control_point, Horseshoe, surface_velocity, surface_forces, surface_moments, nearfield_drag, geometry_to_wind_axes, geometry_to_stability_axes, stability_to_geometry_axes, wind_to_geometry_axes, rate_coefficient, nearfield, farfield, farfield_forces, surface_velocities, surface_forces, surface_dynamics, surface_coefficients, nearfield_coefficients, farfield_coefficients, VortexRing, velocity, kinematic_viscosity, mach_number, stream_velocity, center_of_pressure, freestream_derivatives!, freestream_derivatives, print_coefficients, print_derivatives +## Panel-VLM interface +include("Aerodynamics/vlm_interface.jl") ## Profile drag estimation include("Aerodynamics/profile_drag.jl") diff --git a/src/Aerodynamics/VortexLattice/VortexLattice.jl b/src/Aerodynamics/VortexLattice/VortexLattice.jl index 92607e43..2c18ea5f 100644 --- a/src/Aerodynamics/VortexLattice/VortexLattice.jl +++ b/src/Aerodynamics/VortexLattice/VortexLattice.jl @@ -14,9 +14,6 @@ using DiffResults: JacobianResult, jacobian, value ## Package imports #==========================================================================================# -# Panel geometry -import ..PanelGeometry: Panel3D, panel_area, panel_coordinates, midpoint, normal_vector, transform, p1, p2, p3, p4, average_chord, average_width - # Non-dimensionalization import ..NonDimensional: dynamic_pressure, aerodynamic_coefficients, force_coefficient, moment_coefficient, rate_coefficient @@ -28,9 +25,7 @@ import ..AeroFuse: velocity, solve_system, solve_linear, solve_nonlinear, solve_ ## Vortex types and methods #==========================================================================================# -include("horseshoes.jl") - -include("vortex_rings.jl") +include("vortices.jl") ## Reference frames #==========================================================================================# @@ -43,11 +38,6 @@ struct Body <: AbstractAxisSystem end struct Stability <: AbstractAxisSystem end struct Wind <: AbstractAxisSystem end -Base.show(io :: IO, :: Geometry) = print(io, "Geometry") -Base.show(io :: IO, :: Body) = print(io, "Body") -Base.show(io :: IO, :: Stability) = print(io, "Stability") -Base.show(io :: IO, :: Wind) = print(io, "Wind") - """ velocity(freestream :: Freestream, ::Body) @@ -87,7 +77,6 @@ include("nearfield.jl") # Farfield forces include("farfield.jl") - # System setups #==========================================================================================# @@ -106,5 +95,4 @@ include("printing.jl") # Streamlines include("streamlines.jl") - end \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/farfield.jl b/src/Aerodynamics/VortexLattice/farfield.jl index 24af9961..7a40a8b7 100644 --- a/src/Aerodynamics/VortexLattice/farfield.jl +++ b/src/Aerodynamics/VortexLattice/farfield.jl @@ -6,45 +6,52 @@ Compute the induced velocity at point ``r_i`` of the wake in the Trefftz plane due to the strength ``Γ_j`` at ``r_j``. """ -farfield_velocity(r_i, r_j, Γ_j) = let r = r_i - r_j; Γ_j / 2π * SVector(1., 0., 0.) × r / norm(r)^2 end +function farfield_velocity(r_i, r_j, Γ_j) + T = promote_type(eltype(r_i), eltype(r_j), typeof(Γ_j)) + r = r_i - r_j + + Γ_j / 2π * SVector{3,T}(1, 0, 0) × r / norm(r)^2 +end """ farfield_influence_matrix(centers, normals, points) Compute the aerodynamic influence coefficient matrix of the wake in the Trefftz plane given the center points, normal vectors, and the points of the wake panels. """ -farfield_influence_matrix(centers, normals, points) = [ dot(farfield_velocity(r_i, r_j2, 1.) - farfield_velocity(r_i, r_j1, 1.), n_i) for (r_i, n_i) in zip(centers, normals), (r_j2, r_j1) in zip(points[2:end], points) ] +@views farfield_influence_matrix(centers, normals, points) = [ dot(farfield_velocity(r_i, r_j2, 1.) - farfield_velocity(r_i, r_j1, 1.), n_i) for (r_i, n_i) in zip(centers, normals), (r_j2, r_j1) in zip(points[2:end], points) ] """ doublet_normal_derivatives(wake_points, Δφs, normals) Compute the normal derivative strengths of the doublets given the wake points, net doublet strengths ``Δφ``s, and associated normal vectors. """ -function doublet_normal_derivatives(wake_points, Δφs, normals) +@views function doublet_normal_derivatives(wake_points, Δφs, normals) centers = @. (wake_points[1:end-1] + wake_points[2:end]) / 2 AIC = farfield_influence_matrix(centers, normals, wake_points) ∂φ_∂n = AIC * Δφs + + return ∂φ_∂n end project_vector(vector, U) = vector - dot(U, vector) * U normal(wake_proj_vector, U) = normalize(U × wake_proj_vector) -dihedral(wake_proj_vec) = atan(wake_proj_vec[3], wake_proj_vec[2]) +@views dihedral(wake_proj_vec) = atan(wake_proj_vec[3], wake_proj_vec[2]) """ - trefftz_plane_quantities(horseshoes :: AbstractArray{<: Horseshoe}, α, β) + trefftz_plane_quantities(vortices, α, β) Project `Horseshoe`s into the Trefftz plane aligned with the wind axes angles ``α,~β``, and compute normal vectors, projection angles and lengths. """ -function trefftz_plane_quantities(horseshoes :: AbstractArray{<: Horseshoe}, α, β) +function trefftz_plane_quantities(vortices, α, β) # Reference velocity for broadcasting U_ref = (Ref ∘ SVector)(1, 0, 0) # Transform wake to wind axes - wake_horsies = @views vec(horseshoes[end,:]) # Trailing edge + wake_horsies = @views vec(vortices[end,:]) # Trailing edge wake_points = [ r1.(wake_horsies); [r2(wake_horsies[end])] ] # Trailing edge points wake_wind = geometry_to_wind_axes.(wake_points, α, β) # Wind axis transformation - # Project trailing edge horseshoes' bound legs into Trefftz plane along wind axes + # Project trailing edge vortices' bound legs into Trefftz plane along wind axes wake_vectors = @. wake_wind[2:end] - wake_wind[1:end-1] # Vectors wake_proj_vecs = @. project_vector(wake_vectors, U_ref) # Project to wind axes @@ -53,7 +60,7 @@ function trefftz_plane_quantities(horseshoes :: AbstractArray{<: Horseshoe}, α, wake_dihedrals = @. dihedral(wake_proj_vecs) wake_lengths = @. norm(wake_proj_vecs) - wake_points, wake_normals, wake_dihedrals, wake_lengths + return wake_points, wake_normals, wake_dihedrals, wake_lengths end """ @@ -63,8 +70,8 @@ Compute the aerodynamic forces in the Trefftz plane given cumulative doublet str """ function compute_farfield_forces(Δφs, Δs, ∂φ_∂n, θs, V, ρ) D_i = - 1/2 * ρ * sum(@. Δφs * Δs * ∂φ_∂n) - Y = - ρ * V * sum(@. Δφs * Δs * sin(θs)) - L = ρ * V * sum(@. Δφs * Δs * cos(θs)) + Y = - ρ * V * sum(@. Δφs * Δs * sin(θs)) + L = ρ * V * sum(@. Δφs * Δs * cos(θs)) SVector(D_i, Y, L) end @@ -79,9 +86,9 @@ function farfield_forces(Γs, horseshoes, speed, α, β, ρ) wake_points, normals, dihedrals, Δs = trefftz_plane_quantities(horseshoes, α, β) # Compute directional derivatives of doublets in the normal direction - Δφs = vec(sum(Γs, dims = 1)) + Δφs = vec(sum(Γs, dims = 1)) ∂φ_∂n = doublet_normal_derivatives(wake_points, Δφs, normals) # Compute forces - compute_farfield_forces(Δφs, Δs, ∂φ_∂n, dihedrals, speed, ρ) + return compute_farfield_forces(Δφs, Δs, ∂φ_∂n, dihedrals, speed, ρ) end \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/horseshoes.jl b/src/Aerodynamics/VortexLattice/horseshoes.jl deleted file mode 100644 index 85e25806..00000000 --- a/src/Aerodynamics/VortexLattice/horseshoes.jl +++ /dev/null @@ -1,151 +0,0 @@ -## Horseshoe methods -#==========================================================================================# - -function quarter_point(p1, p2) - μ = SVector(1/4, 0, 1/4) - return @. (1 - μ) * p1 + μ * p2 -end - -function three_quarter_point(p1, p2) - μ = SVector(3/4, 0, 3/4) - return @. (1 - μ) * p1 + μ * p2 -end - -control_point(p1, p2, p3, p4) = ( three_quarter_point(p1, p2) + three_quarter_point(p4, p3) ) / 2 -bound_leg(p1, p2, p3, p4) = SVector(quarter_point(p1, p2), quarter_point(p4, p3)) - -""" - bound_leg(panel :: Panel3D) - -Compute the bound leg for a `Panel3D`, for horseshoes/vortex rings. -""" -bound_leg(panel :: Panel3D) = bound_leg(panel.p1, panel.p2, panel.p3, panel.p4) - -""" - control_point(panel :: Panel3D) - -Compute the control point of a `Panel3D` for horseshoes/vortex rings, which is the 3-quarter point on each side in the ``x``-``z`` plane. -""" -control_point(panel :: Panel3D) = control_point(panel.p1, panel.p2, panel.p3, panel.p4) - -# Velocity kernels -bound_leg_velocity(a, b, Γ) = Γ/4π * (1/norm(a) + 1/norm(b)) * a × b / (norm(a) * norm(b) + dot(a, b)) -trailing_leg_velocity(r, Γ, u) = Γ/4π * normalize(r) × normalize(u) / (norm(r) - dot(r, u)) - -trailing_legs_velocities(a, b, Γ, u) = trailing_leg_velocity(a, Γ, u) - trailing_leg_velocity(b, Γ, u) -total_horseshoe_velocity(a, b, Γ, u) = bound_leg_velocity(a, b, Γ) + trailing_legs_velocities(a, b, Γ, u) - -# Finite-core velocity kernels -function bound_leg_velocity(a, b, Γ, ε) - na, nb, σ = norm(a), norm(b), dot(a, b) - term_1 = (na^2 - σ) / √(na^2 + ε^2) + (nb^2 - σ) / √(nb^2 + ε^2) - term_2 = a × b / (na^2 * nb^2 - σ^2 + ε^2 * (na^2 + nb^2 - 2 * na * nb)) - - Γ/4π * term_1 * term_2 -end - -trailing_leg_velocity(r, Γ, u, ε) = Γ/4π * normalize(r) × u / (norm(r) - dot(r, u) + ε^2 / (norm(r) + dot(r, u))) -trailing_legs_velocities(a, b, Γ, u, ε) = trailing_leg_velocity(a, Γ, u, ε) - trailing_leg_velocity(b, Γ, u, ε) -total_horseshoe_velocity(a, b, Γ, u, ε) = bound_leg_velocity(a, b, Γ, ε) + trailing_legs_velocities(a, b, Γ, u, ε) - -## Arrays of vortex lines -#==========================================================================================# - -abstract type AbstractVortex end - -""" - Horseshoe(r1, r2, rc, normal, chord) - -Define a horseshoe vortex with a start and endpoints ``r₁, r₂`` for the bound leg, a collocation point ``r``, a normal vector ``n̂``, and a finite core size. - -The finite core setup is not implemented for now. -""" -struct Horseshoe{T <: Real} <: AbstractVortex - r1 :: SVector{3,T} - r2 :: SVector{3,T} - rc :: SVector{3,T} - normal :: SVector{3,T} - core :: T -end - -Base.length(::Horseshoe) = 1 - -r1(horseshoe :: Horseshoe) = horseshoe.r1 -r2(horseshoe :: Horseshoe) = horseshoe.r2 - -function Horseshoe(r1, r2, rc, n, c) - T = promote_type(eltype(r1), eltype(r2), eltype(rc), eltype(n), eltype(c)) - Horseshoe{T}(r1, r2, rc, n, c) -end - -""" - bound_leg(horseshoe :: Horseshoe) - -Getter for bound leg field of a `Horseshoe`. -""" -bound_leg(horseshoe :: Horseshoe) = horseshoe.bound_leg -control_point(horseshoe :: Horseshoe) = horseshoe.rc -normal_vector(horseshoe :: Horseshoe) = horseshoe.normal - -r1(r, horseshoe :: Horseshoe) = r - horseshoe.r1 -r2(r, horseshoe :: Horseshoe) = r - horseshoe.r2 - -""" - Horseshoe(panel :: Panel3D, normal, drift = zeros(3)) - -Generate a `Horseshoe` corresponding to a `Panel3D`, an associated normal vector, and a "drift velocity". -""" -function Horseshoe(panel :: Panel3D, normal, drift = SVector(0., 0., 0.); core_size = 0.) - r1, r2 = bound_leg(panel) - rc = control_point(panel) + drift - Horseshoe(r1, r2, rc, normal, core_size) -end - -""" - transform(horseshoe :: Horseshoe, T :: LinearMap) - -Generate a new `Horseshoe` with the points and normal vectors transformed by the `LinearMap` ``T``. -""" -transform(horseshoe :: Horseshoe, T :: LinearMap) = setproperties(horseshoe, - r1 = T(horseshoe.r1), - r2 = T(horseshoe.r2), - rc = T(horseshoe.rc), - normal = T(horseshoe.normal), -) - -transform(horseshoe :: Horseshoe; rotation = I(3), translation = zeros(3)) = transform(horseshoe, Translation(translation) ∘ LinearMap(rotation)) - -""" - bound_leg_center(horseshoe :: Horseshoe) - -Compute the midpoint of the bound leg of a `Horseshoe`. -""" -bound_leg_center(horseshoe :: Horseshoe) = (horseshoe.r1 + horseshoe.r2) / 2 - -""" - bound_leg_vector(horseshoe :: Horseshoe) - -Compute the direction vector of the bound leg of a `Horseshoe`. -""" -bound_leg_vector(horseshoe :: Horseshoe) = horseshoe.r2 - horseshoe.r1 - -""" - velocity(r, horseshoe, Γ, u_hat = [1.,0.,0.]) - -Compute the induced velocity at a point ``r`` of a given `Horseshoe` with a bound leg of constant strength ``Γ`` and semi-infinite trailing legs pointing in a given direction ``û``, by default `û = x̂`. -""" -velocity(r, horseshoe :: Horseshoe, Γ :: Real, V_hat = SVector(1.,0.,0.)) = total_horseshoe_velocity(r1(r, horseshoe), r2(r, horseshoe), Γ, V_hat, horseshoe.core) - -""" - bound_velocity(r, horseshoe, Γ, u_hat) - -Compute the induced velocity at a point ``r`` from the bound leg with constant strength ``Γ`` of a given `Horseshoe`. -""" -bound_velocity(r, horseshoe :: Horseshoe, Γ) = bound_leg_velocity(r1(r, horseshoe), r2(r, horseshoe), Γ, horseshoe.core) - -""" - trailing_velocity(r, horseshoe, Γ, u_hat) - -Compute the induced velocity at a point ``r`` from the semi-infinite trailing legs with constant strength ``Γ`` of a given `Horseshoe`. -""" -trailing_velocity(r, horseshoe :: Horseshoe, Γ, V) = trailing_legs_velocities(r1(r, horseshoe), r2(r, horseshoe), Γ, V, horseshoe.core) \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/prandtl_glauert.jl b/src/Aerodynamics/VortexLattice/prandtl_glauert.jl index 5bd4a3c3..b2919e65 100644 --- a/src/Aerodynamics/VortexLattice/prandtl_glauert.jl +++ b/src/Aerodynamics/VortexLattice/prandtl_glauert.jl @@ -14,5 +14,14 @@ prandtl_glauert_scale_coordinates(horseshoe :: Horseshoe, β) = setproperties(ho normal = prandtl_glauert_scale_normal(horseshoe.normal, β), ) +prandtl_glauert_scale_coordinates(ring :: VortexRing, β) = setproperties(ring, + r1 = prandtl_glauert_scale_coordinates(ring.r1, β), + r2 = prandtl_glauert_scale_coordinates(ring.r2, β), + r3 = prandtl_glauert_scale_coordinates(ring.r3, β), + r4 = prandtl_glauert_scale_coordinates(ring.r4, β), + rc = prandtl_glauert_scale_coordinates(ring.rc, β), + normal = prandtl_glauert_scale_normal(ring.normal, β), +) + prandtl_glauert_inverse_scale_velocity(vx, vy, vz, β) = SVector(vx / β^2, vy / β, vz / β) prandtl_glauert_inverse_scale_velocity(v, β) = @views prandtl_glauert_inverse_scale_velocity(v[1], v[2], v[3], β) \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/printing.jl b/src/Aerodynamics/VortexLattice/printing.jl index 4acd5bd6..b072976f 100644 --- a/src/Aerodynamics/VortexLattice/printing.jl +++ b/src/Aerodynamics/VortexLattice/printing.jl @@ -1,6 +1,37 @@ -## Pretty-printing +## Pretty-printing and show methods #==========================================================================================# +# Axis systems +Base.show(io :: IO, :: Geometry) = print(io, "Geometry") +Base.show(io :: IO, :: Body) = print(io, "Body") +Base.show(io :: IO, :: Stability) = print(io, "Stability") +Base.show(io :: IO, :: Wind) = print(io, "Wind") + +# Reference values +function Base.show(io :: IO, refs :: References) + println(io, "References: ") + for fname in fieldnames(typeof(refs)) + println(io, " ", fname, " = ", getfield(refs, fname)) + end +end + +function Base.show(io :: IO, ring :: AbstractVortex) + println(io, "Vortex: ") + for fname in fieldnames(typeof(ring)) + println(io, " ", fname, " = ", getfield(ring, fname)) + end +end + +# Vortex lattice system +function Base.show(io :: IO, sys :: VortexLatticeSystem) + println(io, "VortexLatticeSystem -") + println(io, length(sys.vortices), " ", eltype(sys.vortices), " Elements") + println(io, "Axes: ", sys.axes) + show(io, sys.freestream) + show(io, sys.reference) +end + + """ print_coefficients(nf_coeffs, ff_coeffs, name = "") diff --git a/src/Aerodynamics/VortexLattice/reference_frames.jl b/src/Aerodynamics/VortexLattice/reference_frames.jl index 6594b059..70326ad8 100644 --- a/src/Aerodynamics/VortexLattice/reference_frames.jl +++ b/src/Aerodynamics/VortexLattice/reference_frames.jl @@ -3,19 +3,15 @@ """ reflect_xz(vector) - reflect_xz(line :: Line) - reflect_xz(horseshoe :: Horseshoe) -Reflect the ``y``-coordinate of a given 3-dimensional vector (or `Line` or `Horseshoe`) about the ``x``-``z`` plane. +Reflect the ``y``-coordinate of a given 3-dimensional vector about the ``x``-``z`` plane. """ reflect_xz(vector) = SVector(vector[1], -vector[2], vector[3]) -reflect_xz(horseshoe :: Horseshoe) = (Horseshoe ∘ reflect_xz ∘ bound_leg)(horseshoe) """ project_yz(vector) - project_yz(line :: Line) -Project a given 3-dimensional vector or `Line` into the ``y``-``z`` plane. +Project a given 3-dimensional vector or into the ``y``-``z`` plane. """ project_yz(vector) = SVector(0, vector[2], vector[3]) @@ -47,23 +43,23 @@ stability_to_geometry_axes(coords, α :: T) where T <: Real = geometry_to_stabil """ geometry_to_wind_axes(coords, α, β) - geometry_to_wind_axes(horseshoe :: Horseshoe, α, β) + geometry_to_wind_axes(vor :: AbstractVortex, α, β) Convert coordinates from geometry axes to wind axes for given angles of attack ``α`` and sideslip ``\\beta.`` """ geometry_to_wind_axes(coords, α, β) = let T = promote_type(eltype(α), eltype(β)); RotZY{T}(β, α) * coords end -function geometry_to_wind_axes(horseshoe :: Horseshoe, α, β) +function geometry_to_wind_axes(vortex :: AbstractVortex, α, β) T = promote_type(eltype(α), eltype(β)) - return transform(horseshoe, LinearMap(RotZY{T}(β, α))) + return transform(vortex, LinearMap(RotZY{T}(β, α))) end geometry_to_wind_axes(coords, fs :: Freestream) = geometry_to_wind_axes(coords, fs.alpha, fs.beta) -geometry_to_wind_axes(horseshoe :: Horseshoe, fs :: Freestream) = geometry_to_wind_axes(horseshoe, fs.alpha, fs.beta) +geometry_to_wind_axes(vor :: AbstractVortex, fs :: Freestream) = geometry_to_wind_axes(vor, fs.alpha, fs.beta) """ wind_to_geometry_axes(coords, α, β) - wind_to_geometry_axes(horseshoe :: Horseshoe, α, β) + wind_to_geometry_axes(vor :: AbstractVortex, α, β) Convert coordinates from wind axes to geometry axes for given angles of attack ``α`` and sideslip \\beta.`` """ @@ -73,9 +69,9 @@ function wind_to_geometry_axes(coords, α, β) return RotYZ{T}(-α, -β) * coords end -function wind_to_geometry_axes(horseshoe :: Horseshoe, α, β) +function wind_to_geometry_axes(vor :: AbstractVortex, α, β) T = promote_type(eltype(α), eltype(β)) - return transform(horseshoe, LinearMap(RotYZ{T}(-α, -β))) + return transform(vor, LinearMap(RotYZ{T}(-α, -β))) end diff --git a/src/Aerodynamics/VortexLattice/system.jl b/src/Aerodynamics/VortexLattice/system.jl index b9681a28..2acf3fd5 100644 --- a/src/Aerodynamics/VortexLattice/system.jl +++ b/src/Aerodynamics/VortexLattice/system.jl @@ -15,8 +15,8 @@ Define reference values with speed ``V``, density ``ρ``, _dynamic_ viscosity `` - `density :: Real = 1.225`: Density (m) - `viscosity :: Real = 1.5e-5`: Dynamic viscosity (m) - `sound_speed :: Real = 330.`: Speed of sound (m/s) -- `span :: Real = 1.`: Span length (m) - `area :: Real = 1.`: Area (m²) +- `span :: Real = 1.`: Span length (m) - `chord :: Real = 1.`: Chord length (m) - `location :: Vector{Real} = [0,0,0]`: Position """ @@ -42,13 +42,6 @@ References(; speed = 1., density = 1.225, viscosity = 1.5e-5, sound_speed = 330. Base.broadcastable(refs :: References) = Ref(refs) -function Base.show(io :: IO, refs :: References) - println(io, "References: ") - for fname in fieldnames(typeof(refs)) - println(io, " ", fname, " = ", getfield(refs, fname)) - end -end - dynamic_pressure(refs :: References) = 1/2 * refs.density * refs.speed^2 kinematic_viscosity(refs :: References) = refs.viscosity / refs.density @@ -100,7 +93,7 @@ function VortexLatticeSystem(components, fs :: Freestream, refs :: References, a # Mach number bound checks M = mach_number(refs) @assert M < 1. "Only compressible subsonic flow conditions (M < 1) are valid!" - if M > 0.7 @warn "Results in transonic flow conditions (0.7 < M < 1) are most likely incorrect!" end + if M > 0.7 @warn "Results in transonic to sonic flow conditions (0.7 < M < 1) are most likely incorrect!" end # (Prandtl-Glauert ∘ Wind axis) transformation β_pg = √(1 - M^2) @@ -116,16 +109,6 @@ function VortexLatticeSystem(components, fs :: Freestream, refs :: References, a return VortexLatticeSystem(components, refs.speed * Γs / β_pg^2, AIC, boco, fs, refs, axes) end -# Made out of annoyance and boredom -function Base.show(io :: IO, sys :: VortexLatticeSystem) - println(io, "VortexLatticeSystem -") - println(io, length(sys.vortices), " ", eltype(sys.vortices), " Elements\n") - println(io, " Axes: ", sys.axes) - show(io, sys.freestream) - println(io, "") - show(io, sys.reference) -end - # Miscellaneous rate_coefficient(system :: VortexLatticeSystem) = rate_coefficient(system.freestream, system.reference) @@ -334,4 +317,4 @@ function center_of_pressure(system :: VortexLatticeSystem) return x_CP end -residual!(R, system :: VortexLatticeSystem) = solve_nonlinear!(R, system.vortices, system.circulations, -velocity(system.freestream), system.freestream.omega) \ No newline at end of file +residual!(R, Γ, system :: VortexLatticeSystem) = solve_nonlinear!(R, system.vortices, Γ, -velocity(system.freestream), system.freestream.omega) \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/vortex_rings.jl b/src/Aerodynamics/VortexLattice/vortex_rings.jl deleted file mode 100644 index 2b7ec395..00000000 --- a/src/Aerodynamics/VortexLattice/vortex_rings.jl +++ /dev/null @@ -1,60 +0,0 @@ -""" - VortexRing(r1, r2, r3, r4, r_c, n̂, ε) - -A vortex ring consisting of four points ``r_i, i = 1,…,4``, a collocation point ``r_c``, a normal vector ``n̂``, and a core size ``ε``. -""" -struct VortexRing{T <: Real} <: AbstractVortex - r1 :: SVector{3,T} - r2 :: SVector{3,T} - r3 :: SVector{3,T} - r4 :: SVector{3,T} - rc :: SVector{3,T} - normal :: SVector{3,T} - core :: T -end - -control_point(ring :: VortexRing) = ring.rc -normal_vector(ring :: VortexRing) = ring.normal - -function vortex_lines(p1, p2, p3, p4) - v1 = quarter_point(p1, p2) - v4 = quarter_point(p4, p3) - v2 = v1 + p2 - p1 - v3 = v4 + p3 - p4 - - v1, v2, v3, v4 -end - -""" -Constructor for vortex rings on a Panel3D using Lines. The following convention is adopted: - -``` - p1 —front leg→ p4 - | | -left leg right leg - ↓ ↓ - p2 —back leg-→ p3 -``` -""" -function VortexRing(panel :: Panel3D{T}, normal = normal_vector(panel), ε = 0.) where T <: Real - r_c = control_point(panel) - VortexRing{T}(panel.p1, panel.p2, panel.p3, panel.p4, r_c, normal, ε) -end - -Base.length(::VortexRing) = 1 - -""" -Computes the induced velocities at a point ``r`` of a VortexRing with constant strength ``Γ``. -""" -function velocity(r, ring :: VortexRing, Γ) - v1 = bound_leg_velocity(r - ring.r1, r - ring.r2, Γ, ring.core) - v2 = bound_leg_velocity(r - ring.r2, r - ring.r3, Γ, ring.core) - v3 = bound_leg_velocity(r - ring.r3, r - ring.r4, Γ, ring.core) - v4 = bound_leg_velocity(r - ring.r4, r - ring.r1, Γ, ring.core) - - return v1 + v2 + v3 + v4 -end - -# influence_matrix(rings) = [ influence_coefficient(ring_j, ring_i) for ring_j in rings, ring_i in rings ] - -# boundary_condition(rings, U, Ω) = map(hs -> dot(U + Ω × control_point(hs), ring_normal(hs)), rings) \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/vortices.jl b/src/Aerodynamics/VortexLattice/vortices.jl new file mode 100644 index 00000000..c7424052 --- /dev/null +++ b/src/Aerodynamics/VortexLattice/vortices.jl @@ -0,0 +1,179 @@ +# Velocity kernels +#==========================================================================================# + +bound_leg_velocity(a, b, Γ) = Γ/4π * (1/norm(a) + 1/norm(b)) * a × b / (norm(a) * norm(b) + dot(a, b)) +trailing_leg_velocity(r, Γ, u) = Γ/4π * normalize(r) × normalize(u) / (norm(r) - dot(r, u)) + +trailing_legs_velocities(a, b, Γ, u) = trailing_leg_velocity(a, Γ, u) - trailing_leg_velocity(b, Γ, u) +total_horseshoe_velocity(a, b, Γ, u) = bound_leg_velocity(a, b, Γ) + trailing_legs_velocities(a, b, Γ, u) + +# Finite-core velocity kernels +function bound_leg_velocity(a, b, Γ, ε) + na, nb, σ = norm(a), norm(b), dot(a, b) + term_1 = (na^2 - σ) / √(na^2 + ε^2) + (nb^2 - σ) / √(nb^2 + ε^2) + term_2 = a × b / (na^2 * nb^2 - σ^2 + ε^2 * (na^2 + nb^2 - 2 * na * nb)) + + Γ/4π * term_1 * term_2 +end + +trailing_leg_velocity(r, Γ, u, ε) = Γ/4π * normalize(r) × u / (norm(r) - dot(r, u) + ε^2 / (norm(r) + dot(r, u))) +trailing_legs_velocities(a, b, Γ, u, ε) = trailing_leg_velocity(a, Γ, u, ε) - trailing_leg_velocity(b, Γ, u, ε) +total_horseshoe_velocity(a, b, Γ, u, ε) = bound_leg_velocity(a, b, Γ, ε) + trailing_legs_velocities(a, b, Γ, u, ε) + +## Arrays of vortex lines +#==========================================================================================# + +abstract type AbstractVortex end + +## Horseshoe type +#==========================================================================================# + +""" + Horseshoe(r1, r2, rc, normal, chord) + +Define a horseshoe vortex with a start and endpoints ``r₁, r₂`` for the bound leg, a collocation point ``r``, a normal vector ``n̂``, and a finite core size. + +The finite core setup is not implemented for now. +""" +struct Horseshoe{T <: Real} <: AbstractVortex + r1 :: SVector{3,T} + r2 :: SVector{3,T} + rc :: SVector{3,T} + normal :: SVector{3,T} + core :: T +end + +Base.length(::Horseshoe) = 1 + +r1(hs :: Horseshoe) = hs.r1 +r2(hs :: Horseshoe) = hs.r2 + +function Horseshoe(r1, r2, rc, n, c) + T = promote_type(eltype(r1), eltype(r2), eltype(rc), eltype(n), eltype(c)) + Horseshoe{T}(r1, r2, rc, n, c) +end + +control_point(hs :: AbstractVortex) = hs.rc +normal_vector(hs :: AbstractVortex) = hs.normal + +r1(r, hs :: Horseshoe) = r - hs.r1 +r2(r, hs :: Horseshoe) = r - hs.r2 + +""" + transform(hs :: Horseshoe, T :: LinearMap) + +Generate a new `Horseshoe` with the points and normal vectors transformed by the `LinearMap` ``T``. +""" +transform(hs :: Horseshoe, T :: LinearMap) = setproperties(hs, + r1 = T(hs.r1), + r2 = T(hs.r2), + rc = T(hs.rc), + normal = T(hs.normal), +) + +transform(hs :: Horseshoe; rotation = I(3), translation = zeros(3)) = transform(hs, Translation(translation) ∘ LinearMap(rotation)) + +""" + bound_leg_center(hs :: Horseshoe) + +Compute the midpoint of the bound leg of a `Horseshoe`. +""" +bound_leg_center(hs :: Horseshoe) = (hs.r1 + hs.r2) / 2 + +""" + bound_leg_vector(hs :: Horseshoe) + +Compute the direction vector of the bound leg of a `Horseshoe`. +""" +bound_leg_vector(hs :: Horseshoe) = hs.r2 - hs.r1 + +""" + velocity(r, hs :: Horseshoe, Γ, u_hat = [1.,0.,0.]) + +Compute the induced velocity at a point ``r`` of a given `Horseshoe` `hs` with a bound leg of constant strength ``Γ`` and semi-infinite trailing legs pointing in a given direction ``û``, by default `û = x̂`. +""" +velocity(r, hs :: Horseshoe, Γ, V_hat = SVector{3,promote_type(eltype(r),eltype(Γ),eltype(hs.core))}(1,0,0)) = total_horseshoe_velocity(r - hs.r1, r - hs.r2, Γ, V_hat, hs.core) + +""" + bound_velocity(r, hs :: Horseshoes, Γ, u_hat) + +Compute the induced velocity at a point ``r`` from the bound leg with constant strength ``Γ`` of a given `Horseshoe` `hs`. +""" +bound_velocity(r, hs :: Horseshoe, Γ) = bound_leg_velocity(r - hs.r1, r - hs.r2, Γ, hs.core) + +""" + trailing_velocity(r, hs :: Horseshoes, Γ, u_hat) + +Compute the induced velocity at a point ``r`` from the semi-infinite trailing legs with constant strength ``Γ`` of a given `Horseshoe` `hs`. +""" +trailing_velocity(r, hs :: Horseshoe, Γ, V) = trailing_legs_velocities(r - hs.r1, r - hs.r2, Γ, V, hs.core) + +## Vortex ring type +#==========================================================================================# + +""" + VortexRing(r1, r2, r3, r4, r_c, n̂, ε) + +A vortex ring consisting of four points ``r_i, i = 1,…,4``, a collocation point ``r_c``, a normal vector ``n̂``, and a core size ``ε``. The following convention is adopted: + +``` + r1 —front leg→ r4 + | | +left leg right leg + ↓ ↓ + r2 —back leg-→ r3 +``` +""" +struct VortexRing{T <: Real} <: AbstractVortex + r1 :: SVector{3,T} + r2 :: SVector{3,T} + r3 :: SVector{3,T} + r4 :: SVector{3,T} + rc :: SVector{3,T} + normal :: SVector{3,T} + core :: T +end + +control_point(ring :: VortexRing) = ring.rc +normal_vector(ring :: VortexRing) = ring.normal + +Base.length(::VortexRing) = 1 + +""" +Computes the induced velocities at a point ``r`` of a VortexRing with constant strength ``Γ``. +""" +function velocity(r, ring :: VortexRing, Γ) + # Compute vectors to evaluation point + r1, r2, r3, r4 = r - ring.r1, r - ring.r2, r - ring.r3, r - ring.r4 + core = ring.core + + # Evaluate bound leg velocities and sum + v1 = bound_leg_velocity(r1, r4, Γ, core) + v2 = bound_leg_velocity(r4, r3, Γ, core) + v3 = bound_leg_velocity(r3, r2, Γ, core) + v4 = bound_leg_velocity(r2, r1, Γ, core) + + return v1 + v2 + v3 + v4 +end + +trailing_velocity(r, ring :: VortexRing, Γ, V) = trailing_legs_velocities(r - ring.r1, r - ring.r4, Γ, V, ring.core) + +""" + transform(ring :: VortexRing, T :: LinearMap) + +Generate a new `VortexRing` with the points and normal vectors transformed by the `LinearMap` ``T``. +""" +transform(ring :: VortexRing, T :: LinearMap) = setproperties(ring, + r1 = T(ring.r1), + r2 = T(ring.r2), + r3 = T(ring.r3), + r4 = T(ring.r4), + rc = T(ring.rc), + normal = T(ring.normal) +) + +bound_leg_center(ring :: VortexRing) = (ring.r1 + ring.r4) / 2 +bound_leg_vector(ring :: VortexRing) = ring.r4 - ring.r1 + +r1(ring :: VortexRing) = ring.r1 +r2(ring :: VortexRing) = ring.r4 \ No newline at end of file diff --git a/src/Aerodynamics/vlm_interface.jl b/src/Aerodynamics/vlm_interface.jl new file mode 100644 index 00000000..5fbae8a2 --- /dev/null +++ b/src/Aerodynamics/vlm_interface.jl @@ -0,0 +1,58 @@ +function quarter_point(p1, p2) + μ = SVector(1/4, 0, 1/4) + return @. (1 - μ) * p1 + μ * p2 +end + +function three_quarter_point(p1, p2) + μ = SVector(3/4, 0, 3/4) + return @. (1 - μ) * p1 + μ * p2 +end + +control_point(p1, p2, p3, p4) = (three_quarter_point(p1, p2) + three_quarter_point(p4, p3)) / 2 +bound_leg(p1, p2, p3, p4) = (quarter_point(p1, p2), quarter_point(p4, p3)) + +""" + bound_leg(panel :: Panel3D) + +Compute the bound leg for a `Panel3D`, for horseshoes/vortex rings. +""" +bound_leg(panel :: Panel3D) = bound_leg(panel.p1, panel.p2, panel.p3, panel.p4) + +""" + control_point(panel :: Panel3D) + +Compute the control point of a `Panel3D` for horseshoes/vortex rings, which is the 3-quarter point on each side in the ``x``-``z`` plane. +""" +control_point(panel :: Panel3D) = control_point(panel.p1, panel.p2, panel.p3, panel.p4) + + +""" + Horseshoe(panel :: Panel3D, normal, drift = zeros(3)) + +Generate a `Horseshoe` corresponding to a `Panel3D`, an associated normal vector, and a "drift velocity". +""" +function Horseshoe(panel :: Panel3D, normal, drift = SVector(0., 0., 0.); core_size = 0.) + r1, r2 = bound_leg(panel) + rc = control_point(panel) + drift + Horseshoe(r1, r2, rc, normal, core_size) +end + +""" +Constructor for vortex rings on a Panel3D using Lines. The following convention is adopted: + +``` + p1 —front leg→ p4 + | | +left leg right leg + ↓ ↓ + p2 —back leg-→ p3 +``` +""" +function VortexRing(panel :: Panel3D{T}, normal = normal_vector(panel), ε = 0.) where T <: Real + # r1 = quarter_point(panel.p1, panel.p2) + # r4 = quarter_point(panel.p4, panel.p3) + # r2 = normalize(panel.p2 - panel.p1) * 0.25 + r1 + # r3 = normalize(panel.p3 - panel.p4) * 0.25 + r2 + rc = control_point(panel) # (r1 + r2 + r3 + r4) / 4 + VortexRing{T}(panel.p1, panel.p2, panel.p3, panel.p4, rc, normal, ε) +end \ No newline at end of file From 0b37d4d96d2c32b34de1018af512840c2293ba40 Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 14:00:07 +0800 Subject: [PATCH 02/13] Updated docs, examples and plots, added benchmark --- .../3d_vortex_lattice_method/vlm_benchmark.jl | 62 ++++++++++++++++--- .../3d_vortex_lattice_method/vlm_wing.jl | 50 +++++---------- src/Aerodynamics/Cases/cases.jl | 32 +++++++--- src/Aerodynamics/VortexLattice/printing.jl | 12 ++-- .../AircraftGeometry/Wings/halfwing.jl | 51 ++++++++++++++- src/Tools/plot_tools.jl | 22 +++---- 6 files changed, 155 insertions(+), 74 deletions(-) diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl index 891f2086..20c419e3 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_benchmark.jl @@ -49,8 +49,7 @@ end ## @time CF, CM, CDiff, byu_sys = vlm_byu(); -## -# @benchmark vlm_byu() +print_coefficients([CF; CM], [CDiff,"—","—"], "BYU") ## AeroFuse.jl using AeroFuse @@ -83,10 +82,10 @@ function vlm_aerofuse() speed = 1., # m/s density = 1.225, # kg/m³ viscosity = 1.5e-5, # ??? - area = Sref, # m² - span = bref, # m - chord = cref, # m - location = rref # m + area = 30.0, # m² + span = 15.0, # m + chord = 2.0, # m + location = [0.50, 0.0, 0.0] # m ) ## Horseshoes @@ -103,8 +102,55 @@ end ## @time nfs, ffs, sys = vlm_aerofuse(); -## -print_coefficients(nfs,ffs) +print_coefficients(nfs, ffs, "AeroFuse") + +## Results +# ─────────────────────────────────── +# BYU Nearfield Farfield +# ─────────────────────────────────── +# CX 0.00246661 CDi 0.00247615 +# CY 0.0 CYff — +# CZ 0.244371 CL — +# Cl 0.0 — +# Cm -0.0208467 — +# Cn 0.0 — +# ─────────────────────────────────── + +# ──────────────────────────────────────── +# AeroFuse Nearfield Farfield +# ──────────────────────────────────────── +# CX 0.0023826 CDi 0.00249797 +# CY -0.0 CYff -0.0 +# CZ 0.245316 CL 0.245319 +# Cl 0.0 — +# Cm -0.0211525 — +# Cn 0.0 — +# ──────────────────────────────────────── + +## Benchmarking +@benchmark vlm_byu() + +# BenchmarkTools.Trial: 1938 samples with 1 evaluation. +# Range (min … max): 2.348 ms … 53.391 ms ┊ GC (min … max): 0.00% … 0.00% +# Time (median): 2.393 ms ┊ GC (median): 0.00% +# Time (mean ± σ): 2.578 ms ± 1.586 ms ┊ GC (mean ± σ): 3.06% ± 6.92% + +# █▃▁ +# ███▇█▆▅▅▅▅▃▃▃▄▁▁▃▃▄▄▁▃▁▄▃▁▃▃▃▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▁▁▁▁▃▁▁▁▃ █ +# 2.35 ms Histogram: log(frequency) by time 8.82 ms < + +# Memory estimate: 941.44 KiB, allocs estimate: 8188. ## @benchmark vlm_aerofuse() + +# BenchmarkTools.Trial: 4353 samples with 1 evaluation. +# Range (min … max): 1.037 ms … 14.134 ms ┊ GC (min … max): 0.00% … 91.51% +# Time (median): 1.085 ms ┊ GC (median): 0.00% +# Time (mean ± σ): 1.147 ms ± 666.701 μs ┊ GC (mean ± σ): 3.80% ± 6.19% + +# █▆▁ +# ▃▄▇███▇▄▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▁▂▂▂▂▁▁▂▁▁▂▂▁▁▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▂▁▁▂ ▃ +# 1.04 ms Histogram: frequency by time 1.67 ms < + +# Memory estimate: 587.05 KiB, allocs estimate: 1230. \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl index 382a56ef..dea01b31 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl @@ -4,13 +4,14 @@ using AeroFuse ## Wing section setup wing = Wing( foils = fill(naca4((2,4,1,2)), 3), - chords = [2.0, 1.6, 0.2], + chords = [2.0, 1.2, 0.6], twists = [0.0, 0.0, 0.0], spans = [5., 0.6], dihedrals = [5., 5.], - sweeps = [20.,20.], + sweeps = [20.,30.], w_sweep = 0.25, # Quarter-chord sweep symmetry = true, + position = [4.0, 0.0, 0.0] # flip = true ) @@ -19,15 +20,15 @@ x_w, y_w, z_w = wing_mac = mean_aerodynamic_center(wing) print_info(wing, "Wing") ## Meshing and assembly -wing_mesh = WingMesh(wing, [24,12], 6, - span_spacing = Cosine() +wing_mesh = WingMesh(wing, [24,6], 6, + # span_spacing = Cosine() ); aircraft = ComponentVector(wing = make_horseshoes(wing_mesh)) # Freestream conditions fs = Freestream( alpha = 2.0, # deg - beta = 2.0, # deg + beta = 0.0, # deg omega = [0.,0.,0.] ) @@ -67,7 +68,7 @@ ref = References( dvs = freestream_derivatives(system; axes = ax, print_components = true, - farfield = false + farfield = true ) end; @@ -92,16 +93,12 @@ CDi_ff, CY_ff, CL_ff = ff = farfield(system) nf_v = [ CDi_nf + CDv; CDv; nf ] ff_v = [ CDi_ff + CDv; CDv; ff ] +print_coefficients(nf_v, ff_v) + ## Plotting using Plots gr() -## Display -horseshoe_panels = chord_panels(wing_mesh) -horseshoe_coords = plot_panels(horseshoe_panels) -horseshoe_points = Tuple.(control_point.(system.vortices)) -ys = getindex.(horseshoe_points, 2) - ## Coordinates plot( aspect_ratio = 1, @@ -113,35 +110,17 @@ plot!(wing_mesh, label = "Wing") plot!(system, wing, dist = 3, num_stream = 50, span = 10, color = :green) ## Compute spanwise loads -CL_loads = vec(sum(system.circulations.wing, dims = 1)) / (0.5 * ref.speed * ref.chord) - -span_loads = spanwise_loading(wing_mesh, CFs.wing, ref.area) +span_loads = spanwise_loading(wing_mesh, ref, CFs.wing, system.circulations.wing) ## Plot spanwise loadings plot_CD = plot(span_loads[:,1], span_loads[:,2], label = :none, ylabel = "CDi") plot_CY = plot(span_loads[:,1], span_loads[:,3], label = :none, ylabel = "CY") plot_CL = begin plot(span_loads[:,1], span_loads[:,4], label = :none, xlabel = "y", ylabel = "CL") - plot!(span_loads[:,1], CL_loads, label = "Normalized", xlabel = "y") + plot!(span_loads[:,1], span_loads[:,5], label = "Normalized", xlabel = "y") end plot(plot_CD, plot_CY, plot_CL, size = (800, 700), layout = (3,1)) -## Lift distribution - -# Exaggerated CF distribution for plot -# hs_pts = vec(Tuple.(bound_leg_center.(system.vortices))) - -# plot(xaxis = "x", yaxis = "y", zaxis = "z", -# aspect_ratio = 1, -# camera = (60, 60), -# zlim = span(wing) .* (-0.5, 0.5), -# title = "Forces (Exaggerated)" -# ) -# plot!.(horseshoe_coords, color = :gray, label = :none) -# # scatter!(cz_pts, zcolor = vec(CLs), marker = 2, label = "CL (Exaggerated)") -# quiver!(hs_pts, quiver=(span_loads[:,2], span_loads[:,3], span_loads[:,4]) .* 10) -# plot!(size = (800, 600)) - ## VARIABLE ANALYSES #=========================================================# @@ -170,7 +149,7 @@ res_αs = combinedimsview( map(αs) do α fst = @set fs.alpha = deg2rad(α) sys = solve_case(aircraft, fst, ref) - [ α; farfield(sys)...; nearfield(sys)... ] + [ α; farfield(sys); nearfield(sys) ] end, (1) ) @@ -187,14 +166,13 @@ res = combinedimsview( refs = @set ref.speed = V fst = @set fs.alpha = deg2rad(α) sys = solve_case(aircraft, fst, refs) - [ mach_number(refs); α; farfield(sys)...; nearfield(sys)... ] + [ mach_number(refs); α; farfield(sys); nearfield(sys) ] end ) -## res_p = permutedims(res, (3,1,2)) -# CDi +## CDi plt_CDi_ff = plot(camera = (60,45)) [ plot!( res_p[:,1,n], res_p[:,2,n], res_p[:,3,n], diff --git a/src/Aerodynamics/Cases/cases.jl b/src/Aerodynamics/Cases/cases.jl index ca424b4d..14d1e737 100644 --- a/src/Aerodynamics/Cases/cases.jl +++ b/src/Aerodynamics/Cases/cases.jl @@ -77,11 +77,27 @@ end spanwise_loading(wing :: WingMesh, CFs, S) = spanwise_loading(chord_panels(wing), CFs, S) +""" + spanwise_loading(wing_mesh :: WingMesh, ref :: References, CFs, Γ) + +Obtain the spanwise aerodynamic loads `(CDi, CY, CL, CL_norm)` for a given `WingMesh`, reference values `ref`, surface coefficients `CFs`, and circulations ``Γ``. + +The chordwise normalized lift coefficient loading `CL_norm` is calculated via the formula ``C_L = 2Γ / ρVc`` based on the reference speed ``V`` and chord ``c`` from `ref`. + +Note that the surface coefficients `CFs` should be computed in wind axes to match `CL_norm`. +""" +function spanwise_loading(wing_mesh :: WingMesh, ref :: References, CFs, Γs) + span_loads = spanwise_loading(wing_mesh, CFs, ref.area) # + CL_loads = vec(sum(Γs, dims = 1)) / (0.5 * ref.speed * ref.chord) # Normalized CL loading + + [ span_loads CL_loads ] +end + ## Mesh connectivities -triangle_connectivities(inds) = @views [ - vec(inds[1:end-1,1:end-1]) vec(inds[1:end-1,2:end]) vec(inds[2:end,2:end]) ; - vec(inds[2:end,2:end]) vec(inds[2:end,1:end-1]) vec(inds[1:end-1,1:end-1]) - ] +@views triangle_connectivities(inds) = [ + vec(inds[1:end-1,1:end-1]) vec(inds[1:end-1,2:end]) vec(inds[2:end,2:end]) ; + vec(inds[2:end,2:end]) vec(inds[2:end,1:end-1]) vec(inds[1:end-1,1:end-1]) + ] ## Extrapolating surface values to neighbouring points function extrapolate_point_mesh(mesh, weight = 0.75) @@ -89,10 +105,10 @@ function extrapolate_point_mesh(mesh, weight = 0.75) points = zeros(eltype(mesh), m + 1, n + 1) # The quantities are measured at the bound leg of a Horseshoe by default (0.25×) - @views @. points[1:end-1,1:end-1] += weight / 2 * mesh - @views @. points[1:end-1,2:end] += weight / 2 * mesh - @views @. points[2:end,1:end-1] += (1 - weight) / 2 * mesh - @views @. points[2:end,2:end] += (1 - weight) / 2 * mesh + @. points[1:end-1,1:end-1] += weight / 2 * mesh + @. points[1:end-1,2:end] += weight / 2 * mesh + @. points[2:end,1:end-1] += (1 - weight) / 2 * mesh + @. points[2:end,2:end] += (1 - weight) / 2 * mesh points end \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/printing.jl b/src/Aerodynamics/VortexLattice/printing.jl index b072976f..acaeab8b 100644 --- a/src/Aerodynamics/VortexLattice/printing.jl +++ b/src/Aerodynamics/VortexLattice/printing.jl @@ -15,12 +15,12 @@ function Base.show(io :: IO, refs :: References) end end -function Base.show(io :: IO, ring :: AbstractVortex) - println(io, "Vortex: ") - for fname in fieldnames(typeof(ring)) - println(io, " ", fname, " = ", getfield(ring, fname)) - end -end +# function Base.show(io :: IO, ring :: AbstractVortex) +# println(io, "Vortex: ") +# for fname in fieldnames(typeof(ring)) +# println(io, " ", fname, " = ", getfield(ring, fname)) +# end +# end # Vortex lattice system function Base.show(io :: IO, sys :: VortexLatticeSystem) diff --git a/src/Geometry/AircraftGeometry/Wings/halfwing.jl b/src/Geometry/AircraftGeometry/Wings/halfwing.jl index 69e1ce2c..2462a0ce 100644 --- a/src/Geometry/AircraftGeometry/Wings/halfwing.jl +++ b/src/Geometry/AircraftGeometry/Wings/halfwing.jl @@ -112,11 +112,21 @@ affine_transformation(wing :: Wing) = wing.affine (f :: AffineMap)(wing :: Wing) = @set wing.affine = f ∘ wing.affine +""" + span(wing :: Wing) + +Compute the planform span of a `Wing`. +""" function span(wing :: Wing) b = sum(wing.spans) ifelse(wing.symmetry, 2b, b) end +""" + taper_ratio(wing :: Wing) + +Compute the taper ratio of a `Wing`, defined as the tip chord length divided by the root chord length, independent of the number of sections. +""" taper_ratio(wing :: Wing) = last(wing.chords) / first(wing.chords) section_projected_area(b, c1, c2, t1, t2) = b * (c1 + c2) / 2 * cosd((t1 + t2) / 2) @@ -126,16 +136,35 @@ function section_projected_areas(wing :: Wing) @views section_projected_area.(spans, wing.chords[1:end-1], wing.chords[2:end], wing.twists[1:end-1], wing.twists[2:end]) end +""" + projected_area(wing :: Wing) + +Compute the projected area (onto the spanwise plane) of a `Wing`. +""" projected_area(wing :: Wing) = sum(section_projected_areas(wing)) section_macs(wing :: Wing) = @views @. mean_aerodynamic_chord(wing.chords[1:end-1], wing.chords[2:end] / wing.chords[1:end-1]) +""" + mean_aerodynamic_chord(wing :: Wing) + +Compute the mean aerodynamic chord of a `Wing`. +""" function mean_aerodynamic_chord(wing :: Wing) areas = section_projected_areas(wing) macs = section_macs(wing) sum(macs .* areas) / sum(areas) end +""" + mean_aerodynamic_center(wing :: Wing, + factor = 0.25; + symmetry = wing.symmetry, + flip = wing.flip + ) + +Compute the mean aerodynamic center of a `Wing`. By default, the factor is assumed to be at 25% from the leading edge, which can be adjusted. Similarly, options are provided to account for symmetry or to flip the location in the ``x``-``z`` plane. +""" function mean_aerodynamic_center(wing :: Wing, factor = 0.25; symmetry = wing.symmetry, flip = wing.flip) # Compute mean aerodynamic chords and projected areas macs = section_macs(wing) @@ -167,8 +196,18 @@ function mean_aerodynamic_center(wing :: Wing, factor = 0.25; symmetry = wing.sy return affine_transformation(wing)(mac) end +""" + camber_thickness(wing :: Wing, num :: Integer) + +Compute the camber-thickness distribution at each spanwise intersection of a `Wing`. A `num` must be specified to interpolate the internal `Foil` coordinates, which affects the accuracy of ``(t/c)ₘₐₓ`` accordingly. +""" camber_thickness(wing :: Wing, num :: Integer) = camber_thickness.(wing.foils, num) +""" + maximum_thickness_to_chord(wing :: Wing, num :: Integer) + +Compute the maximum thickness-to-chord ratio ``(t/c)ₘₐₓ`` at each spanwise intersection of a `Wing`. A `num` must be specified to interpolate the internal `Foil` coordinates, which affects the accuracy of ``(t/c)ₘₐₓ`` accordingly. +""" maximum_thickness_to_chord(wing :: Wing, num :: Integer) = map(maximum_thickness_to_chord, camber_thickness(wing, num)) function max_thickness_to_chord_ratio_sweeps(wing :: Wing, num) @@ -181,7 +220,7 @@ function max_thickness_to_chord_ratio_sweeps(wing :: Wing, num) le = leading_edge(wing) # Determine x-coordinates in geometry frame - xs = @. getindex(le, 1) + wing.chords * getindex(xs_temp, 1) + xs = @. le[:,1] + wing.chords * getindex(xs_temp, 1) # Compute sectional x-coordinates ds = xs[2:end] - xs[1:end-1] @@ -222,11 +261,17 @@ function wing_bounds(wing :: Wing) end """ - trailing_edge(wing :: Wing, flip = false) + trailing_edge(wing :: Wing) -Compute the trailing edge coordinates of a `Wing`, with an option to flip the signs of the ``y``-coordinates. +Compute the leading edge coordinates of a `Wing`. """ leading_edge(wing :: Wing) = @views wing_bounds(wing)[:,:,1] + +""" + trailing_edge(wing :: Wing) + +Compute the trailing edge coordinates of a `Wing`. +""" trailing_edge(wing :: Wing) = @views wing_bounds(wing)[:,:,2] diff --git a/src/Tools/plot_tools.jl b/src/Tools/plot_tools.jl index 71ab72f2..4465cebf 100644 --- a/src/Tools/plot_tools.jl +++ b/src/Tools/plot_tools.jl @@ -83,7 +83,7 @@ end @series begin ls := :dash label := foil.name * " Camber" - xcam[:,1], xcam[:,2] + @views xcam[:,1], xcam[:,2] end end if thickness @@ -91,7 +91,7 @@ end @series begin ls := :dot label := foil.name * " Thickness" - xthicc[:,1], xthicc[:,2] + @views xthicc[:,1], xthicc[:,2] end end @@ -107,12 +107,12 @@ end # zlim --> span(wing) .* (-0.5, 0.5) @series begin - wing_plan[:,1], wing_plan[:,2], wing_plan[:,3] + @views wing_plan[:,1], wing_plan[:,2], wing_plan[:,3] end @series begin seriestype := :scatter - [wing_mac[1]], [wing_mac[2]], [wing_mac[3]] + @views [wing_mac[1]], [wing_mac[2]], [wing_mac[3]] end end @@ -125,7 +125,7 @@ end primary := false # linecolor := :lightgray # fillcolor := :lightgray - coords[:,1], coords[:,2], coords[:,3] + @views coords[:,1], coords[:,2], coords[:,3] end end end @@ -143,14 +143,14 @@ end primary := false linecolor := :lightgray # fillcolor := :lightgray - coords[:,1], coords[:,2], coords[:,3] + @views coords[:,1], coords[:,2], coords[:,3] end end @series begin seriestype := :scatter # label --> "Mean Aerodynamic Chord" - [wing_mac[1]], [wing_mac[2]], [wing_mac[3]] + @views [wing_mac[1]], [wing_mac[2]], [wing_mac[3]] end @series begin @@ -199,16 +199,12 @@ end end end -# get_span_points(wing :: Wing, pts) = (wing.right.affine).(chop_leading_edge(wing, pts)) -get_span_points(wing :: Wing, pts) = affine_transformation(wing).(chop_leading_edge(wing, pts)) - - @recipe function streamline_plot(system :: VortexLatticeSystem, wing :: AbstractWing; dist = 5 * mean_aerodynamic_chord(wing), num = 100, span = 20, linecolor = :green) # ys = LinRange(-span(wing) / 2, span(wing) / 2, span_points) # init = SVector.(0., ys, -0.5) - init = get_span_points(wing, span) dx, dy, dz = 0, 0, 1e-3 - seed = init .+ Ref([dx, dy, dz]); + init = chop_leading_edge(wing, span) + seed = init .+ Ref([dx, dy, dz]); streams = plot_streamlines(system, seed, dist, num) From 34a4a9b866915a6f5a199a1f6a05a60487c1be6f Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 17:19:34 +0800 Subject: [PATCH 03/13] Vortex ring formulation works with correct results --- .../vlm_vortex_rings.jl | 21 +++----- src/AeroFuse.jl | 51 ++----------------- src/Aerodynamics/VortexLattice/vortices.jl | 25 +++++++-- src/Aerodynamics/vlm_interface.jl | 42 +++++++++++++-- 4 files changed, 69 insertions(+), 70 deletions(-) diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl index 0fe17ce4..0c7e10ea 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl @@ -5,7 +5,7 @@ using StructArrays ## Wing section setup wing = Wing( - foils = fill(naca4((0,0,1,2)), 2), + foils = fill(naca4((2,4,1,2)), 2), chords = [2.2, 1.8], twists = [0.0, 0.0], spans = [7.5], @@ -15,11 +15,11 @@ wing = Wing( symmetry = true, ) -wing_mesh = WingMesh(wing, [24], 6, span_spacing = Uniform()); +wing_mesh = WingMesh(wing, [24], 8, span_spacing = Uniform()); # Freestream conditions fs = Freestream( - alpha = 1.0, # deg + alpha = 2.0, # deg beta = 0.0, # deg omega = [0.,0.,0.] ) @@ -32,24 +32,15 @@ ref = References( area = projected_area(wing), # m² span = span(wing), # m chord = mean_aerodynamic_chord(wing), # m - location = [0.5,0.,0.] # m + location = mean_aerodynamic_center(wing) # m ) ## Horseshoes ac_hs = ComponentVector(wing = make_horseshoes(wing_mesh)) system = VortexLatticeSystem(ac_hs, fs, ref) -print_coefficients(nearfield(system), farfield(system)) +print_coefficients(nearfield(system ), farfield(system)) ## Vortex rings ac_vs = ComponentVector(wing = make_vortex_rings(wing_mesh)) sys = VortexLatticeSystem(ac_vs, fs, ref) -print_coefficients(nearfield(sys), farfield(sys)) - -## Manual -A = influence_matrix(ac_vs) -b = boundary_condition(ac_vs, -velocity(fs), fs.omega) -gam = A \ b - -## -sys = VortexLatticeSystem(ac_vs, gam * ref.speed, A, b, fs, ref, Wind()) -print_coefficients(nearfield(sys), farfield(sys)) +print_coefficients(nearfield(sys), farfield(sys)) \ No newline at end of file diff --git a/src/AeroFuse.jl b/src/AeroFuse.jl index 5e0ed7bb..1d05dec1 100644 --- a/src/AeroFuse.jl +++ b/src/AeroFuse.jl @@ -125,52 +125,6 @@ export Wing, WingSection, affine_transformation, mean_aerodynamic_chord, span, a # export WingControlSurface -make_horseshoes(wing :: WingMesh) = map((cho,cam) -> Horseshoe(cho, normal_vector(cam)), chord_panels(wing), camber_panels(wing)) - -function make_vortex_rings(wing :: WingMesh) - cam_coo = camber_coordinates(wing) - cams = combinedimsview(cam_coo, (1,2)) - - # Generate vortex ring mesh - vor_cams = similar(cams) - vor_cams = 0.75 * cams[1:end-1,:,:] + 0.25 * cams[2:end,:,:] - vor_cams[end,:,:] = cams[end,:,:] - - @views rings = VortexRing.(make_panels(splitdimsview(vor_cams, (1,2)))) - - # NOTE: JUST ADD A TRAILING TYPE TO THE VORTEX RING AND ADD THE HORSESHOE VELOCITY METHOD - - # Trailing edge semi-infinite vortex lines - cam_pan = camber_panels(wing) - @views horsies = Horseshoe.(cam_pan[end,:], normal_vector.(cam_pan[end,:])) - - # Generate vortex ring mesh - # vor_cams = similar(cams) - # vor_cams[1:end-1,:,:] = 0.75 * cams[1:end-1,:,:] + 0.25 * cams[2:end,:,:] - # vor_cams[end,:,:] = cams[end,:,:] - - # # display(vor_cams) - - # vor_coo = splitdimsview(vor_cams, (1,2)) - # vor_pan = make_panels(vor_coo) - # @views rings = @. VortexRing(vor_pan) - - # # Trailing edge semi-infinite vortex lines - # hor_cams = similar(cams[end-1:end,:,:]) - # hor_cams[end-1,:,:] = vor_cams[end,:,:] - # hor_cams[end,:,:] = (cams[end,:,:] - cams[end-1,:,:]) * 0.25 + cams[end,:,:] - - # hor_coo = splitdimsview(hor_cams, (1,2)) - # hor_pan = make_panels(hor_coo) - # @views horsies = @. Horseshoe(hor_pan[end,:], normal_vector(hor_pan[end,:])) - - - # Assemble - [ rings; permutedims(horsies) ] -end - -export make_horseshoes, make_vortex_rings - ## Aerodynamic analyses #==========================================================================================# @@ -198,13 +152,14 @@ export Horseshoe, VortexLatticeSystem, References, AbstractAxisSystem, Stability ## Panel-VLM interface include("Aerodynamics/vlm_interface.jl") +export make_horseshoes, make_vortex_rings + ## Profile drag estimation include("Aerodynamics/profile_drag.jl") export profile_drag_coefficient, wetted_area_drag, local_dissipation_drag, form_factor ## Cases - include("Aerodynamics/Cases/printing.jl") export print_info @@ -247,6 +202,6 @@ export induced_velocity, induced_speed, inflow_angle, blade_solidity, slipstream include("Tools/plot_tools.jl") -export plot_panel, plot_panels, plot_streamlines, plot_planform, plot_surface +export plot_panel, plot_panels, plot_streamlines, plot_planform, plot_surface, plot_spanload end \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/vortices.jl b/src/Aerodynamics/VortexLattice/vortices.jl index c7424052..f193a80c 100644 --- a/src/Aerodynamics/VortexLattice/vortices.jl +++ b/src/Aerodynamics/VortexLattice/vortices.jl @@ -131,18 +131,26 @@ struct VortexRing{T <: Real} <: AbstractVortex r4 :: SVector{3,T} rc :: SVector{3,T} normal :: SVector{3,T} + trailing :: Bool core :: T end +function VortexRing(r1, r2, r3, r4, rc, n, trailing, c) + T = promote_type(eltype(r1), eltype(r2), eltype(rc), eltype(n), eltype(c)) + VortexRing{T}(r1, r2, r3, r4, rc, n, trailing, c) +end + control_point(ring :: VortexRing) = ring.rc normal_vector(ring :: VortexRing) = ring.normal Base.length(::VortexRing) = 1 """ -Computes the induced velocities at a point ``r`` of a VortexRing with constant strength ``Γ``. + velocity(r, ring :: VortexRing, Γ) + +Computes the velocity at a point ``r`` induced by a `VortexRing` with constant strength ``Γ``. """ -function velocity(r, ring :: VortexRing, Γ) +function velocity(r, ring :: VortexRing, Γ, V_hat = SVector(1,0,0)) # Compute vectors to evaluation point r1, r2, r3, r4 = r - ring.r1, r - ring.r2, r - ring.r3, r - ring.r4 core = ring.core @@ -153,7 +161,14 @@ function velocity(r, ring :: VortexRing, Γ) v3 = bound_leg_velocity(r3, r2, Γ, core) v4 = bound_leg_velocity(r2, r1, Γ, core) - return v1 + v2 + v3 + v4 + v = v1 + v2 + v3 + v4 + + # Add horseshoe velocity contribution if trailing edge panel + if ring.trailing + v += total_horseshoe_velocity(r2, r3, Γ, V_hat, ring.core) + end + + return v end trailing_velocity(r, ring :: VortexRing, Γ, V) = trailing_legs_velocities(r - ring.r1, r - ring.r4, Γ, V, ring.core) @@ -169,7 +184,9 @@ transform(ring :: VortexRing, T :: LinearMap) = setproperties(ring, r3 = T(ring.r3), r4 = T(ring.r4), rc = T(ring.rc), - normal = T(ring.normal) + normal = T(ring.normal), + trailing = ring.trailing, + core = ring.core, ) bound_leg_center(ring :: VortexRing) = (ring.r1 + ring.r4) / 2 diff --git a/src/Aerodynamics/vlm_interface.jl b/src/Aerodynamics/vlm_interface.jl index 5fbae8a2..b1ffa9de 100644 --- a/src/Aerodynamics/vlm_interface.jl +++ b/src/Aerodynamics/vlm_interface.jl @@ -48,11 +48,47 @@ left leg right leg p2 —back leg-→ p3 ``` """ -function VortexRing(panel :: Panel3D{T}, normal = normal_vector(panel), ε = 0.) where T <: Real +function VortexRing(panel :: Panel3D{T}, rc, normal, trailing = false; core_size = 0.) where T <: Real # r1 = quarter_point(panel.p1, panel.p2) # r4 = quarter_point(panel.p4, panel.p3) # r2 = normalize(panel.p2 - panel.p1) * 0.25 + r1 # r3 = normalize(panel.p3 - panel.p4) * 0.25 + r2 - rc = control_point(panel) # (r1 + r2 + r3 + r4) / 4 - VortexRing{T}(panel.p1, panel.p2, panel.p3, panel.p4, rc, normal, ε) + # rc = control_point(panel) # (r1 + r2 + r3 + r4) / 4 + VortexRing{T}(panel.p1, panel.p2, panel.p3, panel.p4, rc, normal, trailing, core_size) +end + + +""" + make_horseshoes(wing :: WingMesh) + +Generate an array of `Horseshoe`s defined by the chord coordinates and normal vectors defined by the camber distribution of a `WingMesh`. +""" +make_horseshoes(wing :: WingMesh) = map((cho,cam) -> Horseshoe(cho, normal_vector(cam)), chord_panels(wing), camber_panels(wing)) + +""" + make_vortex_rings(wing :: WingMesh) + +Generate an array of `VortexRing`s defined by the camber coordinates and normal vectors of a `WingMesh`. +""" +@views function make_vortex_rings(wing_mesh :: WingMesh) + cam_pan = camber_panels(wing_mesh) + cams = combinedimsview(camber_coordinates(wing_mesh), (1,2)) + + # Generate vortex ring mesh + vor_cams = similar(cams) + vor_cams[1:end-1,:,:] = 0.75 * cams[1:end-1,:,:] + 0.25 * cams[2:end,:,:] + vor_cams[end,:,:] = cams[end,:,:] + + # Construct vortex rings with trailing edge identification for boundary condition + vor_pans = make_panels(splitdimsview(vor_cams, (1,2))) + rings = map(CartesianIndices(vor_pans)) do ind + i, j = ind.I + if i == size(vor_pans, 2) + VortexRing(vor_pans[i,j], control_point(cam_pan[i,j]), normal_vector(cam_pan[i,j])) + else + VortexRing(vor_pans[i,j], control_point(cam_pan[i,j]), normal_vector(cam_pan[i,j]), true) + end + end + + return rings end \ No newline at end of file From 876aaf8cc0c44d6dcd961801a5edb72ffb7f9042 Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 17:23:21 +0800 Subject: [PATCH 04/13] VLM compressibility corrections are now optional --- src/Aerodynamics/Cases/cases.jl | 5 ++-- src/Aerodynamics/VortexLattice/stability.jl | 21 +++++++-------- src/Aerodynamics/VortexLattice/system.jl | 30 ++++++++++++++++----- 3 files changed, 36 insertions(+), 20 deletions(-) diff --git a/src/Aerodynamics/Cases/cases.jl b/src/Aerodynamics/Cases/cases.jl index 14d1e737..2b760408 100644 --- a/src/Aerodynamics/Cases/cases.jl +++ b/src/Aerodynamics/Cases/cases.jl @@ -7,14 +7,15 @@ fs :: Freestream, refs :: References; name = :aircraft :: Symbol, + compressible = false :: Boolean, print = false :: Boolean, print_components = false :: Boolean ) Perform a vortex lattice analysis given a vector of `Horseshoe`s, a `Freestream` condition, and `References` values. """ -function solve_case(components :: DenseArray, freestream :: Freestream, refs :: References; name = :aircraft, print = false, print_components = false) - system = VortexLatticeSystem(components, freestream, refs) +function solve_case(components :: DenseArray, freestream :: Freestream, refs :: References; name = :aircraft, compressible = false, print = false, print_components = false) + system = VortexLatticeSystem(components, freestream, refs, compressible) # Printing if needed if print_components diff --git a/src/Aerodynamics/VortexLattice/stability.jl b/src/Aerodynamics/VortexLattice/stability.jl index e3c84e35..cee220f2 100644 --- a/src/Aerodynamics/VortexLattice/stability.jl +++ b/src/Aerodynamics/VortexLattice/stability.jl @@ -26,7 +26,7 @@ scale_freestream(fs :: Freestream, refs :: References) = ] # Closure to generate results with input vector -function freestream_derivatives!(y, x, aircraft, fs, ref, axes = Wind()) +function freestream_derivatives!(y, x, aircraft, fs, ref, compressible = false, axes = Wind()) # Scale Mach number ref = @views setproperties(ref, speed = x[1] * ref.sound_speed @@ -40,17 +40,14 @@ function freestream_derivatives!(y, x, aircraft, fs, ref, axes = Wind()) ) # Solve system - system = VortexLatticeSystem(aircraft, fs, ref) + system = VortexLatticeSystem(aircraft, fs, ref, compressible) # Evaluate nearfield coefficients CFs, CMs = surface_coefficients(system; axes = axes) # Evaluate farfield coefficients - q = dynamic_pressure(system.reference.density, system.reference.speed) - FFs = map(farfield_forces(system)) do F - force_coefficient(F, q, system.reference.area) - end - + FFs = farfield_coefficients(system) + # Create array of nearfield and farfield coefficients as a column, for each component as a row. comp_coeffs = @views mapreduce(name -> [ sum(CFs[name]); sum(CMs[name]); FFs[name] ], hcat, keys(system.vortices)) @@ -63,19 +60,19 @@ function freestream_derivatives!(y, x, aircraft, fs, ref, axes = Wind()) y[:,end] = comp_coeffs end - # return [ comp_coeffs sum_coeffs ] # Append sum of all forces for aircraft + # return nothing end -function freestream_derivatives(aircraft, fs, ref; axes = Wind(), name = :aircraft, print = false, print_components = false, farfield = false) +function freestream_derivatives(aircraft, fs, ref; axes = Wind(), name = :aircraft, compressible = false, print = false, print_components = false, farfield = false) # Reference values and scaling inputs x = scale_freestream(fs, ref) - names = [ reduce(vcat, keys(aircraft)); name ] + names = [ reduce(vcat, keys(aircraft)); name ] num_comps = length(names) y = zeros(eltype(x), 9, num_comps) ∂y∂x = zeros(eltype(x), 9, num_comps, length(x)) - jacobian!(∂y∂x, (y, x) -> freestream_derivatives!(y, x, aircraft, fs, ref, axes), y, x) + jacobian!(∂y∂x, (y, x) -> freestream_derivatives!(y, x, aircraft, fs, ref, compressible, axes), y, x) data = cat(y, ∂y∂x, dims = 3) # Appending outputs with derivatives @@ -103,4 +100,4 @@ end Obtain the force and moment coefficients of a `VortexLatticeSystem` and the derivatives of its components with respect to freestream values. """ -freestream_derivatives(system :: VortexLatticeSystem; axes = Wind(), name = :aircraft, print = false, print_components = false, farfield = false) = freestream_derivatives(system.vortices, system.freestream, system.reference; axes = axes, name = name, print = print, print_components = print_components, farfield = farfield) \ No newline at end of file +freestream_derivatives(system :: VortexLatticeSystem; axes = Wind(), name = :aircraft, print = false, print_components = false, farfield = false) = freestream_derivatives(system.vortices, system.freestream, system.reference; axes = axes, name = name, compressible = system.compressible, print = print, print_components = print_components, farfield = farfield) \ No newline at end of file diff --git a/src/Aerodynamics/VortexLattice/system.jl b/src/Aerodynamics/VortexLattice/system.jl index 2acf3fd5..23e4080d 100644 --- a/src/Aerodynamics/VortexLattice/system.jl +++ b/src/Aerodynamics/VortexLattice/system.jl @@ -85,28 +85,46 @@ struct VortexLatticeSystem{ boundary_vector :: S freestream :: P reference :: Q - axes :: T + compressible :: Bool + axes :: T end -function VortexLatticeSystem(components, fs :: Freestream, refs :: References, axes = Geometry()) +""" + VortexLatticeSystem( + aircraft, + fs :: Freestream, + refs :: References, + compressible = false, + axes = Geometry() + ) + +Construct a `VortexLatticeSystem` for analyzing inviscid aerodynamics of an aircraft (must be a `ComponentArray` of `Horseshoe`s or `VortexRing`s) with `Freestream` conditions and `References` for non-dimensionalization. Options are provided for compressibility corrections via the Prandtl-Glauert transformation (false by default) and axis system for computing velocities and forces (`Geometry` by default). +""" +function VortexLatticeSystem(aircraft, fs :: Freestream, refs :: References, compressible = false, axes = Geometry()) # Mach number bound checks M = mach_number(refs) @assert M < 1. "Only compressible subsonic flow conditions (M < 1) are valid!" if M > 0.7 @warn "Results in transonic to sonic flow conditions (0.7 < M < 1) are most likely incorrect!" end + if M > 0.3 && !compressible @warn "Compressible regime (M > 0.3) but compressibility correction is off, be wary of the analysis!" end # (Prandtl-Glauert ∘ Wind axis) transformation - β_pg = √(1 - M^2) - comp = @. prandtl_glauert_scale_coordinates(geometry_to_wind_axes(components, fs), β_pg) + if compressible + β_pg = √(1 - M^2) + ac = @. prandtl_glauert_scale_coordinates(geometry_to_wind_axes(aircraft, fs), β_pg) + else + β_pg = 1 + ac = @. geometry_to_wind_axes(aircraft, fs) + end # Quasi-steady freestream velocity U = geometry_to_wind_axes(velocity(fs, Body()), fs) Ω = geometry_to_wind_axes(fs.omega, fs) / refs.speed # Solve system - Γs, AIC, boco = solve_linear(comp, U, Ω) + Γs, AIC, boco = solve_linear(ac, U, Ω) - return VortexLatticeSystem(components, refs.speed * Γs / β_pg^2, AIC, boco, fs, refs, axes) + return VortexLatticeSystem(aircraft, refs.speed * Γs / β_pg^2, AIC, boco, fs, refs, compressible, axes) end # Miscellaneous From db0774f811d894a3764e56d7826295a3aaae8922 Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 17:23:35 +0800 Subject: [PATCH 05/13] Set default cosine distribution for single section wing --- src/Geometry/AircraftGeometry/Wings/mesh_wing.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/src/Geometry/AircraftGeometry/Wings/mesh_wing.jl b/src/Geometry/AircraftGeometry/Wings/mesh_wing.jl index af2347f8..ca300019 100644 --- a/src/Geometry/AircraftGeometry/Wings/mesh_wing.jl +++ b/src/Geometry/AircraftGeometry/Wings/mesh_wing.jl @@ -96,8 +96,12 @@ end # Spacing function symmetric_spacing(wing :: Wing) sins = AbstractSpacing[Sine(1); Sine(0)] - if length(wing.spans) == 1 && wing.symmetry - return sins + if length(wing.spans) == 1 + if wing.symmetry + return sins + else + return Cosine() + end else cosins = fill(Cosine(), (length(spans(wing)) - 1)) return [ cosins; sins; cosins ] @@ -169,7 +173,7 @@ end WingMesh(surface, n_span :: Integer, n_chord :: Integer; chord_spacing = Cosine(), span_spacing = symmetric_spacing(surface)) = WingMesh(surface, number_of_spanwise_panels(surface, n_span), n_chord, chord_spacing, span_spacing) # Forwarding functions for Wing type -MacroTools.@forward WingMesh.surface foils, chords, spans, twists, sweeps, dihedrals, mean_aerodynamic_chord, camber_thickness, leading_edge, trailing_edge, wing_bounds, mean_aerodynamic_center, projected_area, span, position, orientation, affine_transformation, maximum_thickness_to_chord +MacroTools.@forward WingMesh.surface foils, chords, spans, twists, sweeps, dihedrals, mean_aerodynamic_chord, camber_thickness, leading_edge, trailing_edge, wing_bounds, mean_aerodynamic_center, projected_area, span, position, orientation, affine_transformation, maximum_thickness_to_chord, chop_leading_edge """ chord_coordinates(wing :: WingMesh, n_span = wing.num_span, n_chord = wing.num_chord) From 1febb7dcd8ae1333bbff3601d68c2490283adc06 Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 18:53:32 +0800 Subject: [PATCH 06/13] Updated runtests with fixed Freestream sign --- src/Tools/Laplace.jl | 4 +- test/aircraft_definition.jl | 58 +++++++++ test/runtests.jl | 245 ++++++++++++------------------------ 3 files changed, 143 insertions(+), 164 deletions(-) create mode 100644 test/aircraft_definition.jl diff --git a/src/Tools/Laplace.jl b/src/Tools/Laplace.jl index 39e7ff80..7aeb821c 100644 --- a/src/Tools/Laplace.jl +++ b/src/Tools/Laplace.jl @@ -252,7 +252,7 @@ Freestream(α_deg, β_deg, Ω) = let T = promote_type(eltype(α_deg), eltype(β_ Freestream(α_deg, β_deg, Ω_x, Ω_y, Ω_z) = Freestream(deg2rad(α_deg), deg2rad(β_deg), SVector(Ω_x, Ω_y, Ω_z)) -Freestream(U, Ω) = let (_, α, β) = cartesian_to_freestream(normalize(U)); Freestream(α, β, Ω) end +Freestream(U, Ω) = let (_, α, β) = cartesian_to_freestream(normalize(U)); Freestream(α, -β, Ω) end Freestream(; alpha = 0., beta = 0., omega = zeros(3)) = Freestream(alpha, beta, omega) @@ -270,7 +270,7 @@ potential(fs :: Freestream, r) = dot(velocity(fs), r) Compute the velocity of a `Freestream`. """ -velocity(fs :: Freestream) = freestream_to_cartesian(1., fs.alpha, fs.beta) +velocity(fs :: Freestream) = freestream_to_cartesian(one(promote_type(eltype(fs.alpha), eltype(fs.beta))), fs.alpha, fs.beta) """ freestream_to_cartesian(r, θ, φ) diff --git a/test/aircraft_definition.jl b/test/aircraft_definition.jl new file mode 100644 index 00000000..245deeb8 --- /dev/null +++ b/test/aircraft_definition.jl @@ -0,0 +1,58 @@ +# Wing +wing = Wing( + foils = fill(naca4((0,0,1,2)), 2), + chords = [1.0, 0.6], + twists = [0.0, 0.0], + spans = [5.0] / 2, + dihedrals = [11.39], + sweeps = [0.], + symmetry = true, +); + +# Horizontal tail +htail = Wing( + foils = fill(naca4((0,0,1,2)), 2), + chords = [0.7, 0.42], + twists = [0.0, 0.0], + spans = [1.25] / 2, + dihedrals = [0.], + sweeps = [6.39], + position = [4., 0, 0], + angle = -2., + axis = [0., 1., 0.], + symmetry = true +) + +# Vertical tail +vtail = Wing( + foils = fill(naca4((0,0,0,9)), 2), + chords = [0.7, 0.42], + twists = [0.0, 0.0], + spans = [1.0], + dihedrals = [0.], + sweeps = [7.97], + position = [4., 0, 0], + angle = 90., + axis = [1., 0., 0.], +) + +# Assembly +wing_mesh = WingMesh(wing, [32], 10; span_spacing = Cosine()) +htail_mesh = WingMesh(htail, [12], 6; span_spacing = Cosine()) +vtail_mesh = WingMesh(vtail, [5], 6; span_spacing = Cosine()) + +## Reference quantities +fs = Freestream( + alpha = 1.0, + beta = 1.0, + omega = zeros(3) +) + +refs = References( + speed = 150.0, + area = projected_area(wing), + span = span(wing), + chord = mean_aerodynamic_chord(wing), + density = 1.225, + location = [0.25 * mean_aerodynamic_chord(wing), 0., 0.] +) diff --git a/test/runtests.jl b/test/runtests.jl index c175c47f..1ab736e6 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -104,7 +104,7 @@ end @testset "Geometry - Two-Section Trapezoidal Wing" begin # Define wing - wing_right = Wing( + wing = Wing( chords = [1.0, 0.6, 0.2], twists = [2.0, 0.0, -0.2], spans = [5.0, 0.5], @@ -114,12 +114,12 @@ end ); # Get wing info - b = span(wing_right) - S = projected_area(wing_right) - c = mean_aerodynamic_chord(wing_right) - AR = aspect_ratio(wing_right) - λ = taper_ratio(wing_right) - wing_mac = mean_aerodynamic_center(wing_right) + b = span(wing) + S = projected_area(wing) + c = mean_aerodynamic_chord(wing) + AR = aspect_ratio(wing) + λ = taper_ratio(wing) + wing_mac = mean_aerodynamic_center(wing) @test b ≈ 5.50000000 atol = 1e-6 @test S ≈ 4.19939047 atol = 1e-6 @@ -129,11 +129,34 @@ end @test wing_mac ≈ [0.4209310, 1.3343524, 0.0] atol = 1e-6 end +@testset "Geometry - Single-Section Trapezoidal Wing" begin + # Define wing section + wing_sec = WingSection(aspect = 6.25, area = 4.0, taper = 0.6, symmetry = true) + + # Using Wing constructor + wing = Wing( + foils = fill(naca4((0,0,1,2)), 2), + chords = [1.0, 0.6], + twists = [0.0, 0.0], + spans = [5.0] / 2, + dihedrals = [11.39], + sweeps = [0.], + symmetry = true, + ); + + @test span(wing) ≈ span(wing_sec) atol = 1e-6 # Span + @test projected_area(wing) ≈ projected_area(wing_sec) atol = 1e-6 # Projected area + @test mean_aerodynamic_chord(wing) ≈ mean_aerodynamic_chord(wing_sec) atol = 1e-6 # Mean aerodynamic chord + @test aspect_ratio(wing) ≈ aspect_ratio(wing_sec) atol = 1e-6 # Aspect ratio + @test taper_ratio(wing) ≈ taper_ratio(wing_sec) atol = 1e-6 # Taper ratio + @test mean_aerodynamic_center(wing) ≈ mean_aerodynamic_center(wing_sec) atol = 1e-6 # Mean aerodynamic center +end + @testset "Geometry - 3D Panel" begin panel = Panel3D( [1.0, -1., 0.0], - [0.0, -1., -0.5], - [0.0, 0.0, -0.5], + [0.0, -1., -0.0], + [0.0, 0.0, -0.0], [1.0, 0.0, 0.0] ) @@ -143,9 +166,24 @@ end T = get_transformation(panel) @test inv(T)(T(mp) + T(p)) ≈ p atol = 1e-6 + @test panel_area(panel) ≈ 1.0 atol = 1e-6 + @test normal_vector(panel) ≈ [0,0,-2.0] atol = 1e-6 +end + +@testset "Freestream 3D Velocity Conversion" begin + φ, θ = 1.0, 1.0 + fs = Freestream(alpha = φ, beta = θ) + V_test = [ cosd(θ) * cosd(φ), -sind(φ), sind(θ) * cosd(φ) ] + V_run = velocity(fs) + + φ_test, θ_test = cartesian_to_freestream(V_run) + + @test V_run ≈ V_test atol = 1e-6 + @test φ_test ≈ φ atol = 1e-6 + @test θ_test ≈ -θ atol = 1e-6 end -@testset "Vortex Lattice Method (Incompressible) - NACA 0012 Tapered Wing" begin +@testset "Vortex Lattice Method (Horseshoes, Incompressible) - NACA 0012 Tapered Wing" begin # Define wing wing = Wing( foils = [ naca4((0,0,1,2)) for i ∈ 1:2 ], @@ -171,7 +209,7 @@ end aircraft = ComponentArray(wing = make_horseshoes(WingMesh(wing, [20], 5, span_spacing = [Sine(1); Sine()]))) # Evaluate stability case - system = VortexLatticeSystem(aircraft, fs, refs) + system = VortexLatticeSystem(aircraft, fs, refs, false) dv_data = freestream_derivatives(system) dcf = dv_data.wing @@ -182,13 +220,13 @@ end # Test values nf_tests = [0.0013533, 0.0002199, 0.1159468, 0.0009056, 0.000387, -9.45e-5] ff_tests = [0.0014281, 0.0001743, 0.1159415] - dv_tests = [ 0.0797419 -0.0004936 -0.0039018 0.0637380 0.0004044 - 0.0066033 0.0024036 0.0809152 0.0302107 -0.0121575 - 3.4110222 -0.0060238 -0.185695 5.179196 0.0065573 - 0.0084357 0.0064223 0.285074 0.1109742 -0.0385723 - -0.0088809 0.0004686 0.0618588 -0.6604790 -0.0051181 - -0.0009645 -0.0018438 0.0054144 -0.0044697 0.0013694] - + dv_tests = [ 0.0797418 -0.0004935 -0.0039018 0.06373702665595818 0.0004043938437192077; + 0.0066032 0.0024036 0.0809152 0.030210659199330615 -0.012157465112956763; + 3.4110144 -0.0060237 -0.1856945 5.1791598961593985 0.006557296956889664; + 0.0084356 0.0064222 0.2850743 0.11097309781645154 -0.03857231877442366; + -0.0088809 0.0004686 0.0618588 -0.6604736234860606 -0.005118057094076678; + -0.0009645 -0.0018438 0.0054143 -0.004469695166159752 0.001369403002390149] + # Nearfield coefficients test [ @test nf_c ≈ nf_t atol = 1e-6 for (nf_c, nf_t) in zip(nfs, nf_tests) ] # Farfield coefficients test @@ -197,181 +235,64 @@ end [ @test dv_c ≈ dv_t atol = 1e-6 for (dv_c, dv_t) in zip(dvs, dv_tests) ] end -@testset "Vortex Lattice Method (Incompressible) - Vanilla Aircraft" begin - ## Wing - wing = Wing( - foils = fill(naca4((0,0,1,2)), 2), - chords = [1.0, 0.6], - twists = [0.0, 0.0], - spans = [5.0] / 2, - dihedrals = [11.39], - sweeps = [0.], - symmetry = true, - ); - - # Horizontal tail - htail = Wing( - foils = fill(naca4((0,0,1,2)), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.25] / 2, - dihedrals = [0.], - sweeps = [6.39], - position = [4., 0, 0], - angle = -2., - axis = [0., 1., 0.], - symmetry = true - ) - - # Vertical tail - vtail = Wing( - foils = fill(naca4((0,0,0,9)), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.0], - dihedrals = [0.], - sweeps = [7.97], - position = [4., 0, 0], - angle = 90., - axis = [1., 0., 0.], - ) - - ## Assembly - wing_mesh = WingMesh(wing, [32], 10; span_spacing = Cosine()) - htail_mesh = WingMesh(htail, [12], 6; span_spacing = Cosine()) - vtail_mesh = WingMesh(vtail, [5], 6; span_spacing = Cosine()) +# Include common aircraft definition +include(string(@__DIR__, "/aircraft_definition.jl")) +@testset "Vortex Lattice Method (Horseshoes, Compressible) - Vanilla Aircraft" begin aircraft = ComponentArray( wing = make_horseshoes(wing_mesh), htail = make_horseshoes(htail_mesh), vtail = make_horseshoes(vtail_mesh), ) - - ## Reference quantities - fs = Freestream( - alpha = 1.0, - beta = 1.0, - omega = zeros(3) - ) - - refs = References( - speed = 1.0, - area = projected_area(wing), - span = span(wing), - chord = mean_aerodynamic_chord(wing), - density = 1.225, - location = [0.25 * mean_aerodynamic_chord(wing), 0., 0.] - ) - ## Stability case - system = VortexLatticeSystem(aircraft, fs, refs) + system = VortexLatticeSystem(aircraft, fs, refs, true) dv_data = freestream_derivatives(system) dcf = dv_data.aircraft nfs = @views dcf[1:6,1] ffs = @views dcf[7:9,1] - dvs = @views dcf[1:6,3:end] + dvs = @views dcf[1:6,2:end] - nf_tests = [0.0003809, -0.0092002, 0.0653607, -0.0026762, 0.0601839, 0.0059991] - ff_tests = [0.0005143, -0.0092706, 0.0653322] - dv_tests = [ 0.0268317 0.0059901 0.0012994 0.0767768 -0.0023994 - 0.0048546 -0.517466 0.30401 -0.0744164 -0.899784 - 4.7952172 0.0598145 -0.0557912 11.9681911 0.0348196 - 0.0413913 -0.155083 0.489466 0.1469693 -0.0990389 - -1.5608632 -0.0360026 -0.0616955 -25.5170944 -0.124964 - 0.010348 0.336414 0.0155232 0.0906727 0.693076 ] + nf_tests = [0.0004394, -0.0096616, 0.0710165, -0.0027487, 0.0643144, 0.0063322] + ff_tests = [0.0005934, -0.0097369, 0.0709867] + dv_tests = [ 0.000303126 0.0318917 0.00657018 0.00105799 0.0840928 -0.000608665 + -0.00222327 0.00632911 -0.542542 0.31718 -0.0952224 -1.07032 + 0.0285508 5.164033 0.0671273 -0.0606896 14.2909359 0.0398133 + -0.000318489 0.0468793 -0.15929 0.514803 0.184925 -0.11513 + 0.0205454 -1.5507 -0.0366944 -0.0862627 -29.7754714 -0.149386 + 0.00161937 0.0117921 0.35447 0.0190155 0.119573 0.825448 ] # Nearfield coefficients test [ @test nf_c ≈ nf_t atol = 1e-6 for (nf_c, nf_t) in zip(nfs, nf_tests) ] # Farfield coefficients test [ @test ff_c ≈ ff_t atol = 1e-6 for (ff_c, ff_t) in zip(ffs, ff_tests) ] - # Stability derivatives' coefficients test + # # Stability derivatives' coefficients test [ @test dv_c ≈ dv_t atol = 1e-6 for (dv_c, dv_t) in zip(dvs, dv_tests) ] end -@testset "Vortex Lattice Method (Compressible) - Vanilla Aircraft" begin - ## Wing - wing = Wing( - foils = fill(naca4((0,0,1,2)), 2), - chords = [1.0, 0.6], - twists = [0.0, 0.0], - spans = [5.0] / 2, - dihedrals = [11.39], - sweeps = [0.], - symmetry = true, - ); - - # Horizontal tail - htail = Wing( - foils = fill(naca4((0,0,1,2)), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.25] / 2, - dihedrals = [0.], - sweeps = [6.39], - position = [4., 0, 0], - angle = -2., - axis = [0., 1., 0.], - symmetry = true - ) - - # Vertical tail - vtail = Wing( - foils = fill(naca4((0,0,0,9)), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.0], - dihedrals = [0.], - sweeps = [7.97], - position = [4., 0, 0], - angle = 90., - axis = [1., 0., 0.], - ) - - ## Assembly - wing_mesh = WingMesh(wing, [32], 10; span_spacing = Cosine()) - htail_mesh = WingMesh(htail, [12], 6; span_spacing = Cosine()) - vtail_mesh = WingMesh(vtail, [5], 6; span_spacing = Cosine()) - +@testset "Vortex Lattice Method (Rings, Compressible) - Vanilla Aircraft" begin aircraft = ComponentArray( - wing = make_horseshoes(wing_mesh), - htail = make_horseshoes(htail_mesh), - vtail = make_horseshoes(vtail_mesh), - ) - - ## Reference quantities - fs = Freestream( - alpha = 1.0, - beta = 1.0, - omega = zeros(3) - ) - - refs = References( - speed = 150.0, - area = projected_area(wing), - span = span(wing), - chord = mean_aerodynamic_chord(wing), - density = 1.225, - location = [0.25 * mean_aerodynamic_chord(wing), 0., 0.] + wing = make_vortex_rings(wing_mesh), + htail = make_vortex_rings(htail_mesh), + vtail = make_vortex_rings(vtail_mesh), ) ## Stability case - system = VortexLatticeSystem(aircraft, fs, refs) - dv_data = freestream_derivatives(system) + dv_data = freestream_derivatives(aircraft, fs, refs, name = :aircraft, compressible = true) dcf = dv_data.aircraft nfs = @views dcf[1:6,1] ffs = @views dcf[7:9,1] dvs = @views dcf[1:6,2:end] - nf_tests = [0.0004394, -0.0096616, 0.0710165, -0.0027487, 0.0643144, 0.0063322] - ff_tests = [0.0005934, -0.0097369, 0.0709867] - dv_tests = [ 0.000303126 0.0318917 0.00657018 0.00105799 0.0840928 -0.000608665 - -0.00222327 0.00632911 -0.542542 0.31718 -0.0952224 -1.07032 - 0.0285508 5.164033 0.0671273 -0.0606896 14.2909359 0.0398133 - -0.000318489 0.0468793 -0.15929 0.514803 0.184925 -0.11513 - 0.0205454 -1.5507 -0.0366944 -0.0862627 -29.7754714 -0.149386 - 0.00161937 0.0117921 0.35447 0.0190155 0.119573 0.825448 ] + nf_tests = [0.0004485, -0.0102081, 0.0706159, -0.0028898, 0.0653328, 0.0067422] + ff_tests = [0.0006049, -0.0102852, 0.0705857] + dv_tests = [ 0.000306494 0.0316461 0.00767591 0.00118087 0.0837572 0.000839109 + -0.00232925 0.00389612 -0.565213 0.315484 -0.0809676 -1.13528256 + 0.0283188 5.1464928 0.0345146 -0.04648 14.2846332 0.0138955 + -0.000401224 0.0417969 -0.167955 0.515156 0.16723 -0.120098 + 0.021139 -1.502605 0.0895636 -0.0791555 -29.800452 -0.0345743 + 0.00166363 0.0116769 0.370111 0.0211191 0.102458 0.882128 ] # Nearfield coefficients test [ @test nf_c ≈ nf_t atol = 1e-6 for (nf_c, nf_t) in zip(nfs, nf_tests) ] From ab737d644376a63d1bdec66a7bdbb623b58348fa Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 18:57:16 +0800 Subject: [PATCH 07/13] Updated VLM examples with vortex rings --- .../3d_vortex_lattice_method/vlm_aircraft.jl | 343 +++++++----------- .../vlm_aircraft_makie_plot.jl | 115 ++++++ .../3d_vortex_lattice_method/vlm_implicit.jl | 120 ------ .../3d_vortex_lattice_method/vlm_sand.jl | 123 ------- .../vlm_vortex_rings.jl | 10 +- .../3d_vortex_lattice_method/vlm_wing.jl | 7 +- 6 files changed, 259 insertions(+), 459 deletions(-) create mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft_makie_plot.jl delete mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl delete mode 100644 examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft.jl index ab4e510f..2195e1f8 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft.jl @@ -51,20 +51,20 @@ print_info(htail, "Horizontal Tail") print_info(vtail, "Vertical Tail") ## WingMesh type -wing_mesh = WingMesh(wing, [24], 6, +wing_mesh = WingMesh(wing, [24], 6, # span_spacing = Cosine() ) htail_mesh = WingMesh(htail, [24], 6, # span_spacing = Cosine() ) vtail_mesh = WingMesh(vtail, [24], 6, - span_spacing = Cosine() + # span_spacing = Cosine() ) aircraft = ComponentVector( - wing = make_horseshoes(wing_mesh), - htail = make_horseshoes(htail_mesh), - vtail = make_horseshoes(vtail_mesh) + wing = make_vortex_rings(wing_mesh), + htail = make_vortex_rings(htail_mesh), + vtail = make_vortex_rings(vtail_mesh) ); ## Case @@ -84,28 +84,30 @@ ref = References( location = mean_aerodynamic_center(wing) ) -## -@time system = solve_case( +## Solve +@time sys = solve_case( aircraft, fs, ref; + compressible = true, # Compressibility correction flag print = true, # Prints the results for only the aircraft print_components = true, # Prints the results for all components ); -## Compute dynamics +## Compute forces, moments and velocities over each surface ax_sys = Wind() # Axis systems: Geometry(), Stability(), Body() -@time CFs, CMs = surface_coefficients(system; axes = ax_sys) -# Fs, Ms = surface_dynamics(system; axes = ax) -# Fs = surface_forces(system; axes = ax) -# vels = surface_velocities(system) +@time CFs, CMs = surface_coefficients(sys; axes = ax_sys) # Coefficients +# Fs, Ms = surface_dynamics(sys; axes = ax) # Forces and moments +# Fs = surface_forces(sys; axes = ax) # Forces only +# vels = surface_velocities(sys) # Velocities -nf = nearfield(system) -ff = farfield(system) +## Aerodynamic coefficients +nf = nearfield(sys) +ff = farfield(sys) -nfs = nearfield_coefficients(system) -ffs = farfield_coefficients(system) +nfs = nearfield_coefficients(sys) +ffs = farfield_coefficients(sys) ## Force/moment coefficients and derivatives -@time dvs = freestream_derivatives(system; +@time dvs = freestream_derivatives(sys; axes = ax_sys, print = true, print_components = true, @@ -115,21 +117,21 @@ ffs = farfield_coefficients(system) ## Viscous drag prediction # Equivalent flat-plate skin friction estimation -CDv_wing = profile_drag_coefficient(wing, [0.8, 0.8], system.reference) -CDv_htail = profile_drag_coefficient(htail, [0.6, 0.6], system.reference) -CDv_vtail = profile_drag_coefficient(vtail, [0.6, 0.6], system.reference) +CDv_wing = profile_drag_coefficient(wing, [0.8, 0.8], sys.reference) +CDv_htail = profile_drag_coefficient(htail, [0.6, 0.6], sys.reference) +CDv_vtail = profile_drag_coefficient(vtail, [0.6, 0.6], sys.reference) CDv_plate = CDv_wing + CDv_htail + CDv_vtail ## Local dissipation form factor friction estimation import LinearAlgebra: norm -edge_speeds = norm.(surface_velocities(system)); # Inviscid speeds on the surfaces +edge_speeds = norm.(surface_velocities(sys)); # Inviscid speeds on the surfaces # Drag coefficients -CDvd_wing = profile_drag_coefficient(wing_mesh, [0.8, 0.8], edge_speeds.wing, system.reference) -CDvd_htail = profile_drag_coefficient(htail_mesh, [0.6, 0.6], edge_speeds.htail, system.reference) -CDvd_vtail = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], edge_speeds.vtail, system.reference) +CDvd_wing = profile_drag_coefficient(wing_mesh, [0.8, 0.8], edge_speeds.wing, sys.reference) +CDvd_htail = profile_drag_coefficient(htail_mesh, [0.6, 0.6], edge_speeds.htail, sys.reference) +CDvd_vtail = profile_drag_coefficient(vtail_mesh, [0.6, 0.6], edge_speeds.vtail, sys.reference) CDv_diss = CDvd_wing + CDvd_htail + CDvd_vtail @@ -137,166 +139,84 @@ CDv_diss = CDvd_wing + CDvd_htail + CDvd_vtail CDv = CDv_diss ## Total force coefficients with empirical viscous drag prediction -CDi_nf, CY_nf, CL_nf, Cl, Cm, Cn = nf = nearfield(system) -CDi_ff, CY_ff, CL_ff = ff = farfield(system) +CDi_nf, CY_nf, CL_nf, Cl, Cm, Cn = nf = nearfield(sys) +CDi_ff, CY_ff, CL_ff = ff = farfield(sys) nf_v = [ CDi_nf + CDv; CDv; nf ] ff_v = [ CDi_ff + CDv; CDv; ff ] print_coefficients(nf_v, ff_v) -## Spanwise forces/lifting line loads -wing_ll = spanwise_loading(wing_mesh, CFs.wing, S) -htail_ll = spanwise_loading(htail_mesh, CFs.htail, S) -vtail_ll = spanwise_loading(vtail_mesh, CFs.vtail, S); - ## Plotting -using CairoMakie -CairoMakie.activate!() +#=========================================================# -set_theme!( - # theme_black() - # theme_light() - ) +using Plots +gr() using LaTeXStrings const LS = LaTeXString -## Streamlines -# Spanwise distribution -span_points = 20 -init = chop_leading_edge(wing, span_points) -dx, dy, dz = 0, 0, 1e-3 -seed = init .+ Ref([dx, dy, dz]) - # init .+ Ref([dx, dy, -dz]) ]; - -distance = 5 -num_stream_points = 100 -streams = streamlines(system, seed, distance, num_stream_points); - -wing_cam_connec = triangle_connectivities(LinearIndices(wing_mesh.camber_mesh)) -htail_cam_connec = triangle_connectivities(LinearIndices(htail_mesh.camber_mesh)) -vtail_cam_connec = triangle_connectivities(LinearIndices(vtail_mesh.camber_mesh)); - -## Surface velocities -vels = surface_velocities(system); -sps = norm.(vels) - -wing_sp_points = extrapolate_point_mesh(sps.wing) -htail_sp_points = extrapolate_point_mesh(sps.htail) -vtail_sp_points = extrapolate_point_mesh(sps.vtail) - -## Surface pressure coefficients -cps = norm.(CFs) * S - -wing_cp_points = extrapolate_point_mesh(cps.wing) -htail_cp_points = extrapolate_point_mesh(cps.htail) -vtail_cp_points = extrapolate_point_mesh(cps.vtail) - -## Figure plot -fig1 = Figure(resolution = (1280, 720)) - -scene = LScene(fig1[1:4,1]) -ax = fig1[1:4,2] = GridLayout() - -ax_cd = Axis(ax[1,1:2], ylabel = L"C_{D_i}", title = LS("Spanwise Loading")) -ax_cy = Axis(ax[2,1:2], ylabel = L"C_Y",) -ax_cl = Axis(ax[3,1:2], xlabel = L"y", ylabel = L"C_L") - -# Spanload plot -function plot_spanload!(ax, ll_loads, name = "Wing") - @views lines!(ax[1,1:2], ll_loads[:,1], ll_loads[:,2], label = name,) - @views lines!(ax[2,1:2], ll_loads[:,1], ll_loads[:,3], label = name,) - @views lines!(ax[3,1:2], ll_loads[:,1], ll_loads[:,4], label = name,) - - nothing +## Coordinates +Plots.plot( + aspect_ratio = 1, + camera = (30, 30), + zlim = span(wing) .* (-0.5, 0.5), + size = (800, 600) +) +Plots.plot!(wing_mesh, label = "Wing") +Plots.plot!(htail_mesh, label = "Wing") +Plots.plot!(vtail_mesh, label = "Wing") +Plots.plot!(sys, wing_mesh, dist = 10, num_stream = 50, span = 10, color = :green) +Plots.plot!(sys, htail_mesh, dist = 2, num_stream = 50, span = 10, color = :green) +Plots.plot!(sys, vtail_mesh, dist = 2, num_stream = 50, span = 10, color = :green) + +## Compute spanwise loads +wing_ll = spanwise_loading(wing_mesh, ref, CFs.wing, sys.circulations.wing) +htail_ll = spanwise_loading(htail_mesh, ref, CFs.htail, sys.circulations.htail) +vtail_ll = spanwise_loading(vtail_mesh, ref, CFs.vtail, sys.circulations.vtail); + +## Plot spanwise loadings +plot_CD = begin + Plots.plot(label = :none, ylabel = "CDi") + Plots.plot!(wing_ll[:,1], wing_ll[:,2], label = "Wing") + Plots.plot!(htail_ll[:,1], htail_ll[:,2], label = "Horizontal Tail") + Plots.plot!(vtail_ll[:,1], vtail_ll[:,2], label = "Vertical Tail") end -plot_spanload!(ax, wing_ll, LS("Wing")) -plot_spanload!(ax, htail_ll, LS("Horizontal Tail")) -plot_spanload!(ax, vtail_ll, LS("Vertical Tail")) - -# Legend -Legend(ax[4,1:2], ax_cl) - -# Surface pressure meshes -m1 = poly!(scene, vec(wing_mesh.camber_mesh), wing_cam_connec, color = vec(wing_cp_points)) -m2 = poly!(scene, vec(htail_mesh.camber_mesh), htail_cam_connec, color = vec(htail_cp_points)) -m3 = poly!(scene, vec(vtail_mesh.camber_mesh), vtail_cam_connec, color = vec(vtail_cp_points)) - -Colorbar(fig1[4,1], m1.plots[1], label = L"Pressure Coefficient, $C_p$", vertical = false) - -fig1[0, :] = Label(fig1, LS("Vortex Lattice Analysis"), textsize = 20) - - -# Airfoil meshes -# wing_surf = surface_coordinates(wing_mesh, wing_mesh.n_span, 60) -# surf_connec = triangle_connectivities(LinearIndices(wing_surf)) -# wing_surf_mesh = mesh(vec(wing_surf), surf_connec) -# w1 = wireframe!(scene, wing_surf_mesh.plot[1][], color = :grey, alpha = 0.1) - -# Borders -lines!(scene, plot_planform(wing)) -lines!(scene, plot_planform(htail)) -lines!(scene, plot_planform(vtail)) - -# l1 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(wing_mesh)) ] -# l2 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(htail_mesh)) ] -# l3 = [ lines!(scene, pts, color = :grey) for pts in plot_panels(camber_panels(vtail_mesh)) ] - -# Streamlines -[ lines!(scene, Point3f.(stream[:]), color = :green) for stream in eachcol(streams) ] - -fig1.scene - -## Save figure -# save("plots/VortexLattice.pdf", fig1, px_per_unit = 1.5) - -## Animation settings -pts = [ Point3f[stream] for stream in streams[1,:] ] - -[ lines!(scene, pts[i], color = :green, axis = (; type = Axis3)) for i in eachindex(pts) ] - -# Recording -fps = 30 -nframes = length(streams[:,1]) +plot_CY = begin + Plots.plot(label = :none, ylabel = "CY") + Plots.plot!(wing_ll[:,1], wing_ll[:,3], label = "Wing") + Plots.plot!(htail_ll[:,1], htail_ll[:,3], label = "Horizontal Tail") + Plots.plot!(vtail_ll[:,1], vtail_ll[:,3], label = "Vertical Tail") +end -# record(fig, "plots/vlm_animation.mp4", 1:nframes) do i -# for j in eachindex(streams[1,:]) -# pts[j][] = push!(pts[j][], Point3f0(streams[i,j])) -# end -# sleep(1/fps) # refreshes the display! -# notify(pts[i]) -# end +plot_CL = begin + Plots.plot(label = :none, ylabel = "CL") + Plots.plot!(wing_ll[:,1], wing_ll[:,4], label = "Wing") + Plots.plot!(htail_ll[:,1], htail_ll[:,4], label = "Horizontal Tail") + Plots.plot!(vtail_ll[:,1], vtail_ll[:,4], label = "Vertical Tail") + Plots.plot!(wing_ll[:,1], wing_ll[:,5], label = "Wing Normalized", xlabel = "y") +end -## Arrows -# hs_pts = vec(Tuple.(bound_leg_center.(horses))) -# arrows!(scene, getindex.(hs_pts, 1), getindex.(hs_pts, 2), getindex.(hs_pts, 3), -# vec(CDis),vec(CYs),vec( Ls), -# arrowsize = Vec3f.(0.3, 0.3, 0.4), -# lengthscale = 10, -# label = "Forces (Exaggerated)") +Plots.plot(plot_CD, plot_CY, plot_CL, size = (800, 700), layout = (3,1)) ## VARIABLE ANALYSES #=========================================================# using Accessors using Base.Iterators: product -using Plots -pyplot() - -## Speed sweep +## Speed variation Vs = 1.0:10:300 res_Vs = combinedimsview( map(Vs) do V ref1 = @set ref.speed = V - sys = solve_case(aircraft, fs, ref1) - [ mach_number(ref1); farfield(sys)...; nearfield(sys)... ] + sys = solve_case(aircraft, fs, ref1, compressible = true) + [ mach_number(ref1); farfield(sys); nearfield(sys) ] end, (1) ) -## +## Plot Plots.plot( res_Vs[:,1], res_Vs[:,2:end], layout = (3,3), size = (900, 800), @@ -304,7 +224,65 @@ Plots.plot( labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"] ) -## +## Alpha variation +αs = -5:0.5:5 +res_αs = combinedimsview( + map(αs) do α + fst = @set fs.alpha = deg2rad(α) + sys = solve_case(aircraft, fst, ref, compressible = true) + [ α; farfield(sys); nearfield(sys) ] + end, (1) +) + +## Plot +Plots.plot( + res_αs[:,1], res_αs[:,2:end], + layout = (3,3), size = (900, 800), + xlabel = "α", + labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"] +) + + +## (Speed, alpha) variation +res = combinedimsview( + map(product(Vs, αs)) do (V, α) + ref1 = @set ref.speed = V + fst = @set fs.alpha = deg2rad(α) + sys = solve_case(aircraft, fst, ref1, compressible = true) + [ mach_number(ref1); α; farfield(sys); nearfield(sys) ] + end +) + +res_p = permutedims(res, (3,1,2)) + +## CDi +plt_CDi_ff = Plots.plot(camera = (75, 30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{D_i}") +[ Plots.plot!(res_p[:,1,n], res_p[:,2,n], res_p[:,3,n], label = "", c = :black) for n in axes(res_p,3) ] + +# CL +plt_CL_ff = Plots.plot(camera = (75,15), +ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{L}", +label = "", c = :black) +[ Plots.plot!( + res_p[:,1,n], res_p[:,2,n], res_p[:,8,n], label = "", c = :black) for n in axes(res_p,3) ] + +# Cm +plt_Cm_ff = Plots.plot(camera = (100,30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_m", label = "") +[ Plots.plot!( + res_p[:,1,n], res_p[:,2,n], res_p[:,10,n], label = "", c = :black) for n in axes(res_p,3) ] + +# Plot +p = Plots.plot(plt_CDi_ff, plt_CL_ff, plt_Cm_ff, layout = (1,3), size = (1300, 400)) + +## Save +savefig(p, "coeffs.png") + +## Fancy plotting using Makie, if needed +#==================================================================# + +include("vlm_aircraft_makie_plot.jl") + +## Speed variation f2 = Figure() ax1 = Axis(f2[1,1], xlabel = L"M", ylabel = L"C_{D_{ff}}") lines!(res_Vs[:,1], res_Vs[:,2]) # CD_ff @@ -327,24 +305,7 @@ lines!(res_Vs[:,1], res_Vs[:,10]) f2 -## Alpha sweep -αs = -5:0.5:5 -res_αs = combinedimsview( - map(αs) do α - fst = @set fs.alpha = deg2rad(α) - sys = solve_case(aircraft, fst, ref) - [ α; farfield(sys)...; nearfield(sys)... ] - end, (1) -) - -Plots.plot( - res_αs[:,1], res_αs[:,2:end], - layout = (3,3), size = (900, 800), - xlabel = "α", - labels = ["CD_ff" "CY_ff" "CL_ff" "CD" "CY" "CL" "Cl" "Cm" "Cn"] -) - -## Makie +## Alpha variation f3 = Figure() ax1 = Axis(f3[1,1], xlabel = L"α", ylabel = L"C_{D_{ff}}") lines!(res_αs[:,1], res_αs[:,2]) # CD_ff @@ -365,38 +326,4 @@ lines!(res_αs[:,1], res_αs[:,9]) ax9 = Axis(f3[3,3], xlabel = L"M", ylabel = L"C_n") lines!(res_αs[:,1], res_αs[:,10]) -f3 - -## (Speed, alpha) sweep -res = combinedimsview( - map(product(Vs, αs)) do (V, α) - ref1 = @set ref.speed = V - fst = @set fs.alpha = deg2rad(α) - sys = solve_case(aircraft, fst, ref1) - [ mach_number(ref1); α; farfield(sys)...; nearfield(sys)... ] - end -) - -## -res_p = permutedims(res, (3,1,2)) - -# CDi -plt_CDi_ff = Plots.plot(camera = (75, 30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{D_i}") -[ Plots.plot!(res_p[:,1,n], res_p[:,2,n], res_p[:,3,n], label = "", c = :black) for n in axes(res_p,3) ] - -# CL -plt_CL_ff = Plots.plot(camera = (75,15), -ylabel = L"\alpha", xlabel = "V", zlabel = L"C_{L}", -label = "", c = :black) -[ Plots.plot!( - res_p[:,1,n], res_p[:,2,n], res_p[:,8,n], label = "", c = :black) for n in axes(res_p,3) ] - -plt_Cm_ff = Plots.plot(camera = (100,30), ylabel = L"\alpha", xlabel = "V", zlabel = L"C_m", label = "") -[ Plots.plot!( - res_p[:,1,n], res_p[:,2,n], res_p[:,10,n], label = "", c = :black) for n in axes(res_p,3) ] - -## -p = Plots.plot(plt_CDi_ff, plt_CL_ff, plt_Cm_ff, layout = (1,3), size = (1300, 400)) - -## -savefig(p, "coeffs.png") \ No newline at end of file +f3 \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft_makie_plot.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft_makie_plot.jl new file mode 100644 index 00000000..7cac038d --- /dev/null +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_aircraft_makie_plot.jl @@ -0,0 +1,115 @@ +using CairoMakie +CairoMakie.activate!() + +# using WGLMakie +# WGLMakie.activate!() + +set_theme!( + # theme_black() + # theme_light() +) + +## Streamlines +# Spanwise distribution +span_points = 20 +dx, dy, dz = 0, 0, 1e-3 +init = chop_leading_edge(wing, span_points) +seed = init .+ Ref([dx, dy, dz]) +# init .+ Ref([dx, dy, -dz]) ]; + +distance = 5 +num_stream_points = 100 +streams = streamlines(sys, seed, distance, num_stream_points); + +wing_cam_connec = triangle_connectivities(LinearIndices(camber_coordinates(wing_mesh))) +htail_cam_connec = triangle_connectivities(LinearIndices(camber_coordinates(htail_mesh))) +vtail_cam_connec = triangle_connectivities(LinearIndices(camber_coordinates(vtail_mesh))); + +## Surface velocities +vels = surface_velocities(sys); +sps = norm.(vels) + +wing_sp_points = extrapolate_point_mesh(sps.wing) +htail_sp_points = extrapolate_point_mesh(sps.htail) +vtail_sp_points = extrapolate_point_mesh(sps.vtail) + +## Surface pressure coefficients +cps = norm.(CFs) * S + +wing_cp_points = extrapolate_point_mesh(cps.wing) +htail_cp_points = extrapolate_point_mesh(cps.htail) +vtail_cp_points = extrapolate_point_mesh(cps.vtail) + +## Figure plot +fig1 = Figure(resolution = (1280, 720)) + +scene = LScene(fig1[1:4,1]) +ax = fig1[1:4,2] = GridLayout() + +ax_cd = Axis(ax[1,1:2], ylabel = L"C_{D_i}", title = LS("Spanwise Loading")) +ax_cy = Axis(ax[2,1:2], ylabel = L"C_Y",) +ax_cl = Axis(ax[3,1:2], xlabel = L"y", ylabel = L"C_L") + +# Spanload plot +function plot_spanload!(ax, ll_loads, name = "Wing") + @views lines!(ax[1,1:2], ll_loads[:,1], ll_loads[:,2], label = name,) + @views lines!(ax[2,1:2], ll_loads[:,1], ll_loads[:,3], label = name,) + @views lines!(ax[3,1:2], ll_loads[:,1], ll_loads[:,4], label = name,) + + nothing +end + +# Surface pressure meshes +m1 = poly!(scene, vec(camber_coordinates(wing_mesh)), wing_cam_connec, color = vec(wing_cp_points), colormap = :rainbow1) +m2 = poly!(scene, vec(camber_coordinates(htail_mesh)), htail_cam_connec, color = vec(htail_cp_points), colormap = :rainbow1) +m3 = poly!(scene, vec(camber_coordinates(vtail_mesh)), vtail_cam_connec, color = vec(vtail_cp_points), colormap = :rainbow1) + +Colorbar(fig1[4,1], m1.plots[1], label = L"Pressure Coefficient, $C_p$", vertical = false) + +fig1[0, :] = Label(fig1, LS("Vortex Lattice Analysis"), fontsize = 20) + + +# Airfoil meshes +# wing_surf = surface_coordinates(wing_mesh, wing_mesh.n_span, 60) +# surf_connec = triangle_connectivities(LinearIndices(wing_surf)) +# wing_surf_mesh = mesh(vec(wing_surf), surf_connec) +# w1 = wireframe!(scene, wing_surf_mesh.plot[1][], color = :grey, alpha = 0.1) + +# Borders +lines!(scene, plot_planform(wing)) +lines!(scene, plot_planform(htail)) +lines!(scene, plot_planform(vtail)) + +# Streamlines +[ lines!(scene, Point3f.(stream[:]), color = :lightblue) for stream in eachcol(streams) ] + +plot_spanload!(ax, wing_ll, LS("Wing")) +plot_spanload!(ax, htail_ll, LS("Horizontal Tail")) +plot_spanload!(ax, vtail_ll, LS("Vertical Tail")) + +# Legend +Legend(ax[4,1:2], ax_cl) + +fig1.scene + +## Save figure +# save("plots/VortexLattice.pdf", fig1, px_per_unit = 1.5) + +## Animation +#=======================================================# + +# pts = [ Point3f[stream] for stream in streams[1,:] ] + +# [ lines!(scene, pts[i], color = :green) for i in eachindex(pts) ] + +# Recording (Doesn't work right now) +# fps = 30 +# nframes = length(streams[:,1]) + +# record(fig1, "plots/vlm_animation.mp4", 1:nframes) do i +# for j in eachindex(streams[1,:]) +# pts[j][] = push!(pts[j][], Point3f(streams[i,j])) +# end +# sleep(1/fps) # refreshes the display! +# notify(pts[i]) +# end \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl deleted file mode 100644 index 63a779fa..00000000 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_implicit.jl +++ /dev/null @@ -1,120 +0,0 @@ -## VLM implicit analysis and derivatives -using AeroFuse -using LinearAlgebra -using ImplicitAD -using NLsolve -using StructArrays - -## Surfaces -#================================================# - -function make_surfaces(x) - ## Wing - wing = Wing( - foils = fill(naca4(2,4,1,2), 2), - chords = [1.0, 0.6], - twists = [2.0, 0.0], - spans = [4.0], - dihedrals = [5.], - sweeps = [5.], - symmetry = true - ) - - ## Horizontal tail - htail = Wing( - foils = fill(naca4(0,0,1,2), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.25], - dihedrals = [0.], - sweeps = [6.39], - position = [4., 0, -0.1], - angle = -2., - axis = [0., 1., 0.], - symmetry = true - ) - - ## Vertical tail - vtail = Wing( - foils = fill(naca4(0,0,0,9), 2), - chords = [0.7, 0.42], - twists = [0.0, 0.0], - spans = [1.0], - dihedrals = [0.], - sweeps = [7.97], - position = [4., 0, 0], - angle = 90., - axis = [1., 0., 0.] - ) - - ## Meshing - #================================================# - - wing_mesh = WingMesh(wing, [24], 12, span_spacing = Cosine()) - htail_mesh = WingMesh(htail, [12], 6, span_spacing = Cosine()) - vtail_mesh = WingMesh(vtail, [12], 6, span_spacing = Cosine()) - - # Assemble - aircraft = ComponentArray( - wing = make_horseshoes(wing_mesh), - htail = make_horseshoes(htail_mesh), - vtail = make_horseshoes(vtail_mesh) - ) - - return aircraft -end - -ac = make_surfaces([1.0, 0.6]) - -## VLM analysis -#================================================# - -fs = Freestream( - alpha = 0.0, - beta = 0.0, - omega = [0., 0., 0.] -) - -ref = References( - speed = 150.0, - density = 1.225, - viscosity = 1.5e-5, -) - -function solve(system) - rwrap(R, p) = AeroFuse.VortexLattice.residual!(R, p, system) - res = nlsolve( - rwrap, - zero(system.circulations), - method = :newton, - autodiff = :forward, - show_trace = true - ) - return res.zero -end - -## -function program(x) - ac = make_surfaces(x) - fs = Freestream( - alpha = convert(eltype(x), 5.0), - beta = 0.0, - omega = [0., 0., 0.] - ) - - ref = References( - speed = 150.0, - density = 1.225, - viscosity = 1.5e-5, - ) - - sys = VortexLatticeSystem(ac, fs, ref) - - solve(sys) - # [ velocity(fs) sum(ac.rc) ] -end - -## -using ForwardDiff - -ForwardDiff.jacobian(program, [1.1, 0.6]) diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl deleted file mode 100644 index 1e2662cb..00000000 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_sand.jl +++ /dev/null @@ -1,123 +0,0 @@ -## -using AeroFuse -using JuMP -using Ipopt -using NLsolve - -## Wing section setup -wing = Wing( - foils = fill(naca4((2,4,1,2)), 3), - chords = [2.0, 1.6, 0.2], - twists = [0.0, 0.0, 0.0], - spans = [5., 0.6], - dihedrals = [5., 5.], - sweeps = [20.,20.], - w_sweep = 0.25, # Quarter-chord sweep - symmetry = true, - # flip = true -) - -## -x_w, y_w, z_w = wing_mac = mean_aerodynamic_center(wing) -print_info(wing, "Wing") - -## Meshing and assembly -wing_mesh = WingMesh(wing, [24,12], 6, - span_spacing = Cosine() - ); -aircraft = ComponentVector(wing = make_horseshoes(wing_mesh)) - -# Freestream conditions -fs = Freestream( - alpha = 2.0, # deg - beta = 2.0, # deg - omega = [0.,0.,0.] -) - -# Reference values -ref = References( - speed = 150, # m/s - density = 1.225, # kg/m³ - viscosity = 1.5e-5, # ??? - area = projected_area(wing), # m² - span = span(wing), # m - chord = mean_aerodynamic_chord(wing), # m - location = mean_aerodynamic_center(wing) # m -) - -## Solve system -system = VortexLatticeSystem(aircraft, fs, ref) - -## NLsolve -R = similar(system.circulations) - -## -AeroFuse.solve_nonlinear!(R, x) = solve_nonlinear!(R, system.vortices, x ./ ref.speed, velocity(fs) / ref.speed, fs.omega) - -## -xz = nlsolve( - solve_nonlinear!, - system.circulations, - autodiff = :forward, - show_trace = true, -) - -n = length(system.circulations) - -## THIS IS THE IMPLICIT FUNCTION THEOREM APPLIED, R(x) = R(x, u(x)) -obtain_root(x) = nlsolve(u -> solve_nonlinear!(R, x, u), system.circulations).zero # Rootfinding method (doesn't matter which) - - -## ChainRules adjoint -## Evaluate -x0 = [1.2, -2.0] -obtain_root(x0) -loss_functional(obtain_root(x0)) - -## Define adjoint rule -using ChainRulesCore - -function ChainRulesCore.rrule(::typeof(obtain_root), x) - # Solve residual equation - # R(x,u) =(by implicit function theorem)= u -> R(u(x)) = 0 - u_class = obtain_root(x) - - # Define pullback - function obtain_root_pullback(∂J_∂u) - # Evaluate gradients ∂R/∂x, ∂R/∂u - ∂R_∂x, ∂R_∂u = ForwardDiff.jacobian!(residual, x, u_class) - - # Adjoint variable: - # λᵀ = ∂J/∂R = (∂R/∂u)^(-1) * ∂J/∂u - λᵀ = ∂R_∂u' \ ∂J_∂u - - # Loss sensitivity to input: - # ∂J/∂x = -λ^T * ∂R/∂x - dJ_dx = -λᵀ' * ∂R_∂x - - # Loss sensitivity to residual equation: - # ∂J/∂R = 0 - ∂J_∂R = NoTangent() - - return ∂J_∂R, dJ_dx' - end - - return u_class, obtain_root_pullback -end - - - - -## Create JuMP model, using Ipopt as the solver -model = Model(optimizer_with_attributes(Ipopt.Optimizer)) - -@variables(model, begin - Γs[j = 1:n] == system.circulations[j] -end) - -register(model, :R_A, n, solve_nonlinear; autodiff = true) - -# @NLconstraint(model, [j = 1:n], solve_nonlinear(Γs)[i] == 0) - -## -optimize!(model) \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl index 0c7e10ea..963ad979 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_vortex_rings.jl @@ -5,7 +5,7 @@ using StructArrays ## Wing section setup wing = Wing( - foils = fill(naca4((2,4,1,2)), 2), + foils = fill(naca4((0,0,1,2)), 2), chords = [2.2, 1.8], twists = [0.0, 0.0], spans = [7.5], @@ -20,13 +20,13 @@ wing_mesh = WingMesh(wing, [24], 8, span_spacing = Uniform()); # Freestream conditions fs = Freestream( alpha = 2.0, # deg - beta = 0.0, # deg + beta = 2.0, # deg omega = [0.,0.,0.] ) # Reference values ref = References( - speed = 15., # m/s + speed = 1., # m/s density = 1.225, # kg/m³ viscosity = 1.5e-5, # ??? area = projected_area(wing), # m² @@ -37,10 +37,10 @@ ref = References( ## Horseshoes ac_hs = ComponentVector(wing = make_horseshoes(wing_mesh)) -system = VortexLatticeSystem(ac_hs, fs, ref) +system = VortexLatticeSystem(ac_hs, fs, ref, false) print_coefficients(nearfield(system ), farfield(system)) ## Vortex rings ac_vs = ComponentVector(wing = make_vortex_rings(wing_mesh)) -sys = VortexLatticeSystem(ac_vs, fs, ref) +sys = VortexLatticeSystem(ac_vs, fs, ref, false) print_coefficients(nearfield(sys), farfield(sys)) \ No newline at end of file diff --git a/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl b/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl index dea01b31..765f44b4 100644 --- a/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl +++ b/examples/aerodynamics/3d_vortex_lattice_method/vlm_wing.jl @@ -47,6 +47,7 @@ ref = References( @time begin system = solve_case( aircraft, fs, ref; + compressible = true, # print = true, # Prints the results for only the aircraft # print_components = true, # Prints the results for all components ); @@ -100,14 +101,14 @@ using Plots gr() ## Coordinates -plot( +Plots.plot( aspect_ratio = 1, camera = (30, 30), zlim = span(wing) .* (-0.5, 0.5), size = (800, 600) ) -plot!(wing_mesh, label = "Wing") -plot!(system, wing, dist = 3, num_stream = 50, span = 10, color = :green) +Plots.plot!(wing_mesh, label = "Wing") +Plots.plot!(system, wing, dist = 3, num_stream = 50, span = 10, color = :green) ## Compute spanwise loads span_loads = spanwise_loading(wing_mesh, ref, CFs.wing, system.circulations.wing) From bdb324d028d68ea4ba403dff9d7b994dd750b58f Mon Sep 17 00:00:00 2001 From: Arjit Seth Date: Mon, 20 Mar 2023 18:57:33 +0800 Subject: [PATCH 08/13] Updated plots --- docs/src/assets/VortexLattice.svg | 6 +++--- plots/VortexLattice.pdf | Bin 173639 -> 531040 bytes plots/VortexLattice.svg | 6 +++--- plots/VortexLatticeLegacy.svg | 9 +++++++++ 4 files changed, 15 insertions(+), 6 deletions(-) create mode 100644 plots/VortexLatticeLegacy.svg diff --git a/docs/src/assets/VortexLattice.svg b/docs/src/assets/VortexLattice.svg index 70241b4c..5fcf153c 100644 --- a/docs/src/assets/VortexLattice.svg +++ b/docs/src/assets/VortexLattice.svg @@ -1,9 +1,9 @@ - + - - + + diff --git a/plots/VortexLattice.pdf b/plots/VortexLattice.pdf index 0cdcc3be7a6b454b32019c90369fefee0bba5b4b..1f926fac7a9a0f746e788f1544167d888b6a1fd5 100644 GIT binary patch literal 531040 zcmeF1Wm6r&+vah1cZc8_+}+(Jz`-4YOM<(*yCt}Da1VCi0Kp-+ySptP;Hm$Mt=+9> zx3=o4>H77|)!lc`^mI*6)2d0yuyC+)Bhj{;)xIEorC_IUHnT$#5~5&JvUIX`x1r$v z%+f%jprBxrv2}2_bp1>n0PdDjmgdeDmPjHZNN(<~mHygaSF4Z?y#M#$PRQCJLk{PrAj(C;=EMyCt?wJ zz1-XqUX#E9?}<|%x9KA9=Px%?s_*q5ub_|T$EuGck@v-3k^3m&$Dn|xmxKUsUo!H8 z3mwqNaos!FM3rlE_x`V~?BN2mFR$kt?{`meu^aUPHv!+DRRPs{f~S@65DPoLD9w2#omC&jf4E2pRXASAFr>saX$WQqK+T0cVBcb>i(|2 zHGjO_@m_AW^jcLM+;rBwEY|fTt}p=`Ss6BSHTGA?YZe4?&9Er>HS&o9%1iCbwaWAWv~BT*k2Re-axj)IuqrA z^SZ8AmM4>deeU@@{~yY&9RXL_=Lc6Kiv@Q;SJ43Yo_92!MaF{im5=9nd>P@XHOH5<>jvS2;=8W21y^IyZ!K>w^GOl<5sLybj*~a&odvpW zr$rKP&pfU!4KydN0H(=(}k7;N#cjVBq85d{@bk#~TG|<(M-6cz*k3F6Y}2vn=X( znq+?Hs!w8Q-kQK|Y}p##PEdEVt?}Xm<=?g1pj+D|b*`_QcYw$7bz#MGUig`(#J&2Y zn&s~O#3S#<0a?}fM_jJ05%!JVu+I%#hBabK&lFu9H2{6T;9(P8Q1GFt^b&35m}v{% zbERdc!x0EreS0MCZ1ypB6%agS{&jU@8PK)b)cSM~j#3wwAyJ(*HEzn?d)KMHwF+}u9>yE5jsk7&Io2EKp1(=v79 z`=L3=`t=_afE|I|?S??YHP`&#WFO}Ka1%ct9_^7gI0EFC**PFQ5jAS58*b@1b*lGIu&lOsEhvx|s;bO-Xo=uYKOi45y0PB* z*r^hp&TY`sU2VrVZnY$e%d=~ooHjY>n=o#3ShgPU&u<_%Im_BI@?JKNb_ZHpq8Z{z1Tx_XdOt36t&!pt$&`A&!d5e2`; z>nnEJ2|VPz|0~Y3rg}t*KO^f{hStihUs)UwAb@Xe`sSINmowtbyz)f-H_vdVSXgAk z$-ubgu+jPstvbXiFafjxU2t^ZbvYQiXD()g9{)mp;$iH%6`Ci@-=rT8?}B%$|F7(e zA^%AD>&mJUJfKZWmp`yChg9X-)r)$+G zSIwy*E8ZHAA^bbwm9q1gW8d*3_CBDo1B`k_n#%Ea?l^k|FcQJ%*Sz#B8vYAXaq)bQ zR|tH|1R$8eqYA| z!5|~Go6!5j-^!n9igz)w4f#E0O;ye;$AU~rfbIMM1O8R!ZOHl6eQZZ|$(O(G`Ny_H z(FOV1LKETS&A@%eRV_;a;J$a+wO*~N2g8BXgepLuhQ#@bZ|49~hMz1RHfw}I;A-lM z6HFzl@~ms}^$X1Ul);zF`is0(t|~8~$c&mPww^`sPFe^D_msKf*;V{EA)t#5HQ!BXs zo_S-uDbfd=*|(UwGQiid@6L^3BP}k+&pw?y_@dU}-KgP(r;oT!B|XuypH|Ew92R;{ zJ)9q}mCbTbOX~?cGUOwT@3kIJzb<;Y-6jNI-$NV&u3$)HL9KF=)@%PPuM0%qm)C>4 zuh54)M<8B78v(D22?5V{-}ZYWNndRR94q#kbdRI(H=38bqjt0}>GpfClwaHU44PP6 zwi|YqbNb(AJ4ZWW-kp8%13KnbmPZ|Pyw**<*1fKruiQ@bObkb}5`6v7AN{T0a>cHs zVoX@GlPsV|mUlB7%z0Lx2E-qCq`+*w-m$to(}Ju2pn9587QV8bW+HRho2>Njz>$6B zG!pa-^4xat$!szoe_7Njhdtp@_ZrE05~_NHd9GXW?u>0k0PyZx0{Ws^kF=b>sne0aFP9Z=C_poZYuK2fCLxf{#zR#U`UOY^jI%cDKkZ z&tYV51Cjoe0#p0~GyDRxwEt1u0<-@?ECRFtLBayF|3TCOv;RR{0<-@?Oain2g-DcT zt$n4feLs)?r6QtI_}s9{RWvdM}#7)%OKvm7Fl2CK5WZYbHMkaLgRh{s~jAnZh8zGjqh#CrrC$ih}@1 zED`+8By=;bnbIHt5le*6C;WBIlm`JASt8m$Vb(QM83YhuiFo>iIoC{e5J2J6(?UWw z@0zI%0vLRHK4HN%Qy&Dd`}BOmqHE?i5Wx4-^9f6?nE((V`qR@&LbvRiX$b;ke|kP) z#Wm9w1gQMe8ReGrq5?QpPo0ga7Mc z2e<4h`LDO~(FHvt>Xc7c^;ctEA~xnY-SlBa#s7`<{vRwF`NoB1Vz7lBeev@l`9BZa z!t0;)#(4R}0DWFaf36R>kzL&4Ip5)3+~GOj!eS0H}u;img``~bNnHKPSuSTti1gpji`d@>|%O_9fzC)auJ zHf7E2|KfeE#3s1ixxf|L2|oLG2KkyjYp?n;UGTo33v_p)AT4Owe`}2AK1kQ`)0drU zJp1M>&Z|_Y#v$IwS z@aqAu4Es5}2PbzQbD8B|Kb^Hq{mknedwb;}^Jh`YipwOLKew99AZ^obTsLct64kVX zETkO}$+}G#hH(P52 zw*NCge!x=udbO082E96+Rp~?bd+tVrPs60sBQ_IdO=kQ_yHib&x0T!M(c!Bc1D>~R)@|xH_5i!VE#sQW!6eAGXrvB{$7qV`xhZ$ zxD`~Wpt za?C*pvNr;Ern8oyhI-xBnvSrfFklxk(1~~j-w^MrTCd61y3n*@8NzAErDXe?a?PH3 zgCiD5|Gt?BL`M?f#XM$S;~2LZ=s3XTvL2>?(&S~(@?zjmO#Bf5$9U>}B%jb3Sy(n} zlN;nljvQcEv3WE55 ztu24m9bo+n5T5Geuoz_Jw0Mud;`T{LRZ{x3Y;0Km6H1{cbRUfGmtF$<4K6&6=`=Po zb>tOI02a;+Dv68Kv#=%3lp9iOa(D%y<~Kj7Bjns^D27>O;FS9LkNYl2Nq+bdz{_6$ z^nzNi6UIX*^x?95*6@W=3hoWGc>3VMl2~HnHTCO|A#e@!ryGn4WmUOeLCL7>edh>G z9*gi7d13Tv9AxK(O4VP6>LgZ= z=?m0vWoeE!OdY4!IiHGKQWWM%bq={L2h#fH9VJwaT_9;SjX@quG6{Tx0R*Q+zo&63 zD$289c<{Z(NlPL8t%G?1R%ZUGF%gVL$fDq_#HL`vj&!2%{#$Lp4 z&}vRd`AsX>1hI7d4Rc`mC=O0L>e07VggGt<-z32yK5n{&qcbba>L3njGFG~|Xgo)x zd=0EokLhF44lN!$^7M;VkL1T>%MREfG@P{nK3LIJ z=HB(|;5Ik2gT76=V^cv=lP@H{p||Kywq+=GbE`{($`V4Q0YW183e?Ha@c}={rZHqf zZ3}8#mDV^6w?L~V>@@5=7G!Y;(W5Lx7azBwiQYJ+S{KPD@f=e;n46Az7?|((=imKC zR$283;ZOKWmYw&#CI1-`_DqFP-aJHl!Cdbeuk{Ps10AXF867Pk=NyKNc_Ggurr1sk3t2{Q$jr#>Hkn~0}RBi4ln zU1zsnNN(NO}Q)(Uz)0MS< zOi&Sxah|DS_E3KI3iaOl)C|XKi}i*l;UUSOVkP2c8EIEA_7y#6i?o=_LBEfF^0FcO zHWmEp0fXA=)nPiPLM+ByK4k+pg^>4U`*1f=?QftvWLtAvP)TTo5Q%IARly$niUXzI zZWE~t;ELsvk})@{rESI(DCoTKs#74bphnRH)@gsN*llqX|K zBb9aZCnr|7_)(twyZtBeYH3cp0JWmxAyQgQiKkd? z+yrt7k6^e0OgMBYifZ&9Vh94szcLe7=o~aKNVcT_RD!-kV0bZPNUk~=pkg^dGGSZI zt{e(d%iN~T@(}#IpL1?t=|KVZMkr?&2}k*1ZSBX$=~C>cUEkJ>(@`?kq1eRrmj}5L zk?TW8JUxAmgMi1wiftCBmGe}d%qkP05C`}5YJ}OHgRx(UF~~?}MuulG^(_UGObid& zip+r_H%WP^Hizy{kFunL%^ZO_!vL{SwUq{2E2np-N~5@~xy|m3MYvb&?}$Z)=#q^&<`$hdV;Bl zXZQnox{|r3qa3^~4^*Z7e5<^Fj{e4&d z!8B6R@Ss`GC$3_s#1=?tX*f~FaPv7cUB9JU)`JaLzkqprVJLC&Z?eoLmIw*Ae|<5- z7%WmiZh$xOq4W zX4x^yz2*mvpi%OKQW82)M0$B&?g|}I<46^Hhuhp)rwOi;8gqx>^)pjaP}w8v(zS0j zV8T!cWRWT1ak_KDhjilG6jLNNaMXvunc;Be^9=p8x@^VEJCow;C;KOb?u+UQf*^xY z*-<>9%O7372u~P_*2-g-&h8@^+3q)o)6hPk5}-KlPuw1ZyPqx;$F?dmE0<;6#w{{~ z)=MN4ujZQ+A+jGHalX6OC`KSE(PP)vejJuF7YH%PduiG*O4yvxV1S5WzASd@HHm(T z@2(-ajBE;~g?&iHakUB$+Rwc)j7qYf>WFPO9kz9E4CUdhH`zFE2P})1MZnN>?FMXpYiN9rO1l7yv=V^WIE&BQ^KTqmPAG(MK&uAWwr1TCcFUsQ~>6T)??UjHAG_5 zxhOxDm&~mw$dXNtDdv-B*e+%#_ZT7YAW@f6unpxXLrZ2>1>nB46wT_MS*k-fkRjQw z#~ihsAde4*NQtBH*kBb-(+Y`MbvE;l((oOp(af``TNJDqS9e%Rt7==VqZM+>*DjO6 z>UF~KQiby)Q=O?;?L!%-vI)EE`8Fma4rmB)FF;4b$V~GdqezFXpuN7uIa#LRpKt%G zuJ@P}GYs}wW8Sy^yk!)pc@Y6743D85a>MGy7AjNY$SEJy;xuVue3pW!v`typpCyDoOud&7_#Pt9 z``40ltv-#v$zwK6#U`HqPE!EV+pNiqr#mK|`4R~*Yb%AS&gl5;`|I!K>Q0@sy=9Uy zW}Ah>y75{uQI{Tp6%NXV=iV^qP2R$v3?)GM|pxq5R+9{hC+S(6OenE5J9X!btO{xO-~)l9zYRtVs@)z3BQ zjRDWgCq!w!g~je?v5}o+%g=St{{6o3e&1uYpOL833U58dW)v+C@;lIK#o=dmlx`yr zhmVzH^JchseT_8-=z%XCW$2FMcznuF_dO_vo{b~%Ujq^F>`{@aoc3A#{xE-xSb8>#sD zlONY07o24!p*!IMq_Uqssbg1 zhVxw>r5d22KkRQ%pGcCDc@GB98G6Ugd4CrBrz~+#kOdux6l$Z7B`RtixvL zcdQ0!&RFwrDu-+rQPDBoHH$Z6RSgLf+CN}axlCEhS?u*QDl!uq(64fiSYEW;87-S> zX;PaMpc(xIk)`x#(`?qkfn@9TqY|r<&7R; zGIi))Rk0!`ygCt0V(-S4IUP8^mE_2Jf^Jr2&H`5a zw9JOqdYeNOWtMPQ+rE_O`lY^Y%*v#kgm`mRXbpRhONqf0ZhF!Wkt$k9_C|#*)n6eh zzKwIt#G~p-LH&g~(=K#h84&P2lZ+NBf!iXF+p1U(@hA-||CN*|VUGfz+6Br=I@ik< z<;qUG3dE7wA2wX_mcctD8x#nvD=SU02KA9g8c@suGy?-i#3&;i4H{GyC@qm-V}bA? z&WQMa(?>^5$|?5H?xj4aKuQGT#Rk8xh1LPHu0`kO7Wnt_rmUPqo zbuG+N;bpG&R#vbXUs)o3B-7g=UjZ>V3Es$7*D0J6n{`Zas;`i*{oBjMDy&RWdw*ZJ5s- z{5fnkELi?(L6$`(Iuwng0+gx>VT(}@ByR-zb1@DQ<{*RV^O{n3IWOf+3%F~PO!5be zlZnnSw3wqrew)-)_x?Ost2ihI#EJ4>l@&mrNCzq6_Rw$#$>=Z3dNH6zC@j_3Jp<7{ z)CrKpX4PI&oY}#6_{Mmaf!v8>B0bPph5yWAsiF_`9Hws~<)U5CaJZs+S1Qn(m_SR zgnY4rhKlX%yKWPgDP~Jco8w+DY%j;00CzFQm>>1UE*D!ttT!EI1g+$hpV?Tn-$p7_V)ZKgL?lMX3{MEf64EX?v+{XZiXbOp!u! zdi91-gNaA$#t~rcq5~Y{U1R61 z+2M9uF#eeHG0VzlMw!8cmI^_HQmg%LK0(t*;r(N68ARZ!T;shvpT@?b6u_@ ze9PHT2D>%+hXj65E|t5Xh2B&y77%I|14X^dx=%ANEZo>uL+PK8lnpF)SN!#bi`Cj? zece-ukgOB1NKm8QjZjGM%yI9TjYasL(*~hDw%e4ODq{?{e6ULa#2H-Lhg2TRriV^q z4ikP$=SWTBjitur?QVpo8VSRkdJp!Xhe#UMw$vaBLb#cU+zhl^o$b}Gj;3XHU=oH9 zxJ<JwcJXsZ#g2Ck zjL~pZ&cj^cVJ|JgdsDMB2NPfAq=_d(ssu<2|n?>bmc; z{ZQRV#rr5pog9%C6%z}b+@V^HMg?-?0)Ds|{ge`BgXu>$k*l`GR;D%xA@k9ET#}MF z56#DI1-EF#Kb(x^{JOb){QCyM7-3jWiP@$z9rr|XljBJ7D`Ri^YxxICp%&0vb~*Tz zoiJ4=*!*K%LTQ-j{mQn^?XRF42*PucB1b$plCn-IXsAfmy)nss%m&GMU&x!xFuc-H zoesM%mj30xp(JBu0zAe;v)A_B^bi9_YtFenjO{7=mn;fx`1-WQq4UfDSW{L8=J z#Gm2KDpa2?D{u1)zXzMwnCfon0$*p{^;vMjB`0uZ{^5j&45g0aPriKB@z<%tsE8bz z7%$M?Mv|NuK4sk!=7T{xuSR4&;>@JY%vh(-%fBs!Lu=5&&1yA5)Y4Mt>b{(u(4vJb z!*?|9gq%&zxKwj@&Y$a$CZg+Sk2GBF#!-w!UGdm*1rLF>yL}HbNzC;o9YBsLm=Arn z#~T||V?%#At32>TBT;rjF#ohgeYMpGCQ@R`FXC7`EBF|e^d0xGD2$112s)*M8s(~~PsUMyaKe=)HbM6~i%;fQm^2g$7P8VxO&*lR0|a`l%a z|H{!!K%Zvm7lp}M$LbqRDIKbpK}qHzi?}*Q3M&+GaVpK>QQyt5-ZOLRG> zFDKd@TQb)}%0Dc6$q*kcP5vRfqm1gbM<_jLDFP8Q%pynO9H?R^vGfWdp>>H0iwH^6 z4lcn~w|FksX$ljc=MFEer|5Faf`6S>Ups2qb)yTmh~OKFJ$0rKGu*DWOP?soX&b<> z^y?NWF0@G=0`#s!6dt!5< z06|)*fOgnNQm4>r{SPyKL@KaJK~z0928v1&XN-bUkUS9Qn?a={@r*$IZ5f7i%N%M8 z3r$5SaihbsnZgdz?p64{t1xw>-!RNyzTa~w;z6!P$%AaFTLk(xKyxPL>v+ibeH3rw zYLFE5o@%suM|x9eo|=|&o_Hh;xjg%vaMx(E6-(%rLI|el=xKI@qe?16drE{F#zoCh zJlYT`G(s^<;qq6-aKi9W;jT;1`7~Q`)t(e&3oMXD4Q1Hp(ud0*PFQO>hIty=OV(tX zACM{ji}&Iq zLV74)Qou;?+D?5>Tyy4sr15l!k2&;Q3RcwW=<~!k2GlO92^zJ>9rZgjH4GDGIZ+Ky zG&GCJVZ3Qgrs;EoVOqLN7$Ro~gdKv6)q`N6A)A-!nL&wEgU=o`*gslpXet0k8wH#s zaNkqeX>pkdAQ3(MyjZdUB&6RH8_7%?@;l`CqBm{TwDty)D~V#tjp)pNO_UCs(fz}y z_LWK!{sM<^bY7pm6w!joDlfI~(7+n3sHzJjU7oQDnM7S*MVGkM0-BL!Z9vVpRf~ClLI;p^diV%m2U-x|CkMK zzQHV|Og7*Z(+)1A{}&KpQ{|IrLW}T?rkE_4q}hY{4g#-2?N6Ykz_|EGgA>w<#kFA{ zmA-~}7mU%eF$~;?%NbLyDhiorhtTGJh$H_+>`F)Y5JZmc2~@fv2`ZJ(#&Zp`W@i=t zp4bXNWchv%Fh_As`jQ80c$-lJKffs9bSDo`#U0X&!-I;`6=<+(O04%mrjqN;u7ttu z)|yHcoI&6aW2CAxm1iokrjbD)(q0MH`P}dABeJQ7`&gBGvj9`9+pncuD5kWs!FrDK z6DyI`Df$uNr}FaM_NcDWduQkm2Qsd-U*fWnkcIlSZ%poyODZ(TPo})9xDHEpheKq# z4i1}lH?@!E$I8ndkmO8RF~5_O%eW6~b;1Irc6F95V7R;xA8Xgb+)!dVI8X8Z1*M#? zb;mae6c#u2Swv#yGUk%0H0PVf-%)Dz2Lb@fL)oQ)b#f~XK`+uhMVq0nu0u|mw|F}! z=xZ|1FyV*YV(=fBuPjsU2q(Xxo1o%OqrPf=M>Ao8i!d$%Tu9DsYrGptge$-G1Sw`4 zbLvEil6Ahy^HmNCUD{;=*YsE#sr5e+T)A@#rM$7b_isq2;{^ znPTfwcjH0aczf&$UV9^~`^uB9lp2SERTjU8MUqWATw@x-r9fI4IaPAru4=!|t$)}fW~wozTU4Ap~$Dv9GEj$;#e z6AK?4_231@OwIWv5wDTP4w$P4UpRg>dXLA%dEc|VdJ*=2S8ciM)AERITQS#Q>R9|$ zu9I4Men$2SQ*$gEu7mRpS?-vTiCOP?1280h*ONk4BTr=b*os-r?!#&!&H}KSkdTrY z_ivO6y)F$-7DlPZhNABT$3{5R+nMEidGKq9BP0$h$}RpGM_t6z;tsOWghL19wV4rjz=dR>P1-7Knvqd!u`zz0}S^u~z zRTIaR3N8)~9nz0=xeVO>>BQg>l8HuaU-HXLg)Xl!`(3+4uMM!dgx6ezs9$=TPdx4uB6bhM6ZwttYo--;yfW2PZwGtay(3x z&Jx2*wWNSAEMD6wLPZ^5{Eno8q2)t2L=N*93Yw{R)IwJ+F_T_d)&5hc5URV~l{;vJx8*}`q7~;Do{m&%cliM~ntoEH z53j|VnoO5v|D^iHfbtei-t@p-X4IudeXq=AwA>D>NGntzwNsiUFA(imXOOurOlH^e z&&H&E(YGa9x|SxgzDxPi8m;RD)C$i+AhP~CTyioa@s%Qe2i-|e130(eoNxbIW27KY zakDu#7_FZG=e>>cs4U6rWMWoqL5wdHy-f&+x>b@0=1$WntyTm|RIdiv)8D^`53N$% z(4|B1dv0d~@5NqNi801>=i}W$6QaMhp@np&aI+a#B=%JKWb+aDzHIH30^7gY+>2qE zQBIr#sddbSOh+@i3{c6g@_U9SAhsFG^=3x4g&dxmW?nKf7xIu!mR>dh{e;LqVBGHP z=tB4a)n|el!9P|>NWR~`>h3@=cRpwLxuz~SiR3e5{S_%)qKhl1vn{9!kM3rFOPflm6)EgSP!teMMPkNU)85u%mnZi*$4xd@$lhF?>q zIGCLrIr?U4+m+r_$h+mLua+}3l33=c`c4alOsgH9V0e_ag@v7{s#BiXI7mI1yH<}5 z5w5I#Ftk5VJ26~g%_`RAUeq8rqNcG8$iV(~GoGiMt;nawXyjrHqSir_?9}I9@)733Hfems|H)A46~b07#A@NwP)#(YiZ7 z=6x#9XF2gaCuQBthN5mxwLr7+jR?I`Y4j$x4&SyqPCF&X>UPek6-s>Mjr9hahs#?n z-}!bt4*gsm7dx54)J+jk!O*tVWQYX44Jjqi#o~nVp@pCB72w+H$?jP8?(DQtwh>aQmYle z6gBBn$6!=EM3*&thgO3XB6dXjN!crN3auTMYqxbY0YhcLzT<53vL3eS-dm%Cu1H_y znKZOTUNXJmp!8x&l%ZCH6?D~}hqFtOZciK=p5H7)_46d0PSL#U?@$FRpL_zjc39nu zmEeiC-Ps^Goz|9T$FsItiN4OvdPUaF%BtZU+9n5>OWA zYJeIrW(4KcW|rU1%q11T*$l8EDs{=GAyZt9LMz;!leKOc!BM8CYDw-KZ@?K0=~$(G z63`84D?w?v_#)N4ZNJf|8NcAkb(W%p9Mr(job5!lhX>u)+*0+{OVNxer92<;ynm!7 z_aQ$LkD1G~PVIzOx{+}$&|3nM`I^Mbkj&CuTESpn0?v!ThYHgQLGS?Xv42+>P&SYt z?%B)Uc`*3xEWz3CHua*)ER>%>l!0=QenPBG%Dc8X#&k(YOX`b0F@@po7iEZ&a+frZ zkIDJ`w1dvV*gx$H3euB@oW_T3Gh@H_Yd=RhCT@@Z>MwCIJ4pZqIzGX8YRN%ihYLe2 zu-VL9`KO04YAfkAtGhK4W+;n(x|)X^|0Nn*|AOiR>30~EFsp&=uG{s4v||f4_3Hei zn{|g-L#Plt?Nb^EZFyy{1^9kT>WqC{)GV|~5?#q`7!l3({A4*YwababL`R2pZf3Y< zy|CF9Jhe!t_}}c*+Kd;v9sY7y#A25eV5-|*j})Nh9%I2V2$R-;G)BXF;XY;@D(N@# zE&C4*-HWx_az_(2*c#iRWf$&4vX`f!T$VAm#ak&b)sCI%Zk#6FyK9Sb?=l7V7qP*y z0X4=AXw*v?uZ}>4*7?B0wQFR;o;N4AuRipFs0`=qY-UiAb50sl-K&mb4BSz7(&Voo!M^GbDm(ETH zt@uhG%yG#4f3!g@5$KrrqIlr!0x9bacvCG*=%EQMrIz_2e%V&g6*WwrnoAnzX51D{ zmF>2vv7VP|@7q)%O6f7{c^D_7GFx3d+zgt~SZ;0vlNzE#4D^J33SU zgUMxF)jF%)YoO6K$*U1~h%zggNZ}_i{nB_>&V!gSuC7u2M{*sl*?^XzS((fd zt16+MZUbKJwn%|9GG^;o+}lcx=R!833D9BThLQ=4fHv!7z_&#JMy~zDJr~4gKq~ts zCHT8rDS4P?-T;NCd4|PE>f`TP7sn&jL3Aw0O80hNF2QhE+R|J5>px+W$sv*yt=vM- zipzm!IJ-KL`C}_T*kz;OvY$Cs?21Sg{cUBw^@H+7jTP1qQ2wl$MLBEnfQaqosM&J% zZFb8Bs6@%&A2d2aBrVLMs(n0sV3|2v5z?*dxZkyQ^{|yeQ9vDlvyG{IrGLKI@M8I5 z^8<^>)YI20LF>(?JR<7^_CCdJbuu{Tdp1b6KW@?cyBNl~`YqBD-pf!z=FW_F5v=(Q zxfq*!Kk=r5a6|AD@UvFx_AOnnil45BN7o0NqPwQktPBKsUsv3%2eH)toU7Mt1`xMR z`1Vg=i}lP7&`OK!LdL!6RKzRcy?G6pM6o_i2SrVjVQ#v|nmja}nRZm&y4+et>`v`l zI1t?x8I=3~d7(DpZ%xG=Oi*Z_!R*M;Kp%OCtQ`f0aeRBKoi#_XeImhYCtT$gI)_so z^Hrq`x?5H!_LY^vRW#&`DhFDB|6cncF`%(b`$h8%4Id;-J7S|?o2ggr2gU+#O!sRe zkRwVW$!7mU5~WHW9c_g+G!P@{$T;rg8OixW)uJWjOG7xKbBA7pCBl*rNvRR+-aMz_ zdW06^&!|{Tu_oS|P8V?1-=IOV1*%|-;(S-p#Z;KudIAmTXelux(ymrd{hGicSn3h6 zOjkQGe)u(n)j<|Qq9`R*S(zxc-RTMo9$)wQ+Zf-(6y@F&-v{?F!MO)6-5V{Ae)*>R zlv-8>&$1X6ZL{oBghY$m(q-EY51{-5aYF}J@l3;V;1QBn`mfD#blSzVmvvr&vP4=v zKPVLD$IF$J+D4{iOE;IyRrtEhnEc$h^h;o>eTE&bRyQnvR8nGru~;_(efySFm33mW zqB=uCbCHkY5(5PRtM5`Z2#x{x%}T+CSIQ!x5W6DnOqToVtA4Tz*`0?#kzZ%6TLXyC zhX0~nKWXjjtLyIxUV;^%!gLqzAe(c)gmMgICb&b43x8O|-v|lji2p7;S~=95p0|hU;bMJ&D&$UZQ2m!BGbR z8Rf81GL^D58c;Xe<03+|s{o03QT`zD*Cw{tS6qzXDa**xM|fdoBltx!@v@W6QdHaZ z*bT)CpKW5e(m|VVxrPrG9NXjLV)D(FRTi)jzz>__0p4G3dCKyO?MvpH;>)I6k*1}` zKla=6lB+eL8h5qawOgW(8s}LWE29@|Dt7N;C22*d9xDwPYqsB>nGExYQGc;xE3Hv4 zOrCI0r;U3z8dmCM;r#v19i|6+DQ}!^6TV#(-d}a4uiVT!YSmxll>*yTX=CT>2!{(E z+t1$m0xKPvx^NZ@r$oOF&!t%&BqvY%DlHJ<3))c6r^IcoQlIH}dU>J_Pjqp*f0W8G z6L$Vw4;qq_qg+t_DqwPOc7#I3Fu00#6JW!wDN!u4&SaCSIhS z@r#J<)oPmO@jU`38|jQ81S+hx8;Rq)*wsX_!#@a###C6IzDvlc`Kf6jD}e?*;IZ(T zT~b8d&Ms!67&#=O^qrGlER!I4)A=tCflym{YO+)a{1ZC;l+_P6wjs^$jyCOB^GiB| z@ah?}K_javl7X{?>xIIXlCWrMLyGMzxHIq+k`G*Z#ehil4Pjm;h@G`Ol!K||#M>_D zypgF<+wHLn&YC_kAYkE z73ZGyPA7sWkzTWxwpIiTaTTqGsr zOZ@@9xJzC(-Vji!KjbrJkOK@uVXD9gt&(4$w{a#69~G`gq0%E8uU zo%qC2N9%9UQ8FXl(9%xqodSj*3J-q@0c1*77V2NaGXcEO-|Z`5!6>EKz1A>bi+u}f z0s0-;6}uqQdo-&x?D%T9YwDPa(`6_^Rlee2EPBmcCCdcpjgCkHtxq_ie^JtSh|ey`X4j+;_(dPG)yS3jB6d()c}bhXCC3sNGT*o$2a zg_HR8&DoxdPP;#J#^&=ja);-+{P(F-#tCg^U8UZ&UHug<{a zo?u5{WPb_ul7WE2LY3}e;Uw;;xeQN%G}!TY4}q`^ym4`u?%04xQR-r#g*;C8n9!VE z5N)g!Th6L=KM)OP@~#IUqcTGDnF6V{Zs{Rl=xQC!sIH2Ucl>9bn_)~>$;#y6Dh4vH zwV}w_K;~U(kND5J^CgEtE7q$xv6cmSgTY}f?e^&;cWgjYyf2K=jkR212Qh%J^! zKxkSLYwYwseiiU^gH4*qGqSzeN0++94hCSV`3npi3PT9~-|3g?;252!xSAU7S5;-B z(biNgDjQ*frpZJaV#Jc+8a|Ea(y}?CCK*N$HJ$*F)#b^|>{MC^2FDaqYPu-{IC&^L zA|8eY(>us57dsjYZU-wv>rJ9>7M=Z)v_#PUMh*XYwgFY&0w3WcHFkw<-wBYw@9A1* z(1en8%d(Flg?1AS3^jWq+1A#2goacVc)6W71v0JN9s)zQRr6DG;A<%4X`~f!!q6bw zCREcKTgk)m+f&Yo{W1+hni$ds zBU5v}M;R4`2KjxVj0I^p;zCpAi7Jm0#t2`yjd;l-E4;yOsn|}Tsf91f#cae0ptn3p z2xFT)_{@luG2mw)Y!td0im03KqGhh-9)rBtg26z9VaTHV%?8vKV9{-VDh^(D35?-q z@uOE39jaQ|?HwauqI5ruWwNHB@+g-_FlT$zFhsA*m=}ae(@#d9?9@@q_CAYKS*fx; zpod!2(6Oeg@|?_p)tg%!gydRpdd&df{H@f)15C3#lMT_N|27hBl^s#Z_MO{)@Z>Sk zSW8KJ`wY$4TDRtu`UIm)ZZ?O;hwZpvrdbRQGN!Saj;Q?QM-wdP-WCKS4>5%%;Zq9) zY2ua+S-agkK{GYBKyV~l?6!&KIifHp3(lr@iz%X}0V*VgCi7jhI{js_&TJGPi2t(I zY=fYpQ@&=t5W>b~P=}8q9xV&gcwqUQD$c}Nn<@_SMB)?~=%6yj?dK^q+vpz8j6w!z z&`q}&B?OVa5v8F{dm<^Z-g=v)bPUTu9nm_)YtDqSf(6~ghlezqt$PPI3<$v}z`ck*DNDInfi$lY|r=FT_M*Tt$8uT^hUq&*04 z#Ap`wrrdtE63%R?#36MiZahp3Cw3NrAmbAF7>`hpB--$T_Db0rN`Jgc17iKXZqUf^ z5$PGJ8Y_hm$j>nLw(xJWwSk_!@+exQXQsDfV}QdG#Bj%nNJgd*6J(F>`ef#X+{puFs6kr3vXSkTh>>WNBF0+^6-$=xbx@9rP5)!aqV^(JaG6;4cSfkiL zS(ml%iDz2RY}kTcJc2Ks76!2aYq!==a)eQxY?jpu4dEmYM7yGkr9-?e)S>&j#4yfV zHno!Cnyz#pF+=U!NUqE!72IH$Mo>h@nmN5cmrd$o6_BdW5Iue9SQ>#M5(hG$#-sPRJbWd}<}^-X{YH{fF8}t*(^;oylb50Ui7@ zDm&97anr7z)T? z*1Qx08oh@aw&ffrT;Jqvk$9R-1Il>%Z87U|sZxnx9)vfF>pF80soi3Ksv!{Ybl1V; zMAA%Do1tFsEX45(B__hH5GUWhq=hLR2f~4Atn4uDs0WPc%8?byW|`jiG#Ny?<+ z4QAwaC#m!DwJnLcX2YGy4)6whK<|<-$^pj>W$G=z2Qn7dTcXJt_%~T(275s7;=S6@ zUUW0)@m`@z(bQX_G6j3`@gE=|Ld$Ynm|5?cO5{;j=q;^D*;h2fsuaT+_*>7L*R+$+ zne;N5s8OHww&+@Xqc)<#Ivv1hV53Gv$oQG7QvGO(OnYNsBfZ6GdYqGz{`O`}8G*;GOOBDVBtS80LPKx!yT_%TvPZQ$6FPu-1? z3!`nf8k5Qq8Kotihk-{H5JK)~JYqr3C8(;#t~<1nOKnu5F=d@ zlx=pGA1;WHvnWdHPJ)sxyh^vRbO^(`DwxS+ED%a6Qnj~Oew+gi?Nf27Y__+J#Vk$1 z;lW~}m~06SS7e}6JRpR-v+J$-V@o}zfjk~;=wE4pl4855av@Iybygn86 zfCe>BV9#r!yDJ-(SjVnDzxy`vcOA^W>Y~Qu&CxTVFqmrM|1A^LJW7;|&1?O;XuaZs z)h0d+sz8yTtT;p_=q(x@)Mk5C28{${Y>t+c!wsiSE$8$AzvL@a9t~zA-Ve3kBA+){ zG9=p+B3azC97o(hU#^L$`u8(J!uOYayhxpEQ8$Zw<_b9Af&D zqepZc`QLiA0kf(Le*id{PZ=pDnU54y;Tjw zX|=Yip5Z|?HlkCz7TrEWxLKSmh2*i`K8#2~oe8vJAGNblT24o3D+9QIXX?ng&=enl zrbks#Wksoxm{3~&3I|g_37BqvZ0Q%6rz$KM%L93=8(Vr)}YL?nzGD*dUhO!c0mEBVFlL zDWDB)B{nhq9ccGzsJC1RGuQ=wU&-H4wfDYFA+0uXu)4Hg^fY0fB-W46lf+u_AgB+e zfyyUE4Xnah``#)CWjsxE0t@&Hu}T@tPVG2KguDPFPd4+R&r5=Ju|kNGg^q5hVOm!P zGbjN_-kgu;F)+QP`hdnI54Bu!$!Ifr;Vq71L-59AQua8K)LV4{ZDl(xq?eTY!=d#L z^GnT}qcTo%qm5o~Qcq`Z^+vgzpZ8%{X}0=7f=S^X$q{6jmZK$g+52iek)JXoQhU6l7uYveQsis z+Pq$g6dJuv0Fh{GkkbY-{iYtXpvg%oNM$V0toMji6*Jg~*UU4_eGqGTG2(#IV#`QC zL~4v-^DYResj%m=_1IlRik?=juw@Ui5{_QG{j${3w6~uP)mjC8|5&^fUDkk4_fCHkCGS z2?m$<|zyjh~y1>CN!7L8+V>qV8rFPDg8o-MeswnQ0lp5 z;c!R^_{F83HIU(fJG>2`?0knM^)^kp0b|l>4Esi<2&{MZg}v%+K_|a(v9ksm`0&Ye z4l)D~L$GQ%h=585lTEhF0cg?FXa5U(brOrJtJ9}|YIaE?<~@1FR!SdAJ#_&q0_T0R zZstZUWuBC`BuNeaLE^}n)XgA!K`}Das%)K}w}cT+#dU^5NeT1V)X3i=)y_g77qm7? zi_;Y=4QB=Wt>9qDs-|;hD%B$J6DioYv~|myds^mF+B&bIa1qA#s&en}(_K??kHq#7!#hD#k<&6!kZAkOl zAPP4mQ_2|a$63Xyns6PHf}qM3J;y1bggadY50gLd+I%W!5mQwURUuX{#Vb~8rz3{A zs(X()O7I%Qd4ylGi_H~9RWX7oO>f(AtN{d8<*i=bjH_;mBitN|tEM*sn4$>SZ}zW>FL(@L6EuAWw5A9 zbsan~#lQ?F9()7-W<+dWc+6;ZA*}M>XBlFxenKf)-GCNEv9X_@f%0*j3@-#m!q?h+ zbjMq>ONM50+!O7jftkgP4*}O$b;B;jlsL1jqPNHvxSHiMm{LC#Gd1R9H8+zFSX`U= zw*5e9bH?qxwtH#~&)luu4;;40)ZG2rZ##g6#ztobd(mF`0uynT-<)tnzcZ~%#W*D> zpy?;u1#xk}wFGKSWxv6M@Hp)yj}llr)Q=TQ+B>Q6#MCCNt!n$gTo!y}oq|2}nQ>b> zuzM2f^9TC{F(yF2&B1=G#4!~9Pvg%FESiU_Q^-d3HB7J0kWlTNn%Dq`6p6B!8$wd) zZMLep*L;(|MXC!W|0Swp|J5O=<#9=7ARgf=`j@}y#7!;wI{~(n`dF7MRo&j?!H+Rr zs^#8q@#^GT-R)Bv)Me&ERkgO-u2jW~T51R+tF*jOcDfK04hNY+FJIj1>*|%zh)t0g zNa;#+SiuK%pD1CmQAQys$n|O=7x$ss1B$Iw{n&@eb0m=eBKl08kca%2GMih{jIQ|k z;fcj03q|RTsruUC@$=}q(ACGA_P0_au7XnZenwsFeiu11Ddp|&h5(D=I&OqQ-BO+ko^V`d{XYG!Zz@>RE?`gtH`$+zZEvZKr`(k>_eDoqL${+M^3LYxk&5-1zXxD^{VuJ%_&)h zUZQoX)HyUSHzAAUW!YZLSa^z$p_O{B74k%n!Ihr~#a_lpb7jwG zV)klrO`L_27GzqQOYK*Aip29ybD^VnQ=Y3h@R?7DJlFCf4Wis_ESq@(C_*izoUL!_ zATcPXEUV#2-Ht_OeH-CQIq=N%!()J@97FP0t2nP@wstc5su$+)8e98x<}V~_7Hlb- zo+yp3?H&O_9#~Qdbj4+%g)$`(hXr40!Jsuh|2UQ&9#lnKsOPj`*PbpqM3SF?a3wyF0rxR2-tnoDcGWXlh+w~2CAsL ziM3}y#>WwlIii~#vrE zD+8Vz>oSgXOY_cZ6}` z7hCbfKw)E1#5;-JkSfR!@7!a8Av$fhUy~&y&8=#~=Db0XIHO+KIzRR1V*lQYn*J| z^wTZcFmI=iThf@R`?NMB$=OvbPNeR=IXKQVseLnL{{@D1X&;z;`%VQa6cbt8184wb zg*>;8+Pub?2MV@_3fRHE-h6A&FKv#H&BF%GjDmO)4fMBaez$ZX79B!WAuOaonOG!n zR6`-<3$r@ZEkA-2iEkW&ZY=MbLPB)rTHeP5oy|HaK$4qk>iqGpwgn@bD=9!OB9!Iq zB5*3gMzob>KLcSS&meM12|DGFOKeJtfkxZxAw-q7s~4aZtd4~2Gj#h`Y?Q`O8R&{2 zmK0MNbi+g)ttS_758>x_PREjJ*1?)*QIM0(SAVXmWvkFxR8ffHT5N;~IGY~xo#32p z9As}-SNzgF18DAs>Yk)_Ejrk5UYsPi`@s7s9^yg3hXb~!qoDe4i%DgxdtjRQ&LOTW z-JoK!>BqZu9cn)|`#8ve;j}PS)T*eFK@J$d+#L(4F}%I|2Z8~r0CAIY==#q1NfyO( zjioB&d%R(IE%(R7dmJ8kmH^|n;4=!QvMx;J7Q3`S~mrn@#V zg?F`9wrVs86b2Z~T0=t{^5&@QjA*DR?Q>z7p`qftTHhkdx;_1A@_m#QSf_#8>l^*8 znxB$~zz+$m?7KxmPk`d~(tO60t-6bnsaU@xAd-vrU{)~0d(Ld-V3??r4zp6N4TmLs zJ{`EWyAN|wsI|QPcVdxX+dhYhxQ|F{Z0I2%{p14!QK6?#Y3e8X$xCmR4O>xfl7&?&fS@LG4tMb76buKip1b+ zs4X{8vg?>|Qk8`T%BTcj>Ve5&W}WD?VOM6hd|<+YDqkHfvA0kX{Xe`U>Y=o%>C`(0 zU&gXP86k%jU1R#b%pPI#o^eMGi8*&xKwW;7f(94KR+Xaya9X=C$%9)(uO_q?CIoGT zhh*iL9;|{%90Rreloi0_W~fxVCoaoWKh-PPL-a7_?tENfFM;)ZE;`r(8V*3thCJ@T1p)QnvGewr5*!d zu!r-WURUcyBTSL`Ont*jMdtdZi$FAwLZVG1%u@A`<%w`>opvx-MvguVjqdzNQjaI^mHmI`G zMb)X01kD5`-*Aqzi&%4s1Nz(8YO!Otk_cB(32k4YEKaF6MDUZ6d8B-WaLIeP4I%>F z+Xr$SV^49~mf#Up1HE}!pp`e8B)@klSrm^{wtvS+Hi&r;BPDOw1~RiR&u3WhVe6pZ zM2W|TeLDr4+3NMFwy>)AuEDJH+;X8KNXB>js7&d6efQ6m`8XG9PBx{b4|N4HU*HyXJL;D5_I}1PYDwmHHDS1|5WOm<E0m8ai6ywMI| zi1Ax?U}{QmFydg(ZXK$|dL_aT4J~3Bc047$fnGtw%bz@;gp9y=b1J?dKTm3T%J4Nc zom%0~l?vCB1vvzBHA(_6%$C*Xp%oo3f#g~w^5o=ioky6P=XFqC!j?gu&7>c#amXj7PBV}kw-X%e!^*95pDYYEWb&2+gK`aX|C8q! zsYXB-d{^o~7GK}}aAYv+_418;t*o?lQkO0)wzItRp*|J$s+z}y*HUp>;H{qAtgC(Z zkW;CXANMzuU8&RG>nM_#-&)v_6=J#m4l(E?GWiJ^OYB%y6Gwy7?)tmoQR=E}iO}ATs z;eTH@pIPX}^ec*~Jky2H#PF&nlvsU&rov~g?_Tc+!xy26KvXCjtZFay2v04z9L-gq z7P%N$8Z>#95s*8ES0rVaojUHJ)+huqM*^kPo@KP@GMBKqxW<;;bnVD*BVJ>KWZ3o{ z>x!-PT!xL%S-rqz=&H^P-rnL81$>h*1@|8-KU3kr`K-0y1CJpkA@bKp6sNKE!#-vW%l2UBft*^_wPIllx?QV6YP3`P($EeT6 z$>nyV*-Yl6P*)*lR$scx5}$9EZU)hwy)uNHo!Z61X_I)t!%50Rp|->JTRhJkVmn%I zi!YnRcKT8ka&E6gSbV$hA|~pe3S%g53b)@yOwDOn_!{O@rYl!9wz5`gI$VuD73f+fluPk_x9EJhZ8vv2`YRN?-!&QYA z6LJ0q`|}%4NkTb_h_`#Z?U0YqZfbBOo-t+6b@7g$xfNwy8qek`HHc0vtZDB~;^?Mt z1l}#*@4S(!PH+S~XgRu?T_Jch^WCp*r179MgQcW0OS_Q*SQ2_Pwt1@kmajI69=an% z6T{e*R$+@cY%e<>t2n{RER*LBd5Dgok|M=HCSoLywrgjBRmF)WuBrjmFp4{ZRX;xJ zCp9V3$K(2;WwW<}fs!8(MUTIZD?%zdNM`u?t~j4~6P28$HMK7;w7vGZu*-u1lQda{ zHR^R{cg3=b4X2-v_XT+Qu}0J`e5bckKOF;W{tR{8#iOd$bRwQRexP;PT%w8+fu1I- zVJnm?aJiBnzt~J`!e)Q24lL5YBFkh4$8*#G7_9Q4m%t;}J;;LHd zZ6ujLCsYtdU1S8^oLsvKQIu&Kp`+taQ6k;QU?G@VAfzs>ABEm70RC$CZ6*!Da)mf-! zKVlx>vhpJ&ssc-j@w%GcBFNI25XwMi&h}diUfD!~BryGJ2-nqp)nz`C;FffGRG#q# zf^T}5xkr933C%bY$3B&EwLrQ@>s7Q|)B|T;P9ib3+4;`)zL^4K;}8Z?N|N%(nOd%= zzIjj&0;%K}YoYOjNUD3@5;bnv>mT2w2xY6)a;q_mln_iJV?U^4(N-T(HTppv6NMOC zx|Vv&vo=>E>X^t7FCH(7I;Ku~)0o{(LeB$Kln_Enmpt}lu4s8+XW50t5%<*>YV#$f z^M-f8?{+Pvj@E{f?uvlmk~K204}muekJ_!R^suRTGdmF7?beh{b^mIUoxs^E!pv*) zGpJZ=IhEpZj73JL?N7Yz5~Ih!8_X>P>?Z!segqaHo^qS=DA(FOIJfgUI^l>%!x+qx z-t65I&7Efy)%W3Ab}A#_@g5%AnUkn}dD9~DwhtglGuD}_HB1s)(V6xzNK$*5PEq@K zpP)=b`cyQ}C-kxdh(S`h)g4InzT5CZst7>x6C4|}%0qR=gkI3rW8HJ9l-nG)RKV0d z_R1{Z?2#q-sEKn!pvkdQX$_4epfncRSlfA9J7{HBoKv-hy#Ie^*Ro^DaUJ*d6}=n* zo?*TpOG`3r7+Ux-kYE^I5YkW>V}>+H8j}3=t`iX_t8VqZ-7S#|vI~`Uv$Co(&Ld77 z1cZ9^70f6N+1IPDWj4Rko;&(dEx*$2+CIE~DQV%m`_%Wr(1N&y0M#LIVs8c66DtE( zqe6ql!*MmJV5!2ig@KgF+|xE*C}cdFkf%tM*NPnI5!%ZiC9h$!jg zxgA!slN>HkDM1}&ld}^uYofsoXQwA5d~rZID-fl~@KZW06RxnqCWlxA$O3yUe_^$M z6seC76gqIXt2vSC3H89x@p3x-r z9u96Do6ynK!kIdOWV&RJvYgf+I-u*=z#f3_R9QqNbWA)nL`!X~=t)PQwz2@tvQ))m zuvz6y`;NZ}9mxTli_IFzOjGvb%VPD|5xE(bCCx0(g*Mt3g#l<(EvW)l^k!GxmRd5B zRMoK%LC#E;vmwT0&sAiprYrrBY2I9Ye{v#7Qc00#a^K2`R(r8Wq?6TFRV4vpHANjN zo+?MdWFev#g`g?>2ItoH!KDWaUTP_xzcC-QGxT>*fxtwvU|<%8QSqALz$N7~S>)tJ z8tgq_JBG^)6b2QiP!=t&m_7vOAa|=#ir6{QZm>|A3Zd2%`x7ZoU2}Lt5u_WnIqpJZH?cw^MLuY7)u;XHlhm62ay61|8@8F#85k=arCekb}p zH7g?NiQtU>Iqcy_iBGU;COLXC{_)xl6?Bg##B{|8dg_QfNMAy(M(GkraxpkyCRkTt zkuQrkId7`G^jz}MpD*LDvHns(DS*(j2||e?L4*a{L+C+8rQl&JN0#7q$S@>7i??;liQd_@IQUe0Tto)Bth&U+y zMXZIbmb`R~ii8UBbr0imD-Z62+xiXSk4noy7AjR~V?rwvQMr0)321zI5n7`plI+`C zVci|pz?#FKssmaaI%o|~vLpH`Gzm_RXzqay*3nj$OHd%l=@A6$B=04W*OlIVr(XJHXiINZ<5w>~&UsK~hzVcp6rEzDZl3Ybm><%{3MA`Km=GovQxu!G8y zz+fN!&kZG;MX*CK7?Pj|?6~4Op4~y2Gx?uF7C9?wSWwv*-t$;+Vr$8|Qnn*1lGvSn z)F(-))w?mFppxYCwXb(-WD$$J#1ad8^c?Q&U;!{9dcp?T2u4dR&+fhz$E&Wevq4dH zu^*I=H*YD04#C}#XcNNXbvGy)UQL;9p{P+DpHNhUy0GN>A=1-Gm06loXm4pyZAerE zy5N2F;%v&Qn71L3Y?3M)(cCE-X|*u+dbCzF&)8o2K3023r6jtvGPi?UJAx~W7Rrm~ z;UxxzeO&VJ9zf8DH9Q?<$u$}3+BRZuiL_jPMVd3x>k?P64pGXp=wJ5Yc1iWI zD#TvyCPs$*x#}9VPF_IFA?3wJBlD^|vzRIhV%es{%`^prC6UDIxcpI7G&*~F+IZ*1U03$) z=)Qw177U3bd4|%rM_Q{Uq3jgzTojvJB{vj2&YX+rOero*7i-vDF=b^zID#j}31x{i zBo|BuN)7f0l(3xF+QRQDPB0?lqPw+_SV@+=+yxp-krWJz`fNm0SUfqb;BOMxMDaZk z>CK{G$j4OFDJKbtfUXi)R^<#`i8kaSl~x3KyGm}Ghiy>=YxT8_!ab4QPR#b{j*dWp z5>(Te18Lx1jzWzZkiOIMKrL8P@(S3Bcw?OJ?%~-e5OX?FRMnN!^r1e`_Rt2bUu+KZL7CW z9~6MfZ80gaomaq>$U~iE8-uIh0_wQWw(7|{<%2fBm0*wZItEwrrHmaBZ{bP;{e7u3 zG^TJi;FTm|8*Ps6fe)lSaiP_)eIenJk8epcGP%RxDNnpW?C~umw1{pW4n=AiVtD0Z zuZ{N#X?XHnER<^>vU&H}lb}emDH={Fgg-pW3cM@E-7KHskyubLSHYLJ z#WKW_NJ>%B&EN;U(E$a&X~qvP#4JHhp&Eo4@L%yDDEOL`9n5NPji!Put8QIFd)Wnp zz?C6OUC|;_A>TB9tHUy;$`#}+~@$TlSlui|1Qkuy{8Bq#2 zsf|Qia&wQAEu!2bC2`%g(~u!nd)r)vj)*ZCk)0k+P2`?Hu7yrljYmHy>=kw-q;qD= z;8LAyHN(R0BAdG`^is+Bm`>;@T@9C(7gk)*)m1LCq<#^ujDo&Gt%oc_9>JzeAB_G% z5g4R(6D5khfeKysTs@7YV#qufQ(eDq02)}}bOO*k+y&blL7V$v?x(_a%bx!sz_<^# zs77L{tFaw2s%rQoq{%L!?vg7AJr4O{$rGwbds+-I5go%~MSTwl=-U69Ee2 zEjD2>Vy3lt+B$f2kMXIE(2G21?Vfg=>Md7SD~>t1-`6S4*V)B6SA>#sr-<&wbr}{B z+jb$LGc-e*WL3z)m8@~|8EzoR1ycP~gKV@qTj$oCXJO((wsEvYf*papD((Cd*DvGmO1eRO-K;CTwXsHCId+KZTDuehzf|KQ5Df&F-?uip% zFUsTkDuIEhE+-4cE9N$PmqscZ*w#C?0BUY{Y|Y zhTn5jo9AsoG=hKGXD*(!OZwQUHlF{_vjs92GqM7d^_ z@9?3_2QXISC=o&pj5G6Aw64J?W8zA>IE+H6H%AeuVinUog?_2t5PZ!_W!qlB5v%ph z$H*+Nb#x|81SzTImS&PFNV|ty+9%avRZ(ckZC8iwo>EN|hn0ns<$PrC)CA#JmC}ak z(O?lM=nsHZn2AMVjArTNmgs0;kzYcQSZ?6z79x*wwTVr@;F)`os9;80e^S`e7P7jJ z8q8<_@0_u7vo(U!N^&>ZunzIP;sv2#>aD1EVsRl-=XG&);f~a@iZscswO7XDxnQ_g z@4X2_3}gC->d{>2O;ztdfOtb?o{8TrdP@^yvL{XvmTBakqYMQ^KVXc!O)3olQN%oI zWT-T?nggNOUhLMEdOH+GHUdg_q(rfM8y;c^UPCpVx^5JrfsI09T$80VJ-MCbY-)3Y zT;xYySudU{nETTy_T_ELNJK12_)t3O1ARjM* z3uxMdyLxo_?dp>;(?ZwGaKX&0^Ce$T0~D2bbo7+nfGGO4##iHzSSRr2ITk9iW-HG{ zb7^clR3Jeyyw9b5`NjC`q7RXjTS9U}qCsJ|1jX&|W0yIZvCAbg>@vlPNHGb4&|n2G zd$aRM_#ncd9IK9q8KXCP!TV5yi|CBF-qmS0xjWcdu$I-df(E9T$UB5yK5|7t+XPV39&3wwK8P~ok+mc% zhd3#A1K+&JUFD$Xu@Ar?J&zq@kOQu)W6IStpUAMIDen+edP<2Y9N3K(nr-sOD70tW zpITbd&4h|aUv_9@`xFqy#he`sE%ikd3!VK!Y{QU1Uyy!GZ|=PZiiE?vkuYt-XV<5T z;089!F-L~RQJ}JGJ{mJlHTfG$+g}QLaAYLPE`2Ub79|y@sH_C?OIi4dDC=1L$)|MN z#^J$6TD++6Bm)+u1E>LYuuK9Or_tg)n03){3rP&KEg^N)%w&}MdJ8ZZ7l+*WU?$%% zC9fu6_OK07>Wb$E>+@!q;$@wbV=WfI%K}Eb%@S6x7&9UEFQr)Os0bH_T^E@>mqW0J zR}6I8<065fXnSPyldoDj9-v<{Iyth@37(99l6Fom;AYRMFmrS$4$X{>RuKZA6 zSmQ_9>PzSsQfMbf^#`_HWuoTUVsiX_!DaUFPjfkl;>?9cL19Sikv-F*m>13a1Q}wo zEAH@Zd;QmNUE zMKo_GdhWARHIEqr$=>BM^r_kbI3!k%Tb;8D_S#Woi_}t~_oIr(0t$0Yr_e5zqOOc# z@mg8r2R5NBL8W@FgIaiG-gY5j8-+~Col-1aTYB&>{H-c5MiefSk_+_a72CI=9O=2A#%6-F< z?nAiIo1SC$sC40`SBDllaC&cSpJH&&3*5Zr0-U<&dP7A(s%)>pW+F9XN$^FePxriJ z*u{J>=OuV4C_-+%2((3y1Npsr8Wh1jgW5KnbqgD+ul74on!s+SBDjg>4q$H#p&v4V zeafq6S7|eOO|)FJAyy{m)tOStgti^06;s8Ej_h2_fiA)ueZ%g)q3s`I@nWS-Z%3y4 zflaYTxyfC*sUhtSq8(J(hO@|5Y37(-|G`$(a1{xEqMuFjnGwU37QtNc+onRT>5c2| z1K^ZEpi)R*F@RDORh=QsG$19f`ld3a8R+c@H_zo!LYZJT??T#XtHy}EFX@*)VNUWl^IfW zcoH1vm0`kWAnkjw4H@?H%Aj-vRL9XJDmWUmLKx~)o6wjviQ(BR3)WI*NkG?C#AnGX zcR>LxNz6zF`I!e=<@U|vbwOo7*iAxToDYI(BC{fh9n(lgMNp3S-jdv`>^rOn#m#|J zG*v`T0X0nmxF!p9~hjpkOp!G#L9^nRf~sa(3^?qgE-Ul}VkC`;!Jc~Oj##XTGX zG0>U~Q?gKu8hh!3%0r53Ei7Ex#F?eps*ox{3bDbmBI`x^kEKz=@XNvIF*mR~NO^w_ zkAiw|aPeNU;x09=m_c^WEO{g{D=JpK9}`4_+dTVxnK4G$h&fSs$8c3XSXSj^9jkyP zNc2_}eE_YXRF*=XtZm;2nJd03%#=Bts64O)m9ncQ8vi1$YESVmB(;IaxqWE{WD4a{ zmK)ApSouoi);tbsQF;IeRTH3%ATDtCDmTF2U^TETT}y0F@-JkJM}2^xLz0wawV$&s zvzSLvVc45U8n!2bZqU0M43R_>V9=0SLxc&|hwliM5kZ`DdJMxPQcDyrWdSK3#t6=s z3$nOIaGNa9tOAik&ZcS>#!?@&iWZM*G7Du;iU5JSp%F<*B23#*DYGn{R})?LcCMQ8pX z^f1Fp8nV`XGoX?G+nH_7uDu9#Szun4>6lc52NbN{ECklu1hv>(f(HQNkAHmn<;QKG zf&KA*|8KPfaOid4XJUgC>D=GYk6*WMzIz|u8;*YbZhQAW{o><)Z}uYW|Lpdpk8W1` z+iTtDcsu&t2Mn-)UkJPTwFaQZIPPEn?fL)tuKwcvzOMCM_{D#H{O9uwdxm!%b~L_@ z`4IL&W5jp+)`DRi@&5d~#h6U7Tih&yHQx8t=KaNWf4{QX?t_tQnfo`Q{~HNcusQd+ zZ{QJf#rt-NOT7EFhFKNtKELe$w|`1D%=AU~?{q#{HQ{bf?O&0N`pl(SH(1QOP0c>d zXfeC}73P^Yk09KC0%bta$DESe!TtnJXF^saUe5~;u3Q6xVn&<0j(Fh)EH@dVt-x{S zi<7<`@OaXf3%_B~weyeGx6S$X{D)J%TsFQs<9qYCJ7IGx;_d^!b&+q+ ze>maGWnZ80?SShFH``c83OiJQ3Ed}8++XXoUh98h-=A@3f7Y&x%SGH{Ly(|CW(P9`eqO@>* zDBoW%+>trk-_py8f3y?-;iI9}exYtJz3tt@?2Ui>$ERmrUd^;}Vq}+}uTM2I?!Uso zwSdVzJRDywr5%PH0N7TB7q+J{B&AgcF_n#3zU*j+F3lEb21lYL>bQ@tR6#3LcN~jE z*NnEbKg=$oZktS#v>|7~UNiG~U~Sru^g`9n<3JCT6FSCBhTHt2LnJr;tQ{!K7p{qo z1`P(z+~zeDb@f{tlXZCO-e+J+!dGV+jlUeLH?5%Al6h-=XJS7>e*`i+_AM=%QDzqI zK6HPMcOPKujEKfeddSzOrv5N7yU=5Uy^8~Iy;^5q%V@1ZKvr%846qG{m@8`LZtu@M zW;easBn1_yhUYq2n-^E7-)@Cq`xddq;=#6V(eZ3p80A@Y69^4FLb(?Ew9Q_#BUq)L zSz^Wm8y+uvbWBk!dHZ((AsMDWXchNg*sg4`%g9_)zVsz9|r?q06 zy-dhMuSM;Xx$9%u`=ofkKbjZp$+Z!j%t};hx7`DpZ!{~{pC{|`_9FMK%VXR!GH-`H zG&0x7ldX9j`Zik=e6+NCCD16|J)GUUdTeZ(sknP7x5lP@+Z7TMa4CmxiT`o4cXG=6 z$`tdbf>lCQRJyj(pRe=dWyTejp>TT*H$m<#vyHQWSl@3TU&noTFjY(Eh24UA@I#L# zc(OPyws@ZLOKyagt?wSUn26Z^wk z<{6=OpDjT2$~-S;zQJ<79Z>&9EJo}2amHfwphXl^U~CN+7 z1ntw-hU(xL>K!BtcaQ~x(dshP%UG@RBLm~K=Uv9VzW;&^UPle*UG=>Gt{L5)cy~rO z+wWySdMgMcR3BQ&_0Eh=1J=ZD&wpX5ZZY^WuDrffXvJRgd~Ky}&U?c`-JIY4Mwo&$ z+w3QpLLKTn!4z7pgVg3Q1xiP|fhh<+-oRAmcDzsu4JwI1>F5APUiyM5us>mCJ4~VM z`v#_{DH%^Ng|cttAQw!X4hg1C!vs^ur5&b@pIzm$U={ zWe^$zdMuput4NY@|=@1$pdF+@-#@!vi7zi*vnch50RdJ{Pgs@r|&oW z3lG-f#MlVN#A*98{^r~8^v9o`zWMVGFi*e!VRJKJe;-hD_owGnfM>t)%y?Ko;=4_3 zmPIRL|1yf1lzo5Gp9eAj<;PE#O)XlQX}7_Q9e6WN;lQ=89s{h$Y#XWFA8@*~Z349z zjqGyf|NeN|iudQ>Bdqyn$YFmN*AbuYPsVz#JerJ2&gK0PJe$mS*U6wh_TII>dvIwB!t^oHvjLwU%--aBVzICJGT`^u zxtOx#-G%x5>A|(3^mT9c*QhDF8}HK&&z|ScfA{p@POXk{2b6pr^4zdvSR!fcvg{Ie zixzN0rk0b92!9I6U+*UxO*tV~R;{XUx5VRkV9-7P?vUcH`wMtG)~wk3zXmc&wuJCy zM0)TGRyWpb3D@b}1QpIpXikU4*PqSlYSCEiC$!pIeKz|CkKTBDA@){Yy|f9)_uSjF zNj#dhRpgaw%#!<;zRm?)&{mAT^vYaNm5|%+;Jo!AWi|*6{W_puFKy)@kvWokw^?;& z#V>>j!u8lKk^KQfTXbK9_-vd0?GHcx{OQk6zhQc#-de23rOfQwFJ+48=S#i&`N2P1F&p}Mbe+FAkLubm-o1B zEhvytC?&evIew(N2-00rTsH_@gYehO1BMvB$DnoO<$*P_@pQBzl7wJT%$)^T6y5v& z?Zob065HJvFL=k*nTS~wHR6tO{R_qo7u`s|OM6nUOFtM;(LF^9rKhK;!hrwOG zf8*V2e|zo2!ntQ>9=+#$=A1KUW@4UAZxlLwP|W5lzH8?0mjpQbJ*vH_c)U^euo0t- zvhTcm@=`P~U`yv6JLXvzarPMVrb>NoP2cebxC-O&16@kYBpnmwMIcDHTQAy#c)-Z7KrWyim7@uPk5{g*Cpn6~(JSn%-F zezg_pLp}^L>A$tl<(SJmKYyC#xy|oD#laP8zDa7bB(&C(@K0a9SywN4c}iBz_Inzf z3;edKQKtWfL*g5AI`vtW_x{Md)p@U`TX!5(Zkf7t@8z4^Iu0CvD9|PNef^v9qh56H?{o;p6D=Z()Z%&(^encll;?OD0aSl=i1151qy z-F5hRZp`Y2=BA6TWYuh!Y80|-f92@epS;gzXSPUaIIUG>%cB`4TRuIGxNd&n*=`dv z@rH;Wi4hsF2Jv=fcX=?bvvJFD+ z3o9NyVLnH>Hq7$L@NF3(K0Y@x<_$bF{6J1EL5)!dtHBc8DkdKrcz?J}z05DxhmZ8J zUis}!`On!;YUKuf8}FeC>Q~Oo&vR7ZNH4eH_?y4`kaBW?NPwRLfj}ykOQf`1swp=! zL&>CTNT7AkK(|0I0V%*THxc|HFzM1+AQqT-xWUg39dKbz6H`qj6iYan`-QiJB!3Bo zW&s6%4bpsFrqN{@ovg*AyIWwO*GNBZxd5d=2`M`ZOafE_mF73KK&|;rEE2$*=wh(| z-kRSf0!UBuyHp^?e>Sc|EEB+6^SfLi#{VX+Tpjth-z5SG-L6C`fVZYxCXmqWO5_4~YswV@3Ei$lDS)@8TqTgu?Ml=Fc$0Fe zNC0om?_zH73l8@fs}4nsuaLmvy4h0rQ4ON1@I>2 zGLZn@n%~6&8Qrc-B7nE1Tq=;!?aE{Vcx%e#0vX+|Od)``rd%nI(e27q0(fi6)dCsa zu3RL5Hz}8k1@PAVE)mG-cI8q5yfx)Aft+qvE*HRCQ?3xm>2~Ex0lYQkDuJACSFRSo zo0Ka=0(fhF7Yh`0y9$W_-kNf$KtZ>wkO|UA000Z%w&Upr+eZs|4`Y zlmj;CcEutD2i>w*1n7Vu%ZWt*4+#1f;6p^WFBSnnAgBfi5z!pPB7g`4)c_+Rnuk~f zAc3G7phU!SL5T4vh8SbH0L8@lBm<6#^I3*4%cmKzOq|a%pqV(IXuvaZKGP6t`BVd@ zi5b_#Ks5|ap1?IRW1ASrhCwyRh4D=cbi<$;PfqfVn4d^Fk1QY}R zFsKH(Fbax+ff!VSTo?((z(EYEK`xAjVqhT#)gTu}L^1FXL!$u^#f*w#U?K+9AQwhP zF>n!sYLE+~qZrtTK{d#Q5mF3%#Go4F!YC;QMq+3*prn|QQVg8Lpc>@DXekC(Vo(io zVZ;;zFEOYFxiD&qfteUogIpLn#lTGrjRxctGkS`FofuSuTo^&cz)uXSK`xA@DxGMqXVo(ioVceAfb1|p}xiIcZfVmhN4VWuo z+?4=xF{p-3T)MwVfVmh{gIpMQCBR$^szEM{yAoh72Gt-J#$5?87ek`~b0v(s5@0R{ z)gTwfT?sH3gKCfqhE`%*CJ@@2DvcqN`Sc-RD)a?cO}4F42=fNl`!r~fVmh{gIpMQCBR$^szEM{ zyAoh72Gt-J#$5?87lUe$3*)W?n2Vv&fVmRJT?sH3gKCfqhE`%*CJ@ zl=wG#W5h!ni8|=3-C{a$($+0CO>@2DvcqN`Sc-RD)a?cO}4F z45~pcjJpzGE`~+}=1LfMCBR$^szEM{yAoh72Gt-J#$5?87lUe$3*)W?n2SL*$c1rN z3e3eIniQCeLH~lzBjc_Vn2SL*=s+^=N`bi;RD(_=qO4VWus z+?4`zF{lQ)Fz!l$xfoP~To`wyz+4QfK`xBDQeZ9y)gTwfT`4dZL!$w6rHs2$U@ivL zuv1Mvk^*xvs0O((?n;5V7*vB?7Tnwr~E{wZUU@nG61LjH@ccs8w45~pcjJr}` zE(X;g7sg#FFc*VrkPG9k6qt)aHOPf=R|?F<&}hJ1DdVmbn2SL*$c1rN3e3fz8sx&b zD+T6aPz`cn+?4`zF{lQ)Fz!l$xfmJ^m@8%6l>&1ys0O((?n;5V7*vB?7@DxGM$bVo(ioVceAhb1|p}xiIcZ zfw>q|gIpMQrNCSajRwq>GVV%&xfoP~To`wyz+4QfK`xBDQeZ9y)gTwfT`4dZgKCfq zq|gIpMQrNCSaszEM{yHa2-2GL}|Tnzdb z44pCV%7D2TRDhE~%*CJ@u`xG#W5h#<(j3=3-C{ za$($+0dp~^2Dvcq%7D2TRD)a?cV)m_45~pcjJq;mE`~+}=E@j%Wx!kvs^JhP-CtzD zTnwr~E{wY}U@ivLAQ#46888=vYLE-#t_+xqq0xZ3GR9pQFc*VrkPG9k448{SHOPf= zR|d?*pc>@DxGMwZVo(ioVceAgb1^g;FjvO7D+A_YPz`cn+?4@yF{lQ)Fz(8LxfoP~ zTo`v{z+4QfK`xBDGGH!-Mg!)`7@DxGMwZVo(ioVceAgb1|p}xiId^ zfVmhN4VWuq+?4@yF{lQ)Fz(8LxfoP~To`v{z+4QfK`xBDGGHzS)gTwfT^TSJL!$w6 zWsJKrU@ivLAQ#46888=vYLE-#t_+xqK{d#QaaRV+#h@DG!ni93=3)>{4$Q@%f5FKE z#$7ou7lUeW0)cT?4$Q@%8k|I6+?4}!F{lP75*T;oz+4QC2F#T+?#h9=7*vB?7iXbFc*Vrcv6LW1m_IlFO~~1 zSI)RA2j*hX??En%yK-PI2Gt-J#$7ou7ek`~bLEV?a$qh7)gTwfT{$opgKCfqhF7%*CJ@q|gIpMQ<-lAFszEM{yK-PI2Gt-J#$7ou7lUe$ z3*)XFn2Vv&fVpzUT{$opgKCfqhF7%*CJ@4RT@Jl>>7z zG#W5h&bTWF=3-C{a$($+19LH`2Dvcq%7M8URD)a?cjdrb45~pcjJtARE`~+}=E@m& z<-lAFszEM{yK-PI2Gt-J#$7ou7lUe$3*)XFn2SL*$c1rN4$Q^SXuw=KhF7%*CJ@4RT@Jl>>7zs0O((?#h9=7(`P5b1~>&aO#wCR{_k$pchE_ z%*CJ@@2DvcqDuB5dRD)a? zcNM@~45~pcjJpb8E(X;g7sg!$Fc(9k0dp0My9!_~2G#HsH{D+pz+4QfK`xBD3Sce< z)gTwfT?H@~gKCfqVo(ioVcb;!b1|p}xiIc3fVmh{gIpMQ6~J5!szEM{ zy9!_~hDHPCDj0Vaz+4QfK`xBD3Sce<)gTwfT?H@~gKCfqhE_%*D`X zz+466t^$~gK{d#QaaRG%#h@DG!nmsd=3-C{a$($60CO>@2DvcqDuB5d8V#7MVBA#z zb1|p}xiIc3fVmh{gIpMQ6~J5!szEM{y9!_~2Gt-J#$5$47ek`~a}|uc3Sce<)gTwf zT?H@~gKCfqhE_%*CJ@q|gIpMQmB3sKjRwqBGVUsYxfoP~To`wiz+4QfK`xBDN?q|gIpMQmB3sKszEM{yGmd#hDHPCDj9c` zz+4QfK`xBDN?hF2%*D`Xz+5Hct`eAwK{d#Q zaaRe<#h@DG!nmsh=3-C{a$($60&_8_2DvcqDuKBe8V#7MWZYE(b1|p}xiIc3fw>q| zgIpMQmB3sKszEM{yGmd#2Gt-J#$6>a7ek`~bCry{N?hF2%*CJ@hE}%*CJ@#Mg!)m7U@ivLAQ#466)+crYLE-#t_qloK{d#QaaRS*#n5QL zTovQ43Yd#QHOPf=R|U+)pc>@DxT^x@Vo(ioVcb;#b1|p}xiIdkfVmhN4VbH9+*JW{ zF{lQ)Fz%{=xfoP~To`v%z+4QfK`xBDDqt=K)gTwfT@^4F12xTz^nab)hj@7eax>CH zR6tq`o^Ivuha#ONKe(N!9kdN`iToll6sZWn^9v~8Y;U6OC3txnJwea7GLt6bgu?u7fi;og; z|5sd*^+^Aq0AAK)8(b;sxQ(8$<&kcFqXKYS9-*|@RuDKc$cz5f1(ve(3ikEza<;a> zLd}1gv)3s9pphP4qXgXl74KiB;*L)A-`3N|&C}PlUOE?kh|f9^miC{l6&qgdPUww*uPjjvoZG%q)F zUbzN&Z$^r0Y{@%JO0+Kz!DAc!|EK#FchLp+U!=JJgA>SqucKSwNMCNx5^`+`k)WLj ze}h8`o_-+zDH-qU} zSCg|NcHR)}8`gE%!=XuTqbJL6jTkif=HRpyFI5ShbCP;*GToZ1A|{iHPD->5$0jQM z|JQ@z7klS5K zExbMK>({;0Dl~h!D{NfkscO4X^Ae?8z-1-N`aZ+&lC^06v`s4}vOSmEbJ3&{Zi`lZ z-!!iGfZD5#dv5N~Hr~i-VXF*c^6<&^zvL|@R+AEK!?~uv9v$cULj|fhpD6}&TrPxZ zYdHGLKcGV$2h?Lf4eL0}yCf%df31!u&8ifDMxE`K{Gm>FkE&0Eea1-E9B+oU+P`U~ zLR+VQtTzU2J@T!>eq0u$T+_+%XnEY|(gW9k=1SM|^Jjck&W(1?Wm|5uJKMQWOw}HN z-W|$J+B6_+QE2o1oqNW2AUG#=+;397p?RxGi4M#0AEZaixdw6p%lT9>nB_7d%qPRk zIbaewJEZ@WuAXIi)PLMZn=wQ0+>sv~Yiy>m zSzN*=?;A0jl<2TrSGV9lY`SqH8+^;TQ6IXG8O(C1#3E>_gXO#fIcuQ9MIA!ua+7hO zU#sNx9;e4$Fa=r`mkx_e{0d=xH=T3LrB zwjjk;UEX>TNG^0eKY!Yf+kM9LRLG~X)e`rPUijvfg{bMI1r5(Rn3Q`HZ!^!kIjP=+ z=-u1emH%O0Sj&wYsG$b@_eXLvs#m~qK2;3nxNffhjpYwF#PvZ7Baa+^a{k%(0rvMF zj#_>3^8vKhfA5F-Xs6w{^*(5CvkpCyaaq(p&Ys%E9+K_F+!ax~bvfaQFu7|Qbe znk`cZVdgzfwd4-pa3v3nBEi*JLb$?J#0|{;S|vvZ@rADK%&D4;E)LCX=8LX2KA~z@ zqk6LwHK$yt-nymzyz&-4g!Uwch@o5J~X75;RaVmpDn>x zO6)Rz-_CV!>hvl~)W*jTI*(R&yq4gN)|0hp|LoRgRD$EEx#OaGJ#}00wEma5#f@vP zJNfXUbJ|BR2Q#r>>IgM$9H9IxOe+koCZItx-4bJdB|&mx_dHxW*gKguvB- ztmzJ8ZE(s8=NQ2inL@Y*Rm2T2|5`0SxwgbWH0ox}Vp&zqYwgLIRW#6ejLAr}Rgq%d z26utb^HXpz#%2H0!ydRSY8Ll;okq*29@@49EvM`G`SW0mX?n{Q6C;y1jWX>aFK#mS z$m8N+udcpWCkagK@F_bujBB~eczyk#;u=z%xPw7T(!oiZ<%(%wrK4M|3rVZbIJ$7d$iYL%=|I9 zENc4?HQaz&=mc|n)aZJC{_GSLU-jzuLEY`X6q9cId}!U|Tl0ra+*Iyc|BAuFLzkz9 zX>9IYq7}5)Ir~Y84$JwSWIb4}HSorrQ8l#XFuuy#>k=WHiogpcc*AXQZM6`tXd$P` zeyx_Hh4=#3VvqNkHV9mFnZ45}4Q=gNX;G8*Cto#J4on#qwBU|oz4)8E68fYEqIUYJ zbR91Hc+_%kK!&VE`{zN6uWkn#Jv}#W>ieT(I-Xy(e%E5D_vwY97jN3t+40!cthx8p z10UT;i;W!7Nr?{2`JH4vSgtkn#+^1dwB?1h-ZTTOf8aQ}hc9$b)w^y2IvE+Y?%8kszhRrd6eF0!>TV#?140vbVux zV`4Kod=s-tiH=&%?2^ZjzIsE$8*cBxQX>Vc%?hFZVtMJfnNR5@V`ii5Rc9u&&wr$zD z7Dlr-&U{^mV4d{SNQn-^u?t22|Nnc#$s9KRFbp@3jc(9TbzGtl!sVOL6$h;6S9cg; zg~{qT!x@j{aF+`*b^X^YKRw(d8!fzLlVe@Q<;lMcFQ`p>7x@-7x}KjuQ?E5|(;{zIv7Mcs4ROjFR$^TLKH@Of^8c*& z8d$MI+|(^$Ttl7I@gU2}*Yj?Y5*?Q7cjS#sR4iaQA1a2joHzI)haMW7F@-r}e^bZN zKzyNT_U}zjpox*s2F?Rdb06NC_s^@`&sjO<$&c>-cysk*bmoZj-?I**Hgj{-L|nG3 zOSew=U9uMKpH6w3yE#=~uy%D`_o(PrO)Zx6{5Utp_sK`Q*l*U(vUg)iy|!%J_vN>E zQoBis4#V|3@J8mM7BHL-6+;;=k;7;kOyGy}qrBevU)6Eh@JBtd&Ng|MGv{Tcm|Ims zKOVo;$aXK9cUt&z8QMAR)O#V0@nXkilRH;aPU_MlxbL;G$I=>a-b?Crc%5S3^S+bXO-giFuHS(-GU2#@<$S0Z z%5q*49p*-`E;&1KhDUO^F#x$+;jdvG?ZX$^XA*I|G}n~sesyg2T`(o>bka2w&kdgRc@-(qVYz+>-pH)-0+#c!Vlc}KbkN{_J|aAU znt!hJ4=hLf@P$^NyZ)>gTE8msLMhx0hRxhv*7;_-sZGBzi5pT@1RD#@|NXqkI|;RK zMz7AG)jy_hT& zPM>ilIxN@kz#F+fq=4motQg92IDo^=OyTkoxD_KCkjX!b0ynM{s^Ok1+?6Q5#&vhg zs&ml6rPq!68xQLE>CCWF@Lb*Tuh6yOg48GIW^|1!mgw&Ja#?F}S#(#B)^`eOohPM~ zL5;5G=g%>NyLOf(zDZjcG5UkA$HPgzW1b(FS=zd0|M%g|x)xtQi!{_pCHLIW`bb`_ z!baD*%Q|R>Ne?|Vt${c0Iv~11LpfeZmz*GwtQ1cQQ6a{?;XZ2oOm8JXSXePaoJX)E8M^vS&Q~htLO>Ru3a4?=9m{5 zJagvCD*oqtHC?e{-EHO5Ema&Hd-#!xt?7t!)mu~5wVrIiUAsda(1Yb#18>~LUUY+o zvRrqk{I6VEQ zQnN`Toz!nqqRn;g0xjx*9xVTP;7!e4(MLCEFw3RDbs#(R#Ch)Hm^$1$5GOw(*x~k( z++|F^2KGIQW;@Wr3+IlUtugMK{WEX3$F0A7T#D9q3VwtAueb?L_Ss;s>!B7KKH##b zneR`*nNM+JL;XghM%VN6XL7&q39XOS=@viv@;29*l2yTLx&}_D>$g5KTfV4M1!WJ? zJ0q5lJyj|??=vaUVL88-{GVG+4Cc>(QgfI27~FCkS%doz!_XVAhyGWVcilI3A6h6J z_oK!0x7QBU>JYx@^5gi4XzhXP(=EVFm+<+0+Jd{W9c|jx#AQ)Sd)~?%t)6_zc^q0# z*YoqI-OHnv8JQ<{cW+{JHg@NNkuyrJEIEGaC=c&Jvj(Le-jhgd9Sk zq*jV8XtAgBaBgN>&9N~hj5f{eTk%KjVs)IV?wm1&(}6v`dJUHC4(;R8-*rdberDTR3*(yhzm{@uA~}|mAvjW^ z!*YHvSr3-~JnIdo9p#Vc8w%{uOow}s;B6Q2yeaLL=btQv8^{acq^XD-JNz}vi?%+= zoj>VSf9hnR;<)FQiWX`8hR;5Q=FN=hDa2hN@$iF&xGP-0KcXVqjam_pWwlT%azZ3T zjjreCPj%cuZ_`d6%Z+TE91%S6MRm1xbEo-B?FTkceX_H$Zx&D5>ZB{)?W{w+AJq#x ziN;+tO*yLv%hBM1wVaO?Ls_o7kH$`q|AFJV=?&_lgRyt24ZhI4)-jh?Wp}J`Fqw_c z28^`siEi3`%v7Sg&BhL|$GMHJFHi1|TAtKrA8K?xKYt!~%A9;;Q*@O%mm*(GcQP|R zzN@*hpR>*W2We|pr?fKePh2J~HsUoY(NW3yon$>Yj{X*OoKF>lIWC7Q67VLoIKZyy zpXZ-A!yX#kXN%nI_}8kq$M?ljn-TIEFTH_k+`SKk!SGt~`KmT17AycgxbZKHt9sd!X(qe25%~H(y zu&aexo3z?iUvE|;tu{@g{duLUd3OsNU+3<9pvNkX_7*h!|H81{VVBf3bDuF?%wv1UcM)1zUTMT(=HJ)gvFLwj4jbLA#Kq3y5kIK4*g?*7&H zQKRem`E%oqDw~h5Q#}l?B@LfrIpXB5uW2bG93IX+Vt4)M?u{dk2XQu&_8Kvpl<27A zdL4OFle?ktM&8ui9Z~2$W~e&OJ9^7xLb%Tj-b9x-(G8yVB7`SOksD9{S}jMr@dd8M zo~kgdBe>|Yu%6Kbw6$-`g@??AvIMhZDQ(8|^BvUC{pOOCnbDO|JD+%PCocQ+`rSVG zU9uMKpQ)7wKRWngp+|Vb^wVKCj8ANQX4TO*(!Ap3z^09o&qPlvJ!k9Cm%GS;oD8{= z5^c-5yO>Z1^nmrBhu+lOZExrX4Q07R45!ZEu9?vN=2dv;i9??pZ%>acceu4GTzd+T z`ZdeD$b7j@`J^Uq8n;}vs7YFtkoi4Kn{z{Nq1X3#V-4Oof9!Iw*Y&cE`B}6ZjjcI2 zog0a5o4)caT2I&W^QTv0>ypZ;;Smw-vi|We_oIrFUB5=NGcVt^?dbm8YPj!9(l^ue z)540R=jDDOGYLv-Y6dVAjo4^_}bljMzl{H&)GZy)(98Aa9wj!qQi22 zA6XBUYYn{>xQg9SmKSmsjm=sA1F)lg_yYS}F{8HUwIu?IyK9o1oRTW{<0d)|0hp|LkRcVa>ua5V zZD$l;BD&<(*7)c|yJfy20&{XMjUZi9qBX$I-E@jpy8qsCa)&bhv7Eb08Qq|vEa%O8 zQo!^Gcw7qH3h(d6VW>%Hn`8UiP0_{`L3^#x)~5-PL*71?JHI+pY31Y7LpFYTJU#t- z{OYf$%@3-tKAo)X3oOzlC^06Oxf_kZck?FCx82_C8KtZe?7xu{8i7XL)IN#+c7K2VrTd8 z_+nM^eh`?GBRVP3VL88(tOv`rhTaNX#%`!u&YKe{wx`4E3=bu>a%P5 zhZ}g2YN4a8yZk)Qrre60HSE@@<^Fq2YEAJpzwJKVLUCBq@MU6`wRJl^&MBD@M8HnQ zTs0-yVCQ#}^MGCNrTsNUca1;$>g6p`-4i3?%@k7xB2@n6R^ULeKx%)R@XhZ zPEWj^v}f5x^IN?Vy+%o%MZau4vv%9FsP(Jd;5IHREK~6%ewVC8`=^cRB1zX&+mZ7} zH}{$qzVK0SmvxO-R*ybm+W2B+$C*wZnhrVFIVaVdl<27B`W<>BHw`UtDCbkfP?fyU zQ!1JZ75<=-Q z>1EIH%!%s~@g?p^KJfIu0iA-bZt8ez-17<9+qyM+*^Z8CSd{A%Z7>tRusY{5iq z#1(oZKh(l{cIx|9>n3!YS7-dQo?W+CpS&|9I>YhG`bT~{orW~`8(wdAUNjkRBZID_ zL`Nmp@7No;`)`3t&ZmmO9A~5Iywhojbu!iM4;<&!aaUPs*X~$veJaMS!8#sQBky%; zgX32Zc;tNuH=VQXPGa6y^WH8~P;2=_^)EVIeN_Fn=saD|&z~KWBcnPtcI`3V*F+s| z5x3r@#g3DUgX1=~u6E;*ovCY2;xajlLcAs=+KzKK?WM<6JsiljhTgclJkt#t%5k30 zaKBbHJah(bWR7h&&8&3|*2%+N|Dck0`WlOseDTB&oP`4pzsl{{t=;#IZCB@oNA>;O zXsJ`)fvg!^CBOV85wA;heg66SM%+zsH_-J;J-R11b1+ak?`%vvsH5A3kAaZvABaSG2mus2)Gidb*yUKRcDp zebBM&lM@rF{Wuz4Eq$(QSogb&2b$GyeQ>kCok?^>QnAUH8>!x;M2F@4*(v{hB`5de zFJL*JDu%LLciZi6DmirsU9HRJ=o)$cpNHjMZESKT zrFh1C)W#H>vlEv!FKU(fzZ`c|cb{nImv(!_$hE~6-=4g8X8Xu~tB=-_bZub2-RVL3 z<)mVhV>qeag`GF$?s7yas|Ux?-vR?~+;bP`1`Xvn?_{YQ##+@9e6#^?(mI}Jfu|&h zxQ8M9KG%=UcRvaCxy%wATY$Sl_r#Wc#D}H_$j=0rc}%Q$zwYDR4VQFY>Wtd%tDXv6 z_SR5IU2wk8=K8))?Tq&?Z{`y`p^sgt#O_k$mVuqTyIiPt=1TnP9c!PJZT9BI!@JWp z=SsQx>!d_SCFggM^_Q>gIlXS_2rbj#peDL#4+ZulmIJMeQ4 zi})?ax!VXNjvV8RF1kN@zX4sX{*T}1axdRqtuSwO_ZqjKZYy=-{AR(~n^RC5DmCY> zTHCoinmZI=#fO*HrAgbJHb%vp+&$WFwWE5#z|@%@&$~A4ojKD#@%_5FTXVW9dp-Dc z>vV6@B`0I-q(q0~{2sC%9M>9m;~t=5XvYf~YGs$7{Xr#1`|yQqR*g=>3s2$#UQJF# z|6;wnu1+5~ZQgOq8`W2&RuDXFGHixV!$k?Gb@fd-gO$5U^3BJ%n~=3=|FjHhn`u$= zU8#ty-aDtXds4x>b(uqv!)*FWmy{S}zstBC8Du5nFPiFI$dzc^eoP(EgXLOt-njd3 z&%6fSDcnRtB*lj>@#Z$7<#1Z$Fyad#;W-qR!x6+{(;Blfk>(=O^YLeIMOy$s3cM*Pol&xBVQM znSDDm;0S7qo_x55PTRg4gO{%qTE~0LG@9Ps%WT#9>UAovIUU-5W!UtQ8}=@E>v-_U z9_QOOUkJ*}Z`riu)lC9)Qli6g{f@hlCtVeA{QttZTOp$`Y{vTU*Kxk%YR_2i(u0ww zoXcTt-P-P3rSmcA{wEu58o1ha#ob4hddB;`^y!Y8?fb4X(Q)IiO^*F9$Ga+CHM2;5 z5Hxf69lNQ%N3({At7I)Jr#vmM9r3}TWWC3vUXz{~sokVR+i~szRTQ#%pt{z$8}~#r zx);kssf z!||mOjny?yndMe+uUdcRqrPF)tE0BCr+W%6d#Y7pU$h?f05*?P~KZrV@2g|j_-MGir z(G42P^8CRzxJNJCa2uX;#Zg^z4h5$^ir_JC++)~&*>N<6JCN`C-ewiLIQ;WJ+$?qD z_mgf^iS|yby}sy^re^0S8}*p*PIgmiiH?UXG`48K@1tG?u&7r~J7&T3_gq-0tS*Ha4<@UQ1G5_nPlJ?(~FXisGwgR@>XjePz*Z z$I9L&6`LHxN%bZrIxN@kxEp!)WC6?hR56t0ypt$!|1Y@Tn~b~tHk-1JmP=JI_a56h zzg5e*KKb$3?>o@GIX%Dl;BL_K$HxZJu{Uj%xiQrSf7-REUY>K=6{8m1Ms2@i`DeK7 zo#s>DX}FGif)?BNU4D1?fXn^9qc=UBHYWCP%JQ$FDc*UrOP&2;c5_GG@nUsNmZZE} zK@R0)EQXZmuw1`0-pB(+3s}y_ioq<;KRM3tps|zVe*kv04_|1XO}BX;(7ySr8fBw@ z*R$jA<=pY9`}l28rjDIUwbBeivT` zmlb@uHU)PVx}KjuyWF3!x3jB<$@C*d11H5DAK-cO-wsRTc2AKE-`dFDwVws4*<=8Q z)NjqWtJY)~?&+_zG4Tz zJJ)k=?5wV;! z+ewLzTF&nz>%nrZkvHy6T67;Xl;u2{Lxq>fBrp;~EyoEF@HZUG#R~4Sx!<=uE&b(e zv~lvv+uZG8y852k*nHQuBSNd`M;(X9Pp{Rq%Y#7U>Go*(x+5iS;j$OqCJp~zmZvPV z&fZ&V(o%opdDj9>O16u5ce;<~;dY&mhCXhaEy zm+P0!%-C;qUg=Y9XKY9*@0H07=7-DCao=*e>9DM8^3?Hz!SO=-XxrXJb=dyuK$E^t zRbl<_j_K+iT5m(|72fHIW^*oWO_x5IT{F*Wdppu#Cj)S#M2F-09eN`Vb}!&KpDKo` zfml( zpRkAfl^|Ij3khe=(Tj76MZqulax#h+C=jL~i zjG3|iNL0|py>=^FD4Vu_cDtElwUrA9%SpS95S^6hupF_aN6WRw-lTl27|ing5jVID z5iTSKhy&N@nKNv>!8PS@SBc;4m2-`^oyP00p?gDYtrKuhXf!-~Mn#vWMXSGi+qZt~ zx#HzQzqD(eyR#%}Bl4R5i_7lLTfl(<_X1t7_z~+|6E|*bFyNu(77s@UGoKN%V*Pq_ zKbfprQDl3|%^qVrmJZJQL^@~Wcuq=mIL_}R>!Fh4|5<_QFa<7RHwE7( zXO~uebpdz~XJ2+Sj>i<*Dev)Qm#F)?KAqO54|7{wF7(cglO^PyrW#|!D8{RP$0J$#{ic8#JP(Y?hT%~qgWcLroWebY3+ z=K10UuO3+4TloI8&8N#ly5B-=_M2JexU5ykv1t4*S&R11?rzP7JKhp)SaaC!=E~h} zc1v$$8?6^l{c1n#!?YgD4*wt|Cr5C?b5f$ias3XwakB&qI?jiRp&aLpsYAmZZtRB- zqtBn%D3OZs;$jK+gvQ_QpT!*wT0I)vbB)OE@Cvm-X57{j+`|_}oSc%3%c>q(^uh0vwP^pe z&vYtkbH;CBWO8(|1xs%mPdT`tZ?xdZQ;UH%+ng%9y(1(iM{vS(Qli6geg|0(j%y9O z$@ow)l;eerw6fXDe^AM(Nobqe&$l(&IP%sd7ex0~d(-esAz|B2N?OFmy#ChjQJrcT zqkTimqvo`7BR+p4u0{G~?!I7jJwJb1wQTjsVxRk*+1-1dni_np-tZ6KpDnAL@<8!5 zH_)L>-#4UUlfF2q-lRl_<@z0V;|wlX$@x?ZMj!zU!cfWPF++G$rPY>5=-8ABC z$USz54+t*^F|DfS+``wjNC^P-@Ur(;U8E3(3ekHi78hqIB&4(@nZ@( zk`ua<5*?QFJIH#dR*p03jIT9&aP>gPjU7;1cgHUSm%Zq(^yOKub)r0F z>5avE`!9F*Z?~g$;6BOm8Ko-cd6luOaCCB+@kvvfl$yFI?_%O3Qp-uBos?)>t`dPK z`WuAdF5N3=IUg$qvpj#I8^440H-70IfxUYlM)?{^Sj7; za9nHLt-uxRhH{*D+Dxq!%0=+NL-Nc^crZCP=qiPNn@r4IDf`P;}IfSPyQ=`y(Njw@c=vutrM(DmB3``L-@ zSHgVzj<`6#f85s$x4DjIdrxRtJ>f+2$6F6oFMDYFi#l^okoGzmeIq409OrkD|9i(d ze+yt89vYi}95IyRya9D+x=WNQyklqnsd0%E&z(!vF#P%lGv4Bot64^Yi>^M+XS7FG z>#NhR8+m-XSM{RwaqSbY4wdy;{>1FQf%WN_05R?;`8LarC#K z<9w(X%<=r$5d7vkn+@>?mZP=!LaXiDY~-#`S^AH|B((7Uq_;=kRdWz!Z(V!sN%Ljy zAFg)KzM+1!7^K|!7rRtIt8IMNPDAU-TC{)mDCM`e+uVDr)*0OnTGb$}{c``xo!38- zzja875cGIk_Y1L^0G*gkN_1Gx?<4EMax}Q0<$S6b%5vU06gA9WmcqnV7+iJhruZD{;dbTR7MK<)~a>@F4OfAm)RzHOU4@Wna5%U^2tF?lH8*a7#3$k00BxUBvD z*{T1_apxx9>D~7%2#T>hzjE!UCRwZQyq~kN;*-H@_jxxx}kl+ z7n7oKSCF)y^|H=l7punO`rArdZk}EFpfJW^|M$zNMTGshkINqIko%g)b=(^Yymp93 z4BMH}sMcWj2RDXp?jRZ+5YT^6_^HGDcB_O59X;>mtUoF7Z$w(_bXd;sBJ06& zt&ul5pDKp3Jby?X&ZViqT$t*~o7adBR)IUiD7Z&9{ygbpR>gMA zy=N^2dFMsnm%Q%YIWM|)xZURv_dGvzd_$zmJ9PZgnYhsZ<@oM$d9U{7-iisR0Z&f& zvBz%a)i__T*y&EkhPujdda@BU1n8#Hm^-^)#)``TpYEv^Vntge9hHY1qU1gWYNTC z(itaalM)@4^Lxm8uv}~4t-uBBhO)fS4w`0S;~#LHnuNAlJUd<;ZJc^>b}O{?k`?984H$JPoy=s3RqC+4kq$CjBqV9m2Nnu{HXTvFuy_kqhf& z$4+;sY#Qxu->S>TG7X%>%@z}ulY=-RIw{e%T&>@MH_qUKwVaO?Ls`yiq$}W&8t~{0 zxCWMY8bu-#31tce%$WGyE*dxe>ENbh0ooVlmYyd2)yi|a_EYiTQe@I<-G zWW@#UM7iZ9wYPEg1w-52AG^Y)Rgi+3u5b5T zG&EuQm6H8VHNWOoDbC2#>~rtFlJN;i4ac?Jw!GuDkmco9=e;1kaZ=Msi4Mp4UF84X zan9cYbzH%RioqObGhlel_52wye?WC~4_{bcJ1#D=AM5M*mPJEwH#oSznC#xsg)=h3 zdZnk$XjUz&!d$V%O%&{66S1v6fv<> z?UbliCzf}bY}B_~^@nE($w{k?@SK$Ba9qCwZ=Amc9p^*EP>vVUCub8I{{ZXMB(%+} zUBB&UV@RJp{m|Crh??ax&mOFoo=`Dy=kjr?C)Z2=JAZI{SJVbsHh+k_hp0)l5xBe1 z_5A#4^D?Wf<@^$l7u?RMIzv$QOlajFW0pVe)7~s1{D8yTZ4C&@H9d0<@uWnD<@_$P z9xA!kz#Dfb2HJ8kl;!y+QKaysLzz-3!n<(8NfZ%u#NqkM5}6!s3MEl+4`Kb~BRM+A zwcQe4T;m>$G^Mw&${3YxDT}s|b(Jn$r+o1FIXJn|b%*!y zb7(!c5*b@|Q%j(@ZD#?|t};h$TheN%qC=6Ye)ENaX8 zMJ3|04#IADZBn7fa%Y>KMY`{7&^~7V#-eM7v=599UfX+PxsmgnOuKLGSyF9HVD5PO z-HUS@2+~Q34$Jv{WIb4}HT0(7W5r;Wv-2o$T``TE~;lmU^}GwChK&WfkrXUA!eHDD(KxK|w{4PDRApEl!Nns<5Ae#iXW&5qL_1PvbT zyy5k7tNCS2W1T-cthXtl8Gzwh>ICSdM3?3KPO=^AJ~fFAGo$DK+~LE^AiYY7l;x zuIJ}ZhesQ>+g9pRW%()lehV8WH;s7uqGGgX{9wybA4@w){@FolHUT=R-=svB<@{c< z9xT@ydMj`RyP+)S4XP`ZLTI*0xY1U2wiJfqgfIh2{rflKMzY5ptRQX({<)5;>(0$r zZr@_jfco`j7HfXY=($Q+y<|1_X@M@z(PhxrIkSE7;uP28$7NgbySN|ddL@??7D