diff --git a/Project.toml b/Project.toml index 66bf4b1b..4b30c58e 100644 --- a/Project.toml +++ b/Project.toml @@ -33,17 +33,22 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] BenchmarkTools = "1" +Colors = "0.13" FLOWMath = "0.4" FastGaussQuadrature = "1" FiniteDiff = "2" FixedPoint = "1" ForwardDiff = "0.10" ImplicitAD = "0.3" +LaTeXStrings = "1" LineSearches = "7" MINPACK = "1" NLsolve = "4" PreallocationTools = "0.4" +ProgressMeter = "1" QuadGK = "2" +RecipesBase = "1" +Roots = "2" SIAMFANLEquations = "1" SimpleNonlinearSolve = "1, 2" SpecialFunctions = "2" diff --git a/src/analysis/setup.jl b/src/analysis/setup.jl index 56dc4e1f..7511fd01 100644 --- a/src/analysis/setup.jl +++ b/src/analysis/setup.jl @@ -118,10 +118,10 @@ function setup_analysis( # copy over operating point for f in fieldnames(typeof(operating_point)) - if f!=:units - solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + if f != :units + solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + end end -end ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- ##### if options.verbose diff --git a/src/postprocess/boundary_layer_utils.jl b/src/postprocess/boundary_layer_utils.jl index 0e80318f..a7cc56f3 100644 --- a/src/postprocess/boundary_layer_utils.jl +++ b/src/postprocess/boundary_layer_utils.jl @@ -57,21 +57,32 @@ function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot dp[1] = dp[2] end - # interpolate the lengths between control points - sum_length = 0.5 * sum(duct_panel_lengths[stag_ids]) - stag_interp = FLOWMath.linear(dp, [0.0, sum_length], 0.0) + if dp[1] == dp[2] + # we're likely in a hover-ish case and there is no stagnation point. - partial_panel_lengths = [stag_interp, sum_length - stag_interp] + s_upper = nothing + s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end:-1:1]) + split_ratio = 1.0 + stag_ids .= length(duct_panel_lengths) + else - s_upper = arc_lengths_from_panel_lengths( - [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]] - ) + # interpolate the lengths between control points + sum_length = 0.5 * sum(duct_panel_lengths[stag_ids]) + stag_interp = FLOWMath.linear(dp, [0.0, sum_length], 0.0) - s_lower = arc_lengths_from_panel_lengths( - [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]] - ) + partial_panel_lengths = [stag_interp, sum_length - stag_interp] + + s_upper = arc_lengths_from_panel_lengths( + [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]] + ) + + s_lower = arc_lengths_from_panel_lengths( + [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]] + ) + split_ratio = stag_interp / sum_length + end - return s_upper, s_lower, stag_ids, stag_interp / sum_length + return s_upper, s_lower, stag_ids, split_ratio end """ diff --git a/src/postprocess/postprocess.jl b/src/postprocess/postprocess.jl index a6267946..1762f738 100644 --- a/src/postprocess/postprocess.jl +++ b/src/postprocess/postprocess.jl @@ -517,7 +517,7 @@ function post_process( body_vortex_panels.tangent[:, 1:Int(body_vortex_panels.npanel[1])], body_vortex_panels.node[2, Int(body_vortex_panels.nnode[1])], operating_point; - verbose=verbose, + verbose=false, ) body_thrust[1] -= duct_viscous_drag diff --git a/src/postprocess/viscous_drag.jl b/src/postprocess/viscous_drag.jl index b790e522..c8f644cb 100644 --- a/src/postprocess/viscous_drag.jl +++ b/src/postprocess/viscous_drag.jl @@ -20,7 +20,7 @@ Note: if no separation occurs, the inputs are simply the final values for the bo """ function squire_young(d2, Ue, Uinf, H12) # note: formula also divides by chord, but we're going to multiply by chord later. - return 2.0 * d2 * (Ue / Uinf)^((5.0 + H12) / 2.0) + return 2.0 * d2 * FLOWMath.abs_smooth(Ue / Uinf, 2.0*eps())^((5.0 + H12) / 2.0) end """ @@ -183,13 +183,15 @@ function compute_viscous_drag_duct( ) # set up boundary layer solve parameters - boundary_layer_functions_upper = setup_boundary_layer_functions_head( - s_upper, - [0.0; Vtan_duct[stag_ids[2]:end]], - operating_point, - boundary_layer_options; - verbose=verbose, - ) + if !isnothing(s_upper) + boundary_layer_functions_upper = setup_boundary_layer_functions_head( + s_upper, + [0.0; Vtan_duct[stag_ids[2]:end]], + operating_point, + boundary_layer_options; + verbose=verbose, + ) + end # set up boundary layer solve parameters boundary_layer_functions_lower = setup_boundary_layer_functions_head( @@ -201,13 +203,15 @@ function compute_viscous_drag_duct( ) # - Set integration steps - # - # upper side - upper_steps = - set_boundary_layer_steps( - boundary_layer_options.n_steps, - boundary_layer_options.first_step_size, - s_upper[end] - boundary_layer_options.offset, - ) .+ boundary_layer_options.offset + if !isnothing(s_upper) + # upper side + upper_steps = + set_boundary_layer_steps( + boundary_layer_options.n_steps, + boundary_layer_options.first_step_size, + s_upper[end] - boundary_layer_options.offset, + ) .+ boundary_layer_options.offset + end # lower side lower_steps = @@ -219,19 +223,26 @@ function compute_viscous_drag_duct( # - Get drag coeffients - # - cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head( - upper_steps, - exit_radius, - operating_point, - boundary_layer_functions_upper, - (; - boundary_layer_options.rk, - boundary_layer_options.separation_criteria, - separation_allowance=boundary_layer_options.separation_allowance_upper, - separation_penalty=boundary_layer_options.separation_penalty_upper, - ); - verbose=verbose, - ) + if !isnothing(s_upper) + cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head( + upper_steps, + exit_radius, + operating_point, + boundary_layer_functions_upper, + (; + boundary_layer_options.rk, + boundary_layer_options.separation_criteria, + separation_allowance=boundary_layer_options.separation_allowance_upper, + separation_penalty=boundary_layer_options.separation_penalty_upper, + ); + verbose=verbose, + ) + else + cdc_upper = 0.0 + usol_upper = -ones(eltype(Vtan_duct), 3) + stepsol_upper = -1 + s_sep_upper = -1.0 + end cdc_lower, usol_lower, stepsol_lower, s_sep_lower = compute_single_side_drag_coefficient_head( lower_steps, diff --git a/src/preprocess/geometry/rotor_geometry.jl b/src/preprocess/geometry/rotor_geometry.jl index b8a85595..1690ad1f 100644 --- a/src/preprocess/geometry/rotor_geometry.jl +++ b/src/preprocess/geometry/rotor_geometry.jl @@ -242,10 +242,10 @@ function get_blade_ends_from_body_geometry!( if !silence_warnings for (irotor, (R, r)) in enumerate(zip(Rtip, Rhub)) if r !== rp_centerbody_coordinates[2, ihub[irotor]] - @warn "Overwriting Rhub for rotor $(irotor) to place it at the centerbody wall. Moving from $(r) to $(rp_centerbody_coordinates[2, ihub[irotor]])" + @warn "Overwriting Rhub for rotor $(ForwardDiff.value(irotor)) to place it at the centerbody wall. Moving from $(ForwardDiff.value(r)) to $(ForwardDiff.value(rp_centerbody_coordinates[2, ihub[irotor]]))" end if R !== rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor] - @warn "Overwriting Rtip for rotor $(irotor) to place it at the correct tip gap relative to the casing wall. Moving from $(R) to $(rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor])" + @warn "Overwriting Rtip for rotor $(ForwardDiff.value(irotor)) to place it at the correct tip gap relative to the casing wall. Moving from $(ForwardDiff.value(R)) to $(ForwardDiff.value(rp_duct_coordinates[2, iduct[irotor]] .- tip_gaps[irotor]))" end end end diff --git a/src/utilities/bookkeeping.jl b/src/utilities/bookkeeping.jl index 8125951e..21919a4a 100644 --- a/src/utilities/bookkeeping.jl +++ b/src/utilities/bookkeeping.jl @@ -65,7 +65,6 @@ function get_problem_dimensions(paneling_constants::PanelingConstants) nbe = nws - 1 # number of body panels - # TODO: need number of casing and centerbody nodes separately as well ncp = nduct_inlet ncbp = ncenterbody_inlet # add rest of panels mutual between centerbody and duct diff --git a/src/utilities/inputs.jl b/src/utilities/inputs.jl index ec82614b..96a919a4 100644 --- a/src/utilities/inputs.jl +++ b/src/utilities/inputs.jl @@ -277,7 +277,7 @@ function Rotor( fliplift; i_know_what_im_doing=false, ) - if length(findall(t -> t > pi / 2, twists)) > 2 + if length(findall(t -> t > 1.75, twists)) > 2 @warn "It looks like your input twist angles may be in degrees. Note that the required units for twist are radians. Converting to radians for you (set the `i_know_what_im_doing` keyword argument to true to disable automatic conversion)." twists .*= pi / 180.0 end diff --git a/src/utilities/misc.jl b/src/utilities/misc.jl index 00e6b0a3..2dfef3ac 100644 --- a/src/utilities/misc.jl +++ b/src/utilities/misc.jl @@ -36,11 +36,21 @@ end Formatted printing when you have lots of debugging to do and want things to line up nicely. """ -function printdebug(variable_name, variable, nspaces=4) - @printf "%*s %-14s %2.12f\n" nspaces " " variable_name variable +function printdebug(variable_name, variable, nspaces=4; colw=16) + if length(variable) == 1 + s = @sprintf "%*s %-*s %2.12f" nspaces " " colw variable_name variable[1] + println(rstrip(s, ['0', ' '])) + else + @printf "%*s %-*s\n" nspaces " " colw variable_name + for v in variable + s = @sprintf "%*s %-*s %2.12f" (nspaces) " " colw "" v + println(rstrip(s, ['0', ' '])) + end + end # println("\t$(variable_name)" * v) return nothing end + """ dot(A, B) = sum(a * b for (a, b) in zip(A, B)) @@ -150,32 +160,33 @@ function reset_containers!(c; exception_keys=[]) end """ - promote_ducted_rotor_type(ducted_rotor) + promote_ducted_rotor_type(ducted_rotor, operating_point) Convenience function for promoting types based on any potential elements of the ducted_rotor object dependent on optimization design variables. # Arguments - `ducted_rotor::DuctedRotor` : the ducted_rotor input +- `operating_point::OperatingPoint` : the operating_point input # Returns - `TP::Type` : the promoted type """ -function promote_ducted_rotor_type(p) +function promote_ducted_rotor_type(d,o) return promote_type( - eltype(p.duct_coordinates), - eltype(p.centerbody_coordinates), - eltype(p.operating_point.Vinf), - eltype(p.operating_point.rhoinf), - eltype(p.operating_point.muinf), - eltype(p.operating_point.asound), - eltype(p.operating_point.Omega), - eltype(p.rotor.B), - eltype(p.rotor.rotorzloc), - eltype(p.rotor.r), - eltype(p.rotor.Rhub), - eltype(p.rotor.Rtip), - eltype(p.rotor.chords), - eltype(p.rotor.twists), + eltype(d.duct_coordinates), + eltype(d.centerbody_coordinates), + eltype(o.Vinf), + eltype(o.rhoinf), + eltype(o.muinf), + eltype(o.asound), + eltype(o.Omega), + eltype(d.rotor.B), + eltype(d.rotor.rotorzloc), + eltype(d.rotor.r), + eltype(d.rotor.Rhub), + eltype(d.rotor.Rtip), + eltype(d.rotor.chords), + eltype(d.rotor.twists), ) end diff --git a/src/utilities/options.jl b/src/utilities/options.jl index fa209747..a7891cf7 100644 --- a/src/utilities/options.jl +++ b/src/utilities/options.jl @@ -660,7 +660,7 @@ end # Fields: - `model_drag::Tb=false` : flag to turn on viscous drag approximation -- `n_steps::Int = Int(2e2)` : number of steps to use in boundary layer integration +- `n_steps::Int = Int(5e2)` : number of steps to use in boundary layer integration - `first_step_size::Float = 1e-6` : size of first step in boundary layer integration - `offset::Float = 1e-3` : size of offset for (where to initialize) boundary layer integration - `rk::Function = RK4` : solver to use for boundary layer integration (RK4 or RK2 available) @@ -672,7 +672,7 @@ end """ @kwdef struct HeadsBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp,Ts} <: BoundaryLayerOptions model_drag::Tb=false - n_steps::Ti = Int(2e2) + n_steps::Ti = Int(5e2) first_step_size::Tf = 1e-6 offset::To = 1e-3 rk::Tfun = RK2 diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl index 9d3637f2..40d129b9 100644 --- a/src/visualization/plot_recipes.jl +++ b/src/visualization/plot_recipes.jl @@ -160,9 +160,8 @@ end ylim --> (0.0, maximum(bvp.node[2, :]) * 1.05) if discrete_labels - xticks --> DuctAPE.determine_geometry_xlabels( - bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits - ) + xticks --> + determine_geometry_xlabels(bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits) yticks --> determine_geometry_ylabels(bvp, rsp; tol=yticktol, tickdigits=tickdigits) xgrid --> true ygrid --> true @@ -406,8 +405,10 @@ function determine_geometry_ylabels(bvp, rsp; tol=1e-2, tickdigits=2) # Rotor for r in 1:Int(rsp.nbodies[]) - push!(yt, rsp.node[2, Int(rsp.endnodeidxs[r, r])]) - push!(yl, @sprintf "%3.*f" tickdigits rsp.node[2, Int(rsp.endnodeidxs[r, r])]) + for i in 1:2 + push!(yt, rsp.node[2, Int(rsp.endnodeidxs[i, r])]) + push!(yl, @sprintf "%3.*f" tickdigits rsp.node[2, Int(rsp.endnodeidxs[i, r])]) + end end # Centerbody @@ -416,18 +417,22 @@ function determine_geometry_ylabels(bvp, rsp; tol=1e-2, tickdigits=2) yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits) end - # Duct front/back - for f in [findmax, findmin] - tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) - tmpy = bvp.node[2, tmp[2]] - yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits) - end + # # Duct front/back + # for f in [findmax, findmin] + # tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) + # tmpy = bvp.node[2, tmp[2]] + # yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits) + # end - # Duct inner/outer - for f in [maximum, minimum] - tmp = f(bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) - yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits) - end + # Duct Exit Radius + tmpy = bvp.node[2, 1] + yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits) + + # # Duct inner/outer + # for f in [maximum, minimum] + # tmp = f(bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) + # yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits) + # end return (yt, yl) end @@ -712,7 +717,7 @@ end # Plot specific values aspect_ratio --> 1 - p, _ = DuctAPE.calculate_streamlines( + p, _ = calculate_streamlines( bvp, bvs, wvp,