diff --git a/Project.toml b/Project.toml index 3183a36c..d2822fff 100644 --- a/Project.toml +++ b/Project.toml @@ -17,6 +17,7 @@ LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" +OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 9dc42366..a904a314 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -44,6 +44,9 @@ using ForwardDiff # used for jacobian for newton solver # For boundary layer stuff using Roots +using OrdinaryDiffEq: + ODEProblem, terminate!, ContinuousCallback, RadauIIA5 +using OrdinaryDiffEq: solve as diffeq_solve # - Utility Packages - # using FLOWMath # used for various items, mostly interpolation @@ -85,6 +88,8 @@ export ChainSolverOptions, FixedPointOptions, CSORSolverOptions export BoundaryLayerOptions +export RK2, RK4 +export RadauIIA5, RadauIIA3 # - Preprocess - # export setup_analysis @@ -190,6 +195,7 @@ include("postprocess/viscous_drag.jl") # include("visualization/plot_recipe_defaults.jl") include("visualization/plot_recipes.jl") include("visualization/calculate_streamlines.jl") + include("visualization/calculate_velocity_field.jl") include("visualization/convenience_plots.jl") ##### ----- DEBUGGING ----- ##### diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index 5c5322b1..427c9bd7 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -66,7 +66,7 @@ function analyze( #TODO: probably just call the post-process function directly and return a reset_container! of the output if return_inputs return [],#zero_outputs(), - (; solve_parameter_tuple..., airfoils, idmaps, panels, problem_dimensions), + (; solve_parameter_tuple..., airfoils, idmaps, panels, problem_dimensions, reference_parameters, multi_point = [operating_point]), false else return [],#zero_outputs(), @@ -410,6 +410,8 @@ function analyze( linsys=(; solve_parameter_tuple.linsys..., A_bb_LU), idmaps, problem_dimensions, + reference_parameters, + multi_point=operating_point, ), options.solver_options.converged else diff --git a/src/postprocess/boundary_layer_green.jl b/src/postprocess/boundary_layer_green.jl index e5535776..59623061 100644 --- a/src/postprocess/boundary_layer_green.jl +++ b/src/postprocess/boundary_layer_green.jl @@ -23,11 +23,10 @@ function calculate_radius_of_curvature(s, controlpoint, ss) xdot = FLOWMath.derivative.(Ref(x_of_s), ss) xdotsp = Akima_smooth(s, FLOWMath.derivative.(Ref(x_of_s), s)) xddot = FLOWMath.derivative.(Ref(xdotsp), ss) - # xddot = FLOWMath.second_derivative.(Ref(x_of_s), ss) + ydot = FLOWMath.derivative.(Ref(y_of_s), ss) ydotsp = Akima_smooth(s, FLOWMath.derivative.(Ref(y_of_s), s)) yddot = FLOWMath.derivative.(Ref(ydotsp), ss) - # yddot = FLOWMath.second_derivative.(Ref(y_of_s), ss) # assemble the numerator and denominator of the radius of curvature expression num = (xdot .^ 2 .+ ydot .^ 2) .^ (1.5) @@ -94,30 +93,12 @@ function setup_boundary_layer_functions_green( # radii of curvature radius_of_curvature(ss) = calculate_radius_of_curvature(s, duct_control_points, ss) - # local mach number - edge_mach(ss) = calculate_mach(edge_velocity(ss), operating_point.asound[]) - - # local density - Pe(ss) = static_pressure(operating_point.Ptot[], edge_mach(ss)) - edge_density(ss) = static_density(Pe(ss), operating_point.asound[]) - - # local viscosity - function Te(ss) - return convert_temperature_to_kelvin.( - Ref(operating_point.units), - static_temperature.(Ref(operating_point.Ttot[]), edge_mach(ss)), - ) - end - function edge_viscosity(ss) - return convert_viscosity.(Ref(operating_point.units), sutherlands_law(Te(ss))) - end - return (; edge_velocity, - edge_mach, + edge_mach=x -> 0.0, edge_acceleration, - edge_density, - edge_viscosity, + edge_density=x -> operating_point.rhoinf[], + edge_viscosity=x -> operating_point.muinf[], r_coords, radial_derivative, radius_of_curvature, @@ -148,187 +129,175 @@ end """ Rename of `calculate_H12bar0` used for initialization to avoid confusion """ -function H12bar_init(Cf0, M) - return calculate_H12bar0(Cf0, M) +function H12bar_init(Cf0) + return calculate_H12bar0(Cf0) end """ Wrapper of `calculate_CEeq` used for initialization to avoid confusion (`lambda` set to 1) """ -function CE_init(Ctaueq0, M, Cf0) - return calculate_CEeq(Ctaueq0, M, 1, Cf0) +function CE_init(Ctaueq0, Cf0) + return calculate_CEeq(Ctaueq0, 1, Cf0) end """ - initialize_turbulent_boundary_layer_states_green( - s_init, r_init, Ue, M, rhoe, mue; verbose=false - ) + initialize_green_states(boundary_layer_functions, s_init; verbose=false) Initialize the states for the turbulent boundary layer solve. # Arguments: +- `boundary_layer_functions::NamedTuple` : - `s_init::Float` : surface length starting point -- `r_init::Float` : surface radial position at starting point -- `Ue::Float` : edge velocity at starting point -- `M::Float` : edge Mach at starting point -- `rhoe::Float` : edge density at starting point -- `mue::Float` : edge dynamic viscosity at starting point # Returns: - `initial_states::Vector{Float}` : Initial values for the states: r*delta2, Hbar12, CE - `Cf_init::Vector{Float}` : initial value for friction coefficietn Cf (used for determining separateion) - `H12_init::Vector{Float}` : initial value for shape factor H12 (used in viscous drag calculation) """ -function initialize_turbulent_boundary_layer_states_green( - s_init, r_init, Ue, M, rhoe, mue; verbose=false -) - verbose && println("INITIALIZATION") +function initialize_green_states(boundary_layer_functions, s_init; verbose=false) + (; r_coords, edge_density, edge_velocity, edge_viscosity) = boundary_layer_functions + + r_init = r_coords(s_init) + Ue = edge_velocity(s_init) + rhoe = edge_density(s_init) + mue = edge_viscosity(s_init) - # println("s_init: ", s_init) - # println("r_init: ", r_init) - # println("Ue: ", Ue) - # println("M: ", M) - # println("rho: ", rhoe) - # println("mue: ", mue) + verbose && println("INITIALIZATION") # - Initialize States - # # initialize momentum thickness (d2) using flat plate schlichting model Rex = calculate_Re(rhoe, Ue, s_init, mue) verbose && printdebug("Rex: ", Rex) + #d2 = ksmax([1e-3;d2_init(s_init, Rex)],10000) d2 = d2_init(s_init, Rex) verbose && printdebug("d2: ", d2) # initialize H12bar (compressible shape factor) using _0 equations Red2 = calculate_Re(rhoe, Ue, d2, mue) verbose && printdebug("Red2: ", Red2) - Cf0 = calculate_Cf0(Red2, M) + Cf0 = calculate_Cf0(Red2) verbose && printdebug("Cf0: ", Cf0) - H12bar0 = H12bar_init(Cf0, M) + H12bar0 = H12bar_init(Cf0) verbose && printdebug("H12bar0: ", H12bar0) # initialize the entrainment coefficient (CE) using equilibrium equations H1 = calculate_H1(H12bar0) verbose && printdebug("H1: ", H1) - H12 = calculate_H12(H12bar0, M) + H12 = H12bar0 verbose && printdebug("H12: ", H12) Cf = calculate_Cf(H12bar0, H12bar0, Cf0) verbose && printdebug("Cf: ", Cf) - d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar0, H12, Cf, M) + d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar0, H12, Cf) verbose && printdebug("d2dUedsUeeq0: ", d2dUedsUeeq0) CEeq0 = calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) verbose && printdebug("CEeq0: ", CEeq0) - Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0, M) + Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0) verbose && printdebug("Ctaueq0: ", Ctaueq0) - CE = CE_init(Ctaueq0, M, Cf0) + CE = ksmax([1.0; CE_init(Ctaueq0, Cf0)]) #TODO: HACK for smooth start verbose && printdebug("CE: ", CE) # initial states - return [r_init * d2, H12bar0, CE], Cf, H12 + return [r_init * d2, H12bar0, CE], Cf end #---------------------------------# # Intermediate functions # #---------------------------------# -""" - Fc(M) = sqrt(1.0 + 0.2 * M^2) -""" -function Fc(M) - return sqrt(1.0 + 0.2 * M^2) -end - -""" - FR(M) = 1.0 + 0.056 * M^2 -""" -function FR(M) - return 1.0 + 0.056 * M^2 -end - """ calculate_Re(ρ, V, L, μ) Calculate Reynolds number """ function calculate_Re(rhoe, Ue, d2, mue) - return FLOWMath.ksmax([21.0; rhoe * Ue * d2 / mue]) + # limit Red2 to 11 since there is a discontinuity in Cf_θ just below 11. + # limit Red2 to 18 since there is a discontinuity in H12bar0 just below 18. + return FLOWMath.ksmax([18.0; rhoe * Ue * d2 / mue]) end """ - calculate_Cf0(Red2, M) + calculate_Cf0(Red2) -Calculate Cf₀ value based on momentum thickness Reynolds number (Red2) and edge mach number (M). +Calculate Cf₀ value based on momentum thickness Reynolds number (Red2). +Incompressible version of eqn A-14 in Green_1977. """ -function calculate_Cf0(Red2, M) - # if (0.01013 / (log(10, FR(M) * Red2) - 1.02) - 0.00075) / Fc(M) < 0.0 - # println("Red2: ", Red2) - # println("M: ", M) - # println("FR(M): ", FR(M)) - # println("Log: ", log(10, FR(M) * Red2)) - # println("Log-1.02: ", log(10, FR(M) * Red2) - 1.02) - # println("Cf0: ", (0.01013 / (log(10, FR(M) * Red2) - 1.02) - 0.00075) / Fc(M)) - # end - - # return FLOWMath.ksmax([0.0;(0.01013 / (log(10, FLOWMath.ksmax([1.0;FR(M) * Red2])) - 1.02) - 0.00075) / Fc(M)]) - - return (0.01013 / (log(10, FR(M) * Red2) - 1.02) - 0.00075) / Fc(M) +function calculate_Cf0(Red2) + #NOTE: Red2 must be > ~11 + return (0.01013 / (log(10, Red2) - 1.02) - 0.00075) end """ - calculate_H12bar0(Cf0, M) + calculate_H12bar0(Cf0) + +Incompressible version of eqn A-15 in Green_1977. """ -function calculate_H12bar0(Cf0, M) - return 1.0 / (1.0 - 6.55 * sqrt(Cf0 / 2.0 * (1.0 + 0.04 * M^2))) +function calculate_H12bar0(Cf0) + #NOTE: Red2 must be > ~17.5 + return 1.0 / (1.0 - 6.55 * sqrt(Cf0 / 2.0)) end """ calculate_Cf(H12bar, H12bar0, Cf0) Calculate friction coefficient. +From eqn A-16 in Green_1977. """ function calculate_Cf(H12bar, H12bar0, Cf0) + #NOTE: H12bar0 CANNOT be zero + #NOTE: H12bar/H12bar0 cannot be 0.4 + # return Cf0 * (0.9 / (FLOWMath.ksmax([0.4; H12bar / H12bar0]) - 0.4) - 0.5) return Cf0 * (0.9 / (H12bar / H12bar0 - 0.4) - 0.5) end """ - calculate_H12(H12bar, M, Pr=1.0) -""" -function calculate_H12(H12bar, M, Pr=1.0) - return (H12bar + 1.0) * (1.0 + Pr^(1.0 / 3.0) * M^2 / 5) - 1.0 -end + calculate_Ctau(CE, Cf0) +Incompressible version of eqn A-20 from Green_1977. """ - calculate_Ctau(CE, Cf0, M) -""" -function calculate_Ctau(CE, Cf0, M) - return (0.024 * CE + 1.2 * CE^2 + 0.32 * Cf0) * (1.0 + 0.1 * M^2) +function calculate_Ctau(CE, Cf0) + return 0.024 * CE + 1.2 * CE^2 + 0.32 * Cf0 end """ calculate_F(CE, Cf0) + +From eqn A-21 in Green_1977. """ function calculate_F(CE, Cf0) - return (0.02 * CE + CE^2 + 0.8 * Cf0 / 3) / (0.01 + CE) + #NOTE: CE CANNOT be -0.01 TODO(is it even possible for it to be negative?) + return (0.02 * CE + CE^2 + 0.8 * Cf0 / 3.0) / (0.01 + CE) end """ calculate_H1(H12bar) + +From eqn A-18 in Green_1977. """ function calculate_H1(H12bar) + #NOTE: H12bar CANNOT be 1.0 return 3.15 + 1.72 / (H12bar - 1.0) - 0.01 * (H12bar - 1.0)^2 end """ calculate_dH12bardH1(H12bar) + +From eqn A-19 in Green_1977. """ function calculate_dH12bardH1(H12bar) + #NOTE: H12bar CANNOT be -3.414004962442103 TODO(can H12bar even be negative?) return -(H12bar - 1.0)^2 / (1.72 + 0.02 * (H12bar - 1.0)^3) end """ calculate_richardson_number(H12bar, d2, H12, H1, R) + +From eqn A-22 in Green_1977. """ function calculate_richardson_number(H12bar, d2, H12, H1, R) + #NOTE: R CANNOT be zero TODO(is it possible to get zero radius of longitudinal curvature?) + #NOTE: H12+H1 CANNOT be zero TODO(is it possible for H1 or H12bar to be < 0?) + #NOTE: H1/H12bar CANNOT be -0.3 TODO(is it possible for H1 or H12bar to be < 0?) return 2.0 * d2 / (3.0 * R) * (H12 + H1) * (H1 / H12bar + 0.3) end @@ -336,10 +305,13 @@ end calculate_beta(Ri; hardness=100.0) Sigmoind blended version of piecewise function: + | 7.0 if Ri > 0 β = | | 4.5 if Ri < 0 +From eqn A-23 in Green_1977. + # Arguments: - `Ri::float` : Richardson number @@ -354,24 +326,23 @@ function calculate_beta(Ri, hardness=100.0) end """ - longitudinal_curvature_influence(M, Ri) + longitudinal_curvature_influence(Ri) + +Incompressible version of eqn A-24 in Green_1977. """ -function longitudinal_curvature_influence(M, Ri) - return 1 + calculate_beta(Ri) * (1.0 + M^2 / 5) * Ri +function longitudinal_curvature_influence(Ri) + return 1 + calculate_beta(Ri) * Ri end """ lateral_strain_influence(H12bar, d2, H12, H1, r, drds) -""" -function lateral_strain_influence(H12bar, d2, H12, H1, r, drds) - return 1.0 - 7.0 / 3.0 * (H1 / H12bar + 1.0) * (H12 + H1) * d2 * drds / r -end +From eqn A-25 in Green_1977. """ - dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) -""" -function dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) - return 1.0 + 7.0 / 3.0 * M^2 * (H1 / H12bar + 1.0) * (H12 + H1) * d2 * dUeds / Ue +function lateral_strain_influence(H12bar, d2, H12, H1, r, drds) + #NOTE: H12bar CANNOT be zero TODO(is it possible?) + #NOTE: r cannot be zero TODO(cannot apply this to body of rotation?) + return 1.0 - 7.0 / 3.0 * (H1 / H12bar + 0.3) * (H12 + H1) * d2 * drds / r end """ @@ -384,40 +355,51 @@ function calculate_lambda(args...) end """ - calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) + calculate_d2dUedsUeeq0(H12bar, H12, Cf) + +Incompressible version of eqn A-28 in Green_1977. """ -function calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) - return 1.25 * (Cf / 2.0 - ((H12bar - 1.0) / (6.432 * H12bar))^2 / (1.0 + 0.04 * M^2)) / - H12 +function calculate_d2dUedsUeeq0(H12bar, H12, Cf) + #NOTE: H12 CANNOT be zero TODO(is it possible?) + #NOTE: H12bar CANNOT be zero TODO(is it possible?) + return 1.25 * (Cf / 2.0 - ((H12bar - 1.0) / (6.432 * H12bar))^2) / H12 end """ calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) -Note: applies smooth-max to makes sure this value stays non-negative +From eqn A-29 in Green_1977. """ function calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) - # apply smooth-max to make sure we stay non-negative. - return FLOWMath.ksmax([0.0; H1 * (Cf / 2.0 - (H12 + 1.0) * d2dUedsUeeq0)]) + return ksmax([0.0; H1 * (Cf / 2.0 - (H12 + 1.0) * d2dUedsUeeq0)]) end """ - calculate_Ctaueq0(CEeq0, Cf0, M) + calculate_Ctaueq0(CEeq0, Cf0) + +Incompressible version of eqn A-30 in Green_1977. """ -function calculate_Ctaueq0(CEeq0, Cf0, M) - return (0.24 * CEeq0 + 1.2 * CEeq0^2 + 0.32 * Cf0) * (1.0 + 0.1 * M^2) +function calculate_Ctaueq0(CEeq0, Cf0) + return 0.24 * CEeq0 + 1.2 * CEeq0^2 + 0.32 * Cf0 end """ - calculate_CEeq(Ctaueq0, M, lambda, Cf0) + calculate_CEeq(Ctaueq0, lambda, Cf0) + +Incompressible version of eqns A-31 and A-32 in Green_1977. """ -function calculate_CEeq(Ctaueq0, M, lambda, Cf0) - c = FLOWMath.ksmax([0.0; Ctaueq0 / ((1.0 + 0.1 * M^2) * lambda^2) - 0.32 * Cf0]) +function calculate_CEeq(Ctaueq0, lambda, Cf0) + #NOTE: c CANNOT be < -1.2e-4 + # c = Ctaueq0 / lambda^2 - 0.32 * Cf0 + #TODO: figure out how to fix this upstream rather than adding this patch here. + c = FLOWMath.ksmax([-1.2e-4; Ctaueq0 / lambda^2 - 0.32 * Cf0]) return sqrt(c / 1.2 + 0.0001) - 0.01 end """ calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) + +From eqn A-33 in Green_1977. """ function calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) return (Cf / 2.0 - CEeq / H1) / (H12 + 1.0) @@ -428,21 +410,21 @@ end #---------------------------------# """ - boundary_layer_residual_green(y, s, parameters) + boundary_layer_residual_green(y, parameters, s) Out-of-place version of `boundary_layer_residual_green!` """ -function boundary_layer_residual_green(y, s, parameters) +function boundary_layer_residual_green(y, parameters, s) dy = similar(y) .= 0 - return boundary_layer_residual_green!(dy, y, s, parameters) + return boundary_layer_residual_green!(dy, y, parameters, s) end """ - boundary_layer_residual_green!(dy, y, s, parameters) + boundary_layer_residual_green!(dy, y, parameters, s) Calculate dy give the current states, y, the input position, s, and various parameters. """ -function boundary_layer_residual_green!(dy, y, s, parameters) +function boundary_layer_residual_green!(dy, y, parameters, s) (; r_coords, edge_velocity, @@ -465,8 +447,6 @@ function boundary_layer_residual_green!(dy, y, s, parameters) verbose && printdebug("rhoe: ", rhoe) mue = edge_viscosity(s) verbose && printdebug("mue: ", mue) - M = edge_mach(s) - verbose && printdebug("M: ", M) dUeds = edge_acceleration(s) verbose && printdebug("dUeds: ", dUeds) R = radius_of_curvature(s) @@ -488,16 +468,16 @@ function boundary_layer_residual_green!(dy, y, s, parameters) Red2 = calculate_Re(rhoe, Ue, d2, mue) verbose && printdebug("Red2: ", Red2) - Cf0 = calculate_Cf0(Red2, M) + Cf0 = calculate_Cf0(Red2) verbose && printdebug("Cf0: ", Cf0) - H12bar0 = calculate_H12bar0(Cf0, M) + H12bar0 = calculate_H12bar0(Cf0) verbose && printdebug("H12bar0: ", H12bar0) Cf = calculate_Cf(H12bar, H12bar0, Cf0) verbose && printdebug("Cf: ", Cf) - H12 = calculate_H12(H12bar, M) + H12 = H12bar verbose && printdebug("H12: ", H12) dH12bardH1 = calculate_dH12bardH1(H12bar) #yes, this should be negative according to example @@ -508,22 +488,22 @@ function boundary_layer_residual_green!(dy, y, s, parameters) F = calculate_F(CE, Cf0) verbose && printdebug("F: ", F) - d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) + d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar, H12, Cf) verbose && printdebug("d2dUedsUeeq0: ", d2dUedsUeeq0) CEeq0 = calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) verbose && printdebug("CEeq0: ", CEeq0) - Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0, M) + Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0) verbose && printdebug("Ctaueq0: ", Ctaueq0) - Ctau = calculate_Ctau(CE, Cf0, M) + Ctau = calculate_Ctau(CE, Cf0) verbose && printdebug("Ctau: ", Ctau) if parameters.lambda if parameters.longitudinal_curvature Ri = calculate_richardson_number(H12bar, d2, H12, H1, r) - l1 = longitudinal_curvature_influence(M, Ri) + l1 = longitudinal_curvature_influence(Ri) else l1 = 1 end @@ -534,19 +514,13 @@ function boundary_layer_residual_green!(dy, y, s, parameters) l2 = 1 end - if parameters.dilation - l3 = dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) - else - l3 = 1 - end - - lambda = calculate_lambda(l1, l2, l3) + lambda = calculate_lambda(l1, l2) else lambda = 1 end verbose && printdebug("lambda: ", lambda) - CEeq = calculate_CEeq(Ctaueq0, M, lambda, Cf0) + CEeq = calculate_CEeq(Ctaueq0, lambda, Cf0) verbose && printdebug("CEeq: ", CEeq) d2dUedsUeeq = calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) @@ -554,24 +528,46 @@ function boundary_layer_residual_green!(dy, y, s, parameters) # - system of equations - # - # momentum - dy[1] = r * Cf / 2.0 - (H12 + 2.0 - M^2) * rd2 * dUeds / Ue + # momentum: Incompressible version of eqn A-8 in Green_1977 + dy[1] = r * Cf / 2.0 - (H12 + 2.0) * rd2 * dUeds / Ue - # entrainment + # entrainment: eqn A-9 in Green_1977 dy[2] = dH12bardH1 * (CE - H1 * (Cf / 2.0 - (H12 + 1.0) * d2 * dUeds / Ue)) / d2 - # lag + # lag: Incompressible version of eqn A-10 in Green_1977 dy[3] = F / d2 * ( 2.8 / (H12 + H1) * (sqrt(Ctaueq0) - lambda * sqrt(Ctau)) + d2dUedsUeeq - - d2 * dUeds / Ue * (1.0 + 0.075 * M^2 * ((1.0 + 0.2 * M^2) / (1.0 + 0.1 * M^2))) + d2 * dUeds / Ue ) verbose && println(" Residuals:") verbose && printdebug("d(rd2)/ds: ", dy[1]) verbose && printdebug("d(H12bar)/ds: ", dy[2]) verbose && printdebug("d(CE)/ds: ", dy[3]) - return dy, [Cf; H12] + return dy +end + +function update_Cf(state, step, parameters) + + # Unpack State + rd2, H12bar, CE = state + + # Unpack parameters + (; r_coords, edge_velocity, edge_density, edge_viscosity) = parameters + + # - Intermediate Calculations - # + + # Get Cf0 + Red2 = calculate_Re( + edge_density(step), edge_velocity(step), rd2 / r_coords(step), edge_viscosity(step) + ) + Cf0 = calculate_Cf0(Red2) + + # Get H12bar0 + H12bar0 = calculate_H12bar0(Cf0) + + return calculate_Cf(H12bar, H12bar0, Cf0) end """ @@ -587,11 +583,12 @@ Integrate the turbulent boundary layer using a Runge-Kutta method. - `parameters::NamedTuple` : boundary layer solve options and other parameters """ function solve_green_boundary_layer!( - f, rk, initial_states, steps, parameters; verbose=false + ::RK, ode, initial_states, steps, parameters; verbose=false ) + f = boundary_layer_residual_green # Unpack States and variables for viscous drag - u0, Cf0, H12_0 = initial_states + u0, Cf0 = initial_states # Initilization separate flags and outputs sep = [false] @@ -601,12 +598,10 @@ function solve_green_boundary_layer!( # TODO; put these in a cache that gets updated in place. us = zeros(eltype(u0), length(u0), length(steps)) Cfs = zeros(eltype(u0), length(steps)) - H12s = zeros(eltype(u0), length(steps)) # Initialize intermediate states and outputs us[:, 1] = u0 Cfs[1] = Cf0 - H12s[1] = H12_0 # Take Runge-Kutta Steps until either the end of the surface or separation is detected for i in 1:(length(steps) - 1) @@ -619,13 +614,11 @@ function solve_green_boundary_layer!( end # take step - us[:, i + 1], aux = rk( - f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters - ) + us[:, i + 1] = ode(f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters) - Cfs[i + 1], H12s[i + 1] = aux + Cfs[i + 1] = update_Cf(us[:, i + 1], steps[i + 1], parameters) - sepid[1] = i + 1 + sepid[1] = i if Cfs[i + 1] <= 0.0 sep[1] = true break @@ -635,26 +628,91 @@ function solve_green_boundary_layer!( if sep[1] == true # - Interpolate to find actual s_sep - # + if sepid[1] == 1 + usep = us[:, 1] + s_sep = steps[1] - usep = [ - FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[1, (sepid[] - 1):sepid[]], 0.0) - FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[2, (sepid[] - 1):sepid[]], 0.0) - FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[3, (sepid[] - 1):sepid[]], 0.0) - ] + usol = us + stepsol = steps - H12sep = FLOWMath.linear( - Cfs[(sepid[] - 1):sepid[]], H12s[(sepid[] - 1):sepid[]], 0.0 - ) - - s_sep = FLOWMath.linear( - Cfs[(sepid[] - 1):sepid[]], steps[(sepid[] - 1):sepid[]], 0.0 - ) + else + usep = [ + FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], us[1, (sepid[] - 1):sepid[]], 0.0 + ) + FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], us[2, (sepid[] - 1):sepid[]], 0.0 + ) + FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], us[3, (sepid[] - 1):sepid[]], 0.0 + ) + ] + + s_sep = FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], steps[(sepid[] - 1):sepid[]], 0.0 + ) + + stepsol = steps[1:sepid[]] + usol = us[:, 1:sepid[]] + end else usep = us[:, end] - H12sep = H12s[end] s_sep = steps[end] + + usol = us + stepsol = steps end # return states at separate, and separation shape factor, and surface length at separation - return usep, FLOWMath.ksmax([1.0; H12sep]), s_sep, sepid + return usep, usep[2], s_sep, usol, stepsol +end + +""" + solve_green_boundary_layer!(::DiffEq, ode, initial_states, steps, parameters; verbose=false) + +Integrate the turbulent boundary layer using a Runge-Kutta method. + +# Arguments: +- `ode::function_handle` : ODE method to use (one of the DifferentialEquations.jl options) +- `initial_states::Float` : initial states +- `steps::Vector{Float}` : steps for integration +- `parameters::NamedTuple` : boundary layer solve options and other parameters +""" +function solve_green_boundary_layer!( + ::DiffEq, ode, initial_states, steps, parameters; verbose=false +) + f = boundary_layer_residual_green! + + prob = ODEProblem(f, initial_states[1], [steps[1], steps[end]], parameters) + + # set up separation termination conditions + function condition(u, t, integrator) + return update_Cf(u, t, parameters) + end + + function affect!(integrator) + return terminate!(integrator) + end + + cb = ContinuousCallback(condition, affect!) + + sol = diffeq_solve( + prob, + ode(); + callback=cb, + abstol=eps(), + dtmax=1e-4, + # alg_hints=[:stiff], + # dt=parameters.first_step_size, + verbose=true, + ) + + usep = sol.u[end] + Hsep = usep[2] + s_sep = sol.t[end] + usol = reduce(hcat, (sol.u)) + stepsol = sol.t + + # return states at separate, and separation shape factor, and surface length at separation + return usep, Hsep, s_sep, usol, stepsol end diff --git a/src/postprocess/boundary_layer_green_compressible.jl b/src/postprocess/boundary_layer_green_compressible.jl new file mode 100644 index 00000000..d979a497 --- /dev/null +++ b/src/postprocess/boundary_layer_green_compressible.jl @@ -0,0 +1,648 @@ +""" + calculate_radius_of_curvature(s, xy, ss) + +Determine the radius of curvature. + +(see https://en.wikipedia.org/wiki/Radius_of_curvature) + +# Arugments: +- `s::Vector{Float}` : vector of surface lengths. +- `controlpoint::Matrix{Float} : control point positions, axial in first row, radial in second row +- `ss::Float` : position along surface length to find radius of curvature. + +# Return: +- `radius_of_curvature::Float` : Radius of curvature at point `ss` along surface. +""" +function calculate_radius_of_curvature(s, controlpoint, ss) + + # spline x and y coordinates with respect to arc length + x_of_s = Akima_smooth(s, controlpoint[1, :]) + y_of_s = Akima_smooth(s, controlpoint[2, :]) + + #get first and second derivatives of coordinates with respect to arc length at integration step locations + xdot = FLOWMath.derivative.(Ref(x_of_s), ss) + xdotsp = Akima_smooth(s, FLOWMath.derivative.(Ref(x_of_s), s)) + xddot = FLOWMath.derivative.(Ref(xdotsp), ss) + + ydot = FLOWMath.derivative.(Ref(y_of_s), ss) + ydotsp = Akima_smooth(s, FLOWMath.derivative.(Ref(y_of_s), s)) + yddot = FLOWMath.derivative.(Ref(ydotsp), ss) + + # assemble the numerator and denominator of the radius of curvature expression + num = (xdot .^ 2 .+ ydot .^ 2) .^ (1.5) + den = xdot .* yddot .- ydot .* xddot + + # convex positive + if abs(den) < eps() + return 0.0 + else + return num ./ den + end +end + +""" + setup_boundary_layer_functions_green( + s, + vtan_duct, + duct_control_points, + operating_point, + boundary_layer_options; + verbose=false + ) + +# Arguments: +- `s::Vector{Float}` : cumulative sum of panel lengths between control points in the given index range, starting from zero. +- `vtan_duct::Vector{Float}` : tangential velocity magnitudes for the entire duct +- `duct_control_points::Matrix{Float}` : Control point coordinates along the duct surface +- `operating_point::OperatingPoint` : OperatingPoint object +- `boundary_layer_options::BoundaryLayerOptions` : BoundaryLayerOptions object + +# Returns: +- `boundary_layer_parameters::NamedTuple` : Namped Tuple containing boundary layer solver parameters: + - `edge_velocity::FLOWMath.Akima` : spline of edge velocities relative to surface length + - `edge_mach::FLOWMath.Akima` : spline of edge Mach numbers relative to surface length + - `edge_acceleration::FLOWMath.Akima` : spline of edge acceleration (dUe/ds) relative to surface length + - `edge_density::FLOWMath.Akima` : spline of edge density relative to surface length + - `edge_viscosity::FLOWMath.Akima` : spline of edge viscosity relative to surface length + - `r_coords::FLOWMath.Akima` : spline radial coordaintes relative to surface length + - `radial_derivative::FLOWMath.Akima` : spline of dr/ds relative to surface length + - `radius_of_curvature::FLOWMath.Akima` : spline of radius of curvature relative to surface length +""" +function setup_boundary_layer_functions_green( + s, + vtan_duct, + duct_control_points, + operating_point, + boundary_layer_options; + verbose=false, +) + + # Edge Velocities + # note: the ss that get's passed in is the ss in the full surface length, so in practice this will be stagnation s ± boundary layer step s + edge_velocity = Akima_smooth(s, vtan_duct) + + # Edge Accelerations (dUe/ds) + edge_acceleration(ss) = FLOWMath.derivative.(Ref(edge_velocity), ss) + + # r's + r_coords = Akima_smooth(s, duct_control_points[2, :]) + + # dr/ds + radial_derivative(ss) = FLOWMath.derivative.(Ref(r_coords), ss) + + # radii of curvature + radius_of_curvature(ss) = calculate_radius_of_curvature(s, duct_control_points, ss) + + # local mach number + edge_mach(ss) = calculate_mach(edge_velocity(ss), operating_point.asound[]) + + # local density + Pe(ss) = static_pressure(operating_point.Ptot[], edge_mach(ss)) + edge_density(ss) = static_density(Pe(ss), operating_point.asound[]) + + # local viscosity + function Te(ss) + return convert_temperature_to_kelvin.( + Ref(operating_point.units), + static_temperature.(Ref(operating_point.Ttot[]), edge_mach(ss)), + ) + end + function edge_viscosity(ss) + return convert_viscosity.(Ref(operating_point.units), sutherlands_law(Te(ss))) + end + + return (; + edge_velocity, + edge_mach, + edge_acceleration, + edge_density, + edge_viscosity, + r_coords, + radial_derivative, + radius_of_curvature, + verbose, + ) +end + +#---------------------------------# +# state initilization functions # +#---------------------------------# + +""" + d2_init(s_init, Rex) + +Initialize momentum thickness state with flat plate Schlichting solution + +# Arguments: +- `x::Float` : length +- `Rex::Float` : Reynolds number at the given length + +# Returns: +- `d2::Float` : momentum thickness +""" +function d2_init(x, Rex) + return 0.036 * x / Rex^0.2 +end + +""" +Rename of `calculate_H12bar0` used for initialization to avoid confusion +""" +function H12bar_init(Cf0, M) + return calculate_H12bar0(Cf0, M) +end + +""" +Wrapper of `calculate_CEeq` used for initialization to avoid confusion (`lambda` set to 1) +""" +function CE_init(Ctaueq0, M, Cf0) + return calculate_CEeq(Ctaueq0, M, 1, Cf0) +end + +""" + initialize_turbulent_boundary_layer_states_green( + s_init, r_init, Ue, M, rhoe, mue; verbose=false + ) + +Initialize the states for the turbulent boundary layer solve. + +# Arguments: +- `s_init::Float` : surface length starting point +- `r_init::Float` : surface radial position at starting point +- `Ue::Float` : edge velocity at starting point +- `M::Float` : edge Mach at starting point +- `rhoe::Float` : edge density at starting point +- `mue::Float` : edge dynamic viscosity at starting point + +# Returns: +- `initial_states::Vector{Float}` : Initial values for the states: r*delta2, Hbar12, CE +- `Cf_init::Vector{Float}` : initial value for friction coefficietn Cf (used for determining separateion) +- `H12_init::Vector{Float}` : initial value for shape factor H12 (used in viscous drag calculation) +""" +function initialize_turbulent_boundary_layer_states_green( + s_init, r_init, Ue, M, rhoe, mue; verbose=false +) + verbose && println("INITIALIZATION") + + # println("s_init: ", s_init) + # println("r_init: ", r_init) + # println("Ue: ", Ue) + # println("M: ", M) + # println("rho: ", rhoe) + # println("mue: ", mue) + + # - Initialize States - # + + # initialize momentum thickness (d2) using flat plate schlichting model + Rex = calculate_Re(rhoe, Ue, s_init, mue) + verbose && printdebug("Rex: ", Rex) + d2 = d2_init(s_init, Rex) + verbose && printdebug("d2: ", d2) + + # initialize H12bar (compressible shape factor) using _0 equations + Red2 = calculate_Re(rhoe, Ue, d2, mue) + verbose && printdebug("Red2: ", Red2) + Cf0 = calculate_Cf0(Red2, M) + verbose && printdebug("Cf0: ", Cf0) + H12bar0 = H12bar_init(Cf0, M) + verbose && printdebug("H12bar0: ", H12bar0) + + # initialize the entrainment coefficient (CE) using equilibrium equations + H1 = calculate_H1(H12bar0) + verbose && printdebug("H1: ", H1) + H12 = calculate_H12(H12bar0, M) + verbose && printdebug("H12: ", H12) + Cf = calculate_Cf(H12bar0, H12bar0, Cf0) + verbose && printdebug("Cf: ", Cf) + d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar0, H12, Cf, M) + verbose && printdebug("d2dUedsUeeq0: ", d2dUedsUeeq0) + CEeq0 = calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) + verbose && printdebug("CEeq0: ", CEeq0) + Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0, M) + verbose && printdebug("Ctaueq0: ", Ctaueq0) + CE = CE_init(Ctaueq0, M, Cf0) + verbose && printdebug("CE: ", CE) + + # initial states + return [r_init * d2, H12bar0, CE], Cf, H12 +end + +#---------------------------------# +# Intermediate functions # +#---------------------------------# + +""" + Fc(M) = sqrt(1.0 + 0.2 * M^2) +""" +function Fc(M) + return sqrt(1.0 + 0.2 * M^2) +end + +""" + FR(M) = 1.0 + 0.056 * M^2 +""" +function FR(M) + return 1.0 + 0.056 * M^2 +end + +""" + calculate_Re(ρ, V, L, μ) + +Calculate Reynolds number +""" +function calculate_Re(rhoe, Ue, d2, mue) + return rhoe * Ue * d2 / mue +end + +""" + calculate_Cf0(Red2, M) + +Calculate Cf₀ value based on momentum thickness Reynolds number (Red2) and edge mach number (M). +""" +function calculate_Cf0(Red2, M) + + return (0.01013 / (log(10, FR(M) * Red2) - 1.02) - 0.00075) / Fc(M) +end + +""" + calculate_H12bar0(Cf0, M) +""" +function calculate_H12bar0(Cf0, M) + return 1.0 / (1.0 - 6.55 * sqrt(Cf0 / 2.0 * (1.0 + 0.04 * M^2))) +end + +""" + calculate_Cf(H12bar, H12bar0, Cf0) + +Calculate friction coefficient. +""" +function calculate_Cf(H12bar, H12bar0, Cf0) + return Cf0 * (0.9 / (H12bar / H12bar0 - 0.4) - 0.5) +end + +""" + calculate_H12(H12bar, M, Pr=1.0) +""" +function calculate_H12(H12bar, M, Pr=1.0) + return (H12bar + 1.0) * (1.0 + Pr^(1.0 / 3.0) * M^2 / 5) - 1.0 +end + +""" + calculate_Ctau(CE, Cf0, M) +""" +function calculate_Ctau(CE, Cf0, M) + return (0.024 * CE + 1.2 * CE^2 + 0.32 * Cf0) * (1.0 + 0.1 * M^2) +end + +""" + calculate_F(CE, Cf0) +""" +function calculate_F(CE, Cf0) + return (0.02 * CE + CE^2 + 0.8 * Cf0 / 3) / (0.01 + CE) +end + +""" + calculate_H1(H12bar) +""" +function calculate_H1(H12bar) + return 3.15 + 1.72 / (H12bar - 1.0) - 0.01 * (H12bar - 1.0)^2 +end + +""" + calculate_dH12bardH1(H12bar) +""" +function calculate_dH12bardH1(H12bar) + return -(H12bar - 1.0)^2 / (1.72 + 0.02 * (H12bar - 1.0)^3) +end + +""" + calculate_richardson_number(H12bar, d2, H12, H1, R) +""" +function calculate_richardson_number(H12bar, d2, H12, H1, R) + return 2.0 * d2 / (3.0 * R) * (H12 + H1) * (H1 / H12bar + 0.3) +end + +""" + calculate_beta(Ri; hardness=100.0) + +Sigmoind blended version of piecewise function: + | 7.0 if Ri > 0 +β = | + | 4.5 if Ri < 0 + +# Arguments: +- `Ri::float` : Richardson number + +# Keyword Arguments: +- `hardness::float` : hardness factor for sigmoid blend + +# Returns: +- `beta::float` : factor used in secondary influence from longitudinal curvature. +""" +function calculate_beta(Ri, hardness=100.0) + return FLOWMath.sigmoid_blend(4.5, 7.0, Ri, 0.0, hardness) +end + +""" + longitudinal_curvature_influence(M, Ri) +""" +function longitudinal_curvature_influence(M, Ri) + return 1 + calculate_beta(Ri) * (1.0 + M^2 / 5) * Ri +end + +""" + lateral_strain_influence(H12bar, d2, H12, H1, r, drds) +""" +function lateral_strain_influence(H12bar, d2, H12, H1, r, drds) + return 1.0 - 7.0 / 3.0 * (H1 / H12bar + 1.0) * (H12 + H1) * d2 * drds / r +end + +""" + dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) +""" +function dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) + return 1.0 + 7.0 / 3.0 * M^2 * (H1 / H12bar + 1.0) * (H12 + H1) * d2 * dUeds / Ue +end + +""" + calculate_lambda(args...) + +λ = *(args...) +""" +function calculate_lambda(args...) + return *(args...) +end + +""" + calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) +""" +function calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) + return 1.25 * (Cf / 2.0 - ((H12bar - 1.0) / (6.432 * H12bar))^2 / (1.0 + 0.04 * M^2)) / + H12 +end + +""" + calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) + +Note: applies smooth-max to makes sure this value stays non-negative +""" +function calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) + return H1 * (Cf / 2.0 - (H12 + 1.0) * d2dUedsUeeq0) +end + +""" + calculate_Ctaueq0(CEeq0, Cf0, M) +""" +function calculate_Ctaueq0(CEeq0, Cf0, M) + return (0.24 * CEeq0 + 1.2 * CEeq0^2 + 0.32 * Cf0) * (1.0 + 0.1 * M^2) +end + +""" + calculate_CEeq(Ctaueq0, M, lambda, Cf0) +""" +function calculate_CEeq(Ctaueq0, M, lambda, Cf0) + c = Ctaueq0 / ((1.0 + 0.1 * M^2) * lambda^2) - 0.32 * Cf0 + return sqrt(c / 1.2 + 0.0001) - 0.01 +end + +""" + calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) +""" +function calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) + return (Cf / 2.0 - CEeq / H1) / (H12 + 1.0) +end + +#---------------------------------# +# Residual functions # +#---------------------------------# + +""" + boundary_layer_residual_green(y, s, parameters) + +Out-of-place version of `boundary_layer_residual_green!` +""" +function boundary_layer_residual_green(y, s, parameters) + dy = similar(y) .= 0 + return boundary_layer_residual_green!(dy, y, s, parameters) +end + +""" + boundary_layer_residual_green!(dy, y, s, parameters) + +Calculate dy give the current states, y, the input position, s, and various parameters. +""" +function boundary_layer_residual_green!(dy, y, s, parameters) + (; + r_coords, + edge_velocity, + edge_density, + edge_viscosity, + edge_mach, + edge_acceleration, + radius_of_curvature, + radial_derivative, + verbose, + ) = parameters + + # - rename for convenience - # + verbose && println(" Inputs:") + r = r_coords(s) + verbose && printdebug("r: ", r) + Ue = edge_velocity(s) + verbose && printdebug("Ue: ", Ue) + rhoe = edge_density(s) + verbose && printdebug("rhoe: ", rhoe) + mue = edge_viscosity(s) + verbose && printdebug("mue: ", mue) + M = edge_mach(s) + verbose && printdebug("M: ", M) + dUeds = edge_acceleration(s) + verbose && printdebug("dUeds: ", dUeds) + R = radius_of_curvature(s) + verbose && printdebug("R: ", R) + drds = radial_derivative(s) + verbose && printdebug("drds: ", drds) + + # - unpack variables - # + rd2, H12bar, CE = y + verbose && println(" States:") + verbose && printdebug("rd2: ", rd2) + verbose && printdebug("H12bar: ", H12bar) + verbose && printdebug("CE: ", CE) + + # - calculate intermediate variables - # + verbose && println(" Intermediates:") + d2 = rd2 / r + verbose && printdebug("d2: ", d2) + Red2 = calculate_Re(rhoe, Ue, d2, mue) + verbose && printdebug("Red2: ", Red2) + + Cf0 = calculate_Cf0(Red2, M) + verbose && printdebug("Cf0: ", Cf0) + + H12bar0 = calculate_H12bar0(Cf0, M) + verbose && printdebug("H12bar0: ", H12bar0) + + Cf = calculate_Cf(H12bar, H12bar0, Cf0) + verbose && printdebug("Cf: ", Cf) + + H12 = calculate_H12(H12bar, M) + verbose && printdebug("H12: ", H12) + + dH12bardH1 = calculate_dH12bardH1(H12bar) #yes, this should be negative according to example + verbose && printdebug("dH12bardH1: ", dH12bardH1) + H1 = calculate_H1(H12bar) + verbose && printdebug("H1: ", H1) + + F = calculate_F(CE, Cf0) + verbose && printdebug("F: ", F) + + d2dUedsUeeq0 = calculate_d2dUedsUeeq0(H12bar, H12, Cf, M) + verbose && printdebug("d2dUedsUeeq0: ", d2dUedsUeeq0) + + CEeq0 = calculate_CEeq0(H1, H12, Cf, d2dUedsUeeq0) + verbose && printdebug("CEeq0: ", CEeq0) + + Ctaueq0 = calculate_Ctaueq0(CEeq0, Cf0, M) + verbose && printdebug("Ctaueq0: ", Ctaueq0) + + Ctau = calculate_Ctau(CE, Cf0, M) + verbose && printdebug("Ctau: ", Ctau) + + if parameters.lambda + if parameters.longitudinal_curvature + Ri = calculate_richardson_number(H12bar, d2, H12, H1, r) + l1 = longitudinal_curvature_influence(M, Ri) + else + l1 = 1 + end + + if parameters.lateral_strain + l2 = lateral_strain_influence(H12bar, d2, H12, H1, r, drds) + else + l2 = 1 + end + + if parameters.dilation + l3 = dilation_influence(H12bar, d2, H12, H1, M, Ue, dUeds) + else + l3 = 1 + end + + lambda = calculate_lambda(l1, l2, l3) + else + lambda = 1 + end + verbose && printdebug("lambda: ", lambda) + + CEeq = calculate_CEeq(Ctaueq0, M, lambda, Cf0) + verbose && printdebug("CEeq: ", CEeq) + + d2dUedsUeeq = calculate_d2dUedsUeeq(H1, H12, Cf, CEeq) + verbose && printdebug("d2dUedsUeeq: ", d2dUedsUeeq) + + # - system of equations - # + + # momentum + dy[1] = r * Cf / 2.0 - (H12 + 2.0 - M^2) * rd2 * dUeds / Ue + + # entrainment + dy[2] = dH12bardH1 * (CE - H1 * (Cf / 2.0 - (H12 + 1.0) * d2 * dUeds / Ue)) / d2 + + # lag + dy[3] = + F / d2 * ( + 2.8 / (H12 + H1) * (sqrt(Ctaueq0) - lambda * sqrt(Ctau)) + d2dUedsUeeq - + d2 * dUeds / Ue * (1.0 + 0.075 * M^2 * ((1.0 + 0.2 * M^2) / (1.0 + 0.1 * M^2))) + ) + verbose && println(" Residuals:") + verbose && printdebug("d(rd2)/ds: ", dy[1]) + verbose && printdebug("d(H12bar)/ds: ", dy[2]) + verbose && printdebug("d(CE)/ds: ", dy[3]) + + return dy, [Cf; H12] +end + +""" + solve_green_boundary_layer!(f, rk, initial_states, steps, parameters; verbose=false) + +Integrate the turbulent boundary layer using a Runge-Kutta method. + +# Arguments: +- `f::function_handle` : Governing residual equations to integrate +- `rk::function_handle` : Runge-Kutta method to use (RK2 or RK4) +- `initial_states::Float` : initial states +- `steps::Vector{Float}` : steps for integration +- `parameters::NamedTuple` : boundary layer solve options and other parameters +""" +function solve_green_boundary_layer!( + f, rk, initial_states, steps, parameters; verbose=false +) + + # Unpack States and variables for viscous drag + u0, Cf0, H12_0 = initial_states + + # Initilization separate flags and outputs + sep = [false] + sepid = [1] + + # Allocate intermediate states and outputs + # TODO; put these in a cache that gets updated in place. + us = zeros(eltype(u0), length(u0), length(steps)) + Cfs = zeros(eltype(u0), length(steps)) + H12s = zeros(eltype(u0), length(steps)) + + # Initialize intermediate states and outputs + us[:, 1] = u0 + Cfs[1] = Cf0 + H12s[1] = H12_0 + + # Take Runge-Kutta Steps until either the end of the surface or separation is detected + for i in 1:(length(steps) - 1) + if verbose + println() + println(i, " of $(length(steps)-1)") + println(" s = ", steps[i]) + println(" ds = ", abs(steps[i + 1] - steps[i])) + println() + end + + # take step + us[:, i + 1], aux = rk( + f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters + ) + + Cfs[i + 1], H12s[i + 1] = aux + + sepid[1] = i + 1 + if Cfs[i + 1] <= 0.0 + sep[1] = true + break + end + end + + if sep[1] == true + + # - Interpolate to find actual s_sep - # + + usep = [ + FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[1, (sepid[] - 1):sepid[]], 0.0) + FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[2, (sepid[] - 1):sepid[]], 0.0) + FLOWMath.linear(Cfs[(sepid[] - 1):sepid[]], us[3, (sepid[] - 1):sepid[]], 0.0) + ] + + H12sep = FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], H12s[(sepid[] - 1):sepid[]], 0.0 + ) + + s_sep = FLOWMath.linear( + Cfs[(sepid[] - 1):sepid[]], steps[(sepid[] - 1):sepid[]], 0.0 + ) + else + usep = us[:, end] + H12sep = H12s[end] + s_sep = steps[end] + end + + # return states at separate, and separation shape factor, and surface length at separation + return usep, H12sep, s_sep, sepid +end diff --git a/src/postprocess/boundary_layer_head.jl b/src/postprocess/boundary_layer_head.jl index 19810d2d..e3c0ea3c 100644 --- a/src/postprocess/boundary_layer_head.jl +++ b/src/postprocess/boundary_layer_head.jl @@ -25,7 +25,7 @@ function setup_boundary_layer_functions_head( s, vtan_duct, - # duct_control_points, + duct_control_points, operating_point, boundary_layer_options; verbose=false, @@ -75,6 +75,15 @@ function calculate_H(H1) return FLOWMath.sigmoid_blend(hlt, hgeq, H1, 5.3) end +""" + limH1(H1) + +Returns a limited H1 to avoid undefined behavior +""" +function limH1(H1) + return FLOWMath.ksmax([H1; 3.3 + 1e-4]) +end + """ calculate_cf(H, Red2) @@ -85,24 +94,24 @@ function calculate_cf(H, Red2) end """ - boundary_layer_residual_head(y, s, parameters) + boundary_layer_residual_head(y, parameters, s) Out of place residual function for Head's method. """ -function boundary_layer_residual_head(y, s, parameters) +function boundary_layer_residual_head(y, parameters, s) dy = similar(y) .= 0 - return boundary_layer_residual_head!(dy, y, s, parameters) + return boundary_layer_residual_head!(dy, y, parameters, s) end """ - boundary_layer_residual_head!(dy, y, s, parameters) + boundary_layer_residual_head!(dy, y, parameters, s) Calculate dy give the current states, y, the input position, s, and various parameters. """ -function boundary_layer_residual_head!(dy, y, s, parameters) +function boundary_layer_residual_head!(dy, y, parameters, s) # - unpack parameters - # - (; verbose) = parameters + (; edge_velocity, edge_acceleration, edge_density, edge_viscosity, verbose) = parameters # - unpack variables - # d2, H1 = y @@ -110,12 +119,9 @@ function boundary_layer_residual_head!(dy, y, s, parameters) verbose && printdebug("H1: ", H1) # limit H1 to be greater than 3.3 - H1lim = FLOWMath.ksmax([H1; 3.3 + 1e-2]) + H1lim = limH1(H1) verbose && printdebug("H1lim: ", H1lim) - # - unpack variables - # - (; edge_velocity, edge_acceleration, edge_density, edge_viscosity) = parameters - # - Intermediate Calculations - # # determine dUedx dUedx = edge_acceleration(s) @@ -146,22 +152,45 @@ function boundary_layer_residual_head!(dy, y, s, parameters) dy[2] = 0.0306 / d2 * (H1lim - 3.0)^(-0.6169) - dUedx * H1lim / Ue - dy[1] * H1lim / d2 verbose && printdebug("dy[2]: ", dy[2]) - return dy, H + return dy +end + +function initialize_head_states(boundary_layer_functions, s_init; verbose=false) + + (; edge_density, edge_velocity, edge_viscosity) = boundary_layer_functions + + # - Initialize Boundary Layer States - # + H10 = 10.6 + d20 = + 0.036 * s_init / calculate_Re( + edge_density(s_init), + edge_velocity(s_init), + s_init, + edge_viscosity(s_init), + )^0.2 + + initial_states = [d20; H10] + H0 = 1.28 + + return initial_states, H0 end """ - solve_head_boundary_layer!(f, rk, initial_states, steps, parameters; verbose=false) + solve_head_boundary_layer!(f, ode, initial_states, steps, parameters; verbose=false) Integrate the turbulent boundary layer using a Runge-Kutta method. # Arguments: - `f::function_handle` : Governing residual equations to integrate -- `rk::function_handle` : Runge-Kutta method to use (RK2 or RK4) +- `ode::function_handle` : ODE method to use (RK2 or RK4) - `initial_states::Float` : initial states - `steps::Vector{Float}` : steps for integration - `parameters::NamedTuple` : boundary layer solve options and other parameters """ -function solve_head_boundary_layer!(f, rk, initial_states, steps, parameters; verbose=false) +function solve_head_boundary_layer!( + ::RK, ode, initial_states, steps, parameters; verbose=false +) + f = boundary_layer_residual_head # Unpack States and variables for viscous drag u0, H0 = initial_states @@ -190,11 +219,11 @@ function solve_head_boundary_layer!(f, rk, initial_states, steps, parameters; ve end # take step - us[:, i + 1], Hs[i + 1] = rk( - f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters - ) + us[:, i + 1] = ode(f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters) + + Hs[i + 1] = calculate_H(limH1(us[2, i + 1])) - sepid[1] = i + 1 + sepid[1] = i if Hs[i + 1] >= parameters.separation_criteria sep[1] = true break @@ -229,8 +258,10 @@ function solve_head_boundary_layer!(f, rk, initial_states, steps, parameters; ve steps[(sepid[] - 1):sepid[]], parameters.separation_criteria, ) + stepsol = steps[1:sepid[]] usol = us[:, 1:sepid[]] + else usep = us[:, end] Hsep = Hs[end] @@ -240,5 +271,56 @@ function solve_head_boundary_layer!(f, rk, initial_states, steps, parameters; ve end # return states at separate, and separation shape factor, and surface length at separation - return usep, Hsep, s_sep, sepid, usol, stepsol + return usep, Hsep, s_sep, usol, stepsol +end + +""" + solve_head_boundary_layer!(::DiffEq, ode, initial_states, steps, parameters; verbose=false) + +Integrate the turbulent boundary layer using a Runge-Kutta method. + +# Arguments: +- `f::function_handle` : Governing residual equations to integrate +- `ode::function_handle` : ODE method to use (one of the DifferentialEquations.jl options) +- `initial_states::Float` : initial states +- `steps::Vector{Float}` : steps for integration +- `parameters::NamedTuple` : boundary layer solve options and other parameters +""" +function solve_head_boundary_layer!( + ::DiffEq, ode, initial_states, steps, parameters; verbose=false +) + f = boundary_layer_residual_head! + + prob = ODEProblem(f, initial_states[1], [steps[1], steps[end]], parameters) + + # set up separation termination conditions + function condition(u, t, integrator) + return calculate_H(limH1(u[2])) - parameters.separation_criteria + end + + function affect!(integrator) + return terminate!(integrator) + end + + cb = ContinuousCallback(condition, affect!) + + sol = diffeq_solve( + prob, + ode(); + callback=cb, + abstol=eps(), + dtmax=1e-4, + # alg_hints=[:stiff], + # dt=parameters.first_step_size, + verbose=verbose, + ) + + usep = sol.u[end] + Hsep = calculate_H(usep[2]) + s_sep = sol.t[end] + usol = reduce(hcat, (sol.u)) + stepsol = sol.t + + # return states at separate, and separation shape factor, and surface length at separation + return usep, Hsep, s_sep, usol, stepsol end diff --git a/src/postprocess/boundary_layer_utils.jl b/src/postprocess/boundary_layer_utils.jl index 67512628..a4f091bb 100644 --- a/src/postprocess/boundary_layer_utils.jl +++ b/src/postprocess/boundary_layer_utils.jl @@ -39,59 +39,83 @@ Split the duct body surface at the leading edge of the duct. - `s_upper::Vector{Float}` : cumulative sum of upper side (nacelle) panel lengths - `s_lower::Vector{Float}` : cumulative sum of lower side (casing) panel lengths """ -function split_at_stagnation_point(duct_panel_lengths, duct_panel_tangents, Vtot_duct) - - # initialize stagnation point finder - stag_ids = [1, 1] - dp = zeros(eltype(Vtot_duct), 2) - dp[1] = dot(duct_panel_tangents[:, 1], Vtot_duct[:, 1]) - - # loop through panels and stop when dot product of panel vector and velocity vector changes sign - # TODO: have to start further than the size of the first step size for the ODE solver - for i in 3:length(duct_panel_lengths) - dp[2] = dot(duct_panel_tangents[:, i], Vtot_duct[:, i]) - - if sign(dp[1]) != sign(dp[2]) - stag_ids[:] .= [i - 1; i] - break +function split_at_stagnation_point( + duct_panel_lengths, duct_panel_tangents, Vtot_duct, Vtan_duct, first_step_size +) + offset = findfirst(x -> x > 1.1 * first_step_size, cumsum(duct_panel_lengths[end:-1:1])) + + s_tot = arc_lengths_from_panel_lengths(duct_panel_lengths[1:(end - offset)]) + + dots = [ + dot(duct_panel_tangents[:, i], Vtot_duct[:, i]) for + i in 1:(length(duct_panel_lengths) - offset) + ] + + if allequal(sign.(dots)) + #= + We're likely in a hover-ish case and there is no stagnation point. + Want to find the minimum Vtan in this case, since either the stagnation point is the trailing edge, or it has lifted off the surface and the most recent point of stagnation could be used. + =# + + _, minvtid = findmin(Vtan_duct) + vtsp = Akima_smooth(s_tot, Vtan_duct) + + if minvtid == length(Vtan_duct) + s_upper = nothing + s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end:-1:1]) + split_ratio = 1.0 + stag_ids = ones(Int, 2) .* length(duct_panel_lengths) + + return s_upper, s_lower, stag_ids, split_ratio, dots + else + stag_point = Roots.find_zero( + x -> FLOWMath.derivative(vtsp, x), + (s_tot[max(minvtid - 1, 1)], s_tot[minvtid + 1]), + Roots.Brent(); + atol=eps(), + ) end - dp[1] = dp[2] - end - if dp[1] == dp[2] - # we're likely in a hover-ish case and there is no stagnation point. + # elseif mindotid == length(duct_panel_lengths) - offset - 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) - - elseif stag_ids[2] == length(duct_panel_lengths) - # things flipped on the last panel - - s_upper = nothing - s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[end-1:-1:1]) - split_ratio = 1.0 - stag_ids .= length(duct_panel_lengths)-1 + # s_upper = nothing + # s_lower = arc_lengths_from_panel_lengths(duct_panel_lengths[(end - 1):-1:1]) + # split_ratio = 1.0 + # stag_ids = ones(Int,2) .* (length(duct_panel_lengths) - offset) else + mindotid = findfirst(x -> sign(x) != sign(dots[1]), dots) + dotsp = Akima_smooth(s_tot, dots) + + # println("mindotid: ", mindotid) + # println("mindot: ", dots[mindotid]) + # println("bracket: ", (s_tot[max(mindotid - 1, 1)], s_tot[mindotid])) + # println( + # "bracket vals: ", dotsp([(s_tot[max(mindotid - 1, 1)], s_tot[mindotid])...]) + # ) + # println("dots: ", dots) + + stag_point = Roots.find_zero( + dotsp, (s_tot[max(mindotid - 1, 1)], s_tot[mindotid]), Roots.Brent(); atol=eps() + ) + 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) + stag_ids = searchsortedfirst(s_tot, stag_point) .+ [-1, 0] + split_ratio = stag_point / s_tot[end] - partial_panel_lengths = [stag_interp, sum_length - stag_interp] + partial_panel_lengths = [ + stag_point - s_tot[stag_ids[1]] + s_tot[stag_ids[2]] - stag_point + ] - s_upper = arc_lengths_from_panel_lengths( - [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]] - ) + 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 = s_lower[end]/(s_lower[end]+s_upper[end]) - end + s_lower = arc_lengths_from_panel_lengths( + [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]] + ) - return s_upper, s_lower, stag_ids, split_ratio + return s_upper, s_lower, stag_ids, split_ratio, dots end """ diff --git a/src/postprocess/postprocess.jl b/src/postprocess/postprocess.jl index 72b615a9..439b6d45 100644 --- a/src/postprocess/postprocess.jl +++ b/src/postprocess/postprocess.jl @@ -145,7 +145,9 @@ Post-process a converged nonlinear solve solution. - `vz_wake` - `vr_wake` - `Cm_avg` - +- `reference_values` + - `Vinf` + - `Vref` """ function post_process( solver_options, @@ -513,19 +515,22 @@ function post_process( boundary_layer_options, Vtan_out[1:Int(body_vortex_panels.npanel[1])], Vtot_out[:, 1:Int(body_vortex_panels.npanel[1])], + body_vortex_panels.controlpoint[:, 1:Int(body_vortex_panels.npanel[1])], body_vortex_panels.influence_length[1:Int(body_vortex_panels.npanel[1])], 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=false, + # body_vortex_panels.node[2, Int(body_vortex_panels.nnode[1])], + Rref[1], + operating_point, + reference_parameters; + verbose=verbose, ) - body_thrust[1] -= duct_viscous_drag + # body_thrust[1] -= duct_viscous_drag else + duct_viscous_drag = [0.0, 0.0] boundary_layer_outputs = nothing end - ### --- TOTAL OUTPUTS --- ### # - Total Thrust - # @@ -571,8 +576,10 @@ function post_process( # panel strengths panel_strengths=gamb[1:(idmaps.body_totnodes)], # body thrust - total_thrust=sum(body_thrust), + body_force_coefficient=body_force_coefficient, + total_thrust=sum(body_thrust) - sum(duct_viscous_drag), thrust_comp=body_thrust, + duct_viscous_drag=duct_viscous_drag, induced_efficiency, # surface pressures cp_in, @@ -603,7 +610,7 @@ function post_process( vtan_centerbody_in, vtan_centerbody_out, # boundary layers - boundary_layers = boundary_layer_outputs, + boundary_layers=boundary_layer_outputs, ), # - Rotor Values - # rotors=(; diff --git a/src/postprocess/viscous_drag.jl b/src/postprocess/viscous_drag.jl index e30564cc..436c6214 100644 --- a/src/postprocess/viscous_drag.jl +++ b/src/postprocess/viscous_drag.jl @@ -24,38 +24,34 @@ function squire_young(d2, Ue, Uinf, H12) end """ - total_viscous_drag_duct(cd_upper, cd_lower, exit_radius, Vref, rhoinf) + total_viscous_drag_duct(cd_upper, cd_lower, rotor_tip_radius, Vref, rhoinf) Calculate the total viscous drag of the duct from Squire-Young drag coefficients, integrated about exit circumference. # Arguments: - `cdc_upper::Float` : upper side drag coefficient times refernce chord - `cdc_lower::Float` : lower side drag coefficient times refernce chord -- `exit_radius::Float` : radius used for integrating circumferentially +- `rotor_tip_radius::Float` : radius used for integrating circumferentially - `Vref::Float` : reference velocity (Vinf) - `rhoinf::Float` : freestream density # Returns: - `viscous_drag::Float` : viscous drag on duct """ -function total_viscous_drag_duct(cdc_upper, cdc_lower, exit_radius, Vref, rhoinf) +function total_viscous_drag_duct(cdc_upper, cdc_lower, rotor_tip_radius, Vref, rhoinf) # note: cdc's are already times chord, so no need to have a separate chord variable # drag per unit length dprime = 0.5 * rhoinf * Vref^2 * (cdc_upper + cdc_lower) # drag of annular airfoil - return dprime * 2.0 * pi * exit_radius + return dprime * 2.0 * pi * rotor_tip_radius end -#---------------------------------# -# HEAD'S METHOD # -#---------------------------------# - """ - compute_single_side_drag_coefficient_head( + compute_single_side_drag_coefficient( steps, - exit_radius, + rotor_tip_radius, operating_point, boundary_layer_options; verbose=false, @@ -65,67 +61,107 @@ Solve integral boundary layer and obtain viscous drag coefficient from Squire-Yo # Arguments: - `steps::Vector{Float}` : positions along surface for integration -- `exit_radius::Float` : radius at duct trailing edge (casing side) -- `operating_point::Float` : OperatingPoint object -- `boundary_layer_functions::NamedTuple` : Various Akima splines and other functions for boundary layer values +- `rotor_tip_radius::Float` : radius at duct trailing edge (casing side) +- `Vref::Float` : reference velocity +- `setup_boundary_layer_functions::NamedTuple` : Various Akima splines and other functions for boundary layer values - `boundary_layer_options::BoundaryLayerOptions` : BoundaryLayerOptions object # Returns: - `cd::Float` : viscous drag coefficient """ -function compute_single_side_drag_coefficient_head( +function compute_single_side_drag_coefficient( + boundary_layer_options::HeadsBoundaryLayerOptions, steps, - exit_radius, - operating_point, + rotor_tip_radius, + Vref, boundary_layer_functions, - boundary_layer_options; + separation_options; verbose=false, ) - (; edge_density, edge_velocity, edge_viscosity) = boundary_layer_functions - - # verbose && printdebug("initial velocity: ", edge_velocity(steps[1])) - # verbose && printdebug("initial density: ", edge_density(steps[1])) - # verbose && printdebug("initial viscosity: ", edge_viscosity(steps[1])) - - # - Initialize Boundary Layer States - # - H10 = 10.6 - d20 = - 0.036 * steps[1] / calculate_Re( - edge_density(steps[1]), - edge_velocity(steps[1]), - steps[1], - edge_viscosity(steps[1]), - ) - initial_states = [d20; H10] - H0 = 1.28 + return compute_single_side_drag_coefficient( + initialize_head_states, + solve_head_boundary_layer!, + steps, + rotor_tip_radius, + Vref, + boundary_layer_functions, + (; + boundary_layer_options.solver_type, + boundary_layer_options.ode, + boundary_layer_options.separation_criteria, + boundary_layer_options.first_step_size, + separation_options..., + ); + verbose=verbose, + ) +end - usep, Hsep, s_sep, sepid, usol, stepsol = solve_head_boundary_layer!( - boundary_layer_residual_head, - boundary_layer_options.rk, - [initial_states, H0], +function compute_single_side_drag_coefficient( + boundary_layer_options::GreensBoundaryLayerOptions, + steps, + rotor_tip_radius, + Vref, + boundary_layer_functions, + separation_options; + verbose=false, +) + return compute_single_side_drag_coefficient( + initialize_green_states, + solve_green_boundary_layer!, steps, - (; boundary_layer_functions..., boundary_layer_options.separation_criteria); - verbose=false, + rotor_tip_radius, + Vref, + boundary_layer_functions, + (; + boundary_layer_options.solver_type, + boundary_layer_options.ode, + boundary_layer_options.separation_criteria, + boundary_layer_options.first_step_size, + boundary_layer_options.lambda, + boundary_layer_options.longitudinal_curvature, + boundary_layer_options.lateral_strain, + boundary_layer_options.dilation, + separation_options..., + ); + verbose=verbose, ) +end - # verbose && println(s_sep / steps[end]) +function compute_single_side_drag_coefficient( + initialize_states, + solve_boundary_layer, + steps, + rotor_tip_radius, + Vref, + boundary_layer_functions, + single_side_boundary_layer_options; + verbose=false, +) + usep, Hsep, s_sep, usol, stepsol = solve_boundary_layer( + single_side_boundary_layer_options.solver_type, + single_side_boundary_layer_options.ode, + initialize_states(boundary_layer_functions, steps[1]; verbose=verbose), + steps, + (; boundary_layer_functions..., single_side_boundary_layer_options...); + verbose=verbose, + ) cdsq = squire_young( - usep[1], boundary_layer_functions.edge_velocity(s_sep), operating_point.Vinf[], Hsep + usep[1], boundary_layer_functions.edge_velocity(s_sep), Vref[], Hsep ) - cd = FLOWMath.ksmin([boundary_layer_options.separation_penalty; cdsq]) + cd = FLOWMath.ksmin([single_side_boundary_layer_options.separation_penalty; cdsq]) cdadd = FLOWMath.ksmax( [ 0.0 FLOWMath.linear( - [0.0; steps[end - boundary_layer_options.separation_allowance]], - [boundary_layer_options.separation_penalty; 0.0], + [0.0; steps[end - single_side_boundary_layer_options.separation_allowance]], + [single_side_boundary_layer_options.separation_penalty; 0.0], s_sep, ) ], - 100, + 1e8, ) cd += cdadd @@ -139,8 +175,9 @@ end Vtan_duct, cp_duct, duct_panel_lengths, - exit_radius, + rotor_tip_radius, operating_point, + reference_parameters, ) Determine total, dimensional viscous drag on the duct. @@ -149,8 +186,9 @@ Determine total, dimensional viscous drag on the duct. - `boundary_layer_options::BoundaryLayerOptions` : BoundaryLayerOptions object, used for dispatch as well - `Vtan_duct::Vector{Float}` : tangential velocity magnitudes for the entire duct - `duct_panel_lengths::Vector{Float}` : panel lengths for the entire duct -- `exit_radius::Float` : radius at duct trailing edge (casing side) -- `operating_point::Float` : OperatingPoint object +- `rotor_tip_radius::Float` : radius at duct trailing edge (casing side) +- `operating_point::OperatingPoint` : OperatingPoint object +- `reference_parameters::ReferenceParameters` : ReferenceParameters object # Returns: - `duct_viscous_drag::Float` : total viscous drag of duct @@ -165,28 +203,91 @@ Determine total, dimensional viscous drag on the duct. - `split_ratio` - `separation_point_ratio_upper` - `separation_point_ratio_lower` + - `cdc_upper` + - `cdc_lower` """ function compute_viscous_drag_duct( boundary_layer_options::HeadsBoundaryLayerOptions, Vtan_duct, Vtot_duct, + duct_control_points, duct_panel_lengths, duct_panel_tangents, - exit_radius, - operating_point; + rotor_tip_radius, + operating_point, + reference_parameters; verbose=false, ) + return compute_viscous_drag_duct( + setup_boundary_layer_functions_head, + boundary_layer_options, + Vtan_duct, + Vtot_duct, + duct_control_points, + duct_panel_lengths, + duct_panel_tangents, + rotor_tip_radius, + operating_point, + reference_parameters; + verbose=verbose, + ) +end + +function compute_viscous_drag_duct( + boundary_layer_options::GreensBoundaryLayerOptions, + Vtan_duct, + Vtot_duct, + duct_control_points, + duct_panel_lengths, + duct_panel_tangents, + rotor_tip_radius, + operating_point, + reference_parameters; + verbose=false, +) + return compute_viscous_drag_duct( + setup_boundary_layer_functions_green, + boundary_layer_options, + Vtan_duct, + Vtot_duct, + duct_control_points, + duct_panel_lengths, + duct_panel_tangents, + rotor_tip_radius, + operating_point, + reference_parameters; + verbose=verbose, + ) +end +function compute_viscous_drag_duct( + setup_boundary_layer_functions, + boundary_layer_options, + Vtan_duct, + Vtot_duct, + duct_control_points, + duct_panel_lengths, + duct_panel_tangents, + rotor_tip_radius, + operating_point, + reference_parameters; + verbose=false, +) # find stagnation point - s_upper, s_lower, stag_ids, split_ratio = split_at_stagnation_point( - duct_panel_lengths, duct_panel_tangents, Vtot_duct + s_upper, s_lower, stag_ids, split_ratio, dots = split_at_stagnation_point( + duct_panel_lengths, + duct_panel_tangents, + Vtot_duct, + Vtan_duct, + boundary_layer_options.first_step_size, ) # set up boundary layer solve parameters if !isnothing(s_upper) - boundary_layer_functions_upper = setup_boundary_layer_functions_head( + boundary_layer_functions_upper = setup_boundary_layer_functions( s_upper, [0.0; Vtan_duct[stag_ids[2]:end]], + duct_control_points, operating_point, boundary_layer_options; verbose=verbose, @@ -194,9 +295,10 @@ function compute_viscous_drag_duct( end # set up boundary layer solve parameters - boundary_layer_functions_lower = setup_boundary_layer_functions_head( + boundary_layer_functions_lower = setup_boundary_layer_functions( s_lower, [0.0; Vtan_duct[stag_ids[1]:-1:1]], + duct_control_points, operating_point, boundary_layer_options; verbose=verbose, @@ -205,33 +307,66 @@ function compute_viscous_drag_duct( # - Set integration steps - # 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 + if boundary_layer_options.solver_type == RK() + if isnothing(boundary_layer_options.upper_step_size) + 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 + else + upper_steps = + range( + 0.0, + s_upper[end] - boundary_layer_options.offset; + step=boundary_layer_options.upper_step_size, + ) .+ boundary_layer_options.offset + end + else + upper_steps = + range( + 0.0, + s_upper[end] - boundary_layer_options.offset; + length=length(s_upper), + ) .+ boundary_layer_options.offset + end end # lower side - lower_steps = - set_boundary_layer_steps( - boundary_layer_options.n_steps, - boundary_layer_options.first_step_size, - s_lower[end] - boundary_layer_options.offset, - ) .+ boundary_layer_options.offset + if boundary_layer_options.solver_type == RK() + if isnothing(boundary_layer_options.lower_step_size) + lower_steps = + set_boundary_layer_steps( + boundary_layer_options.n_steps, + boundary_layer_options.first_step_size, + s_lower[end] - boundary_layer_options.offset, + ) .+ boundary_layer_options.offset + else + lower_steps = + range( + 0.0, + s_lower[end] - boundary_layer_options.offset; + step=boundary_layer_options.lower_step_size, + ) .+ boundary_layer_options.offset + end + else + lower_steps = + range( + 0.0, s_lower[end] - boundary_layer_options.offset; length=length(s_lower) + ) .+ boundary_layer_options.offset + end # - Get drag coeffients - # if !isnothing(s_upper) - cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head( + cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient( + boundary_layer_options, upper_steps, - exit_radius, - operating_point, + rotor_tip_radius, + reference_parameters.Vref, 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, ); @@ -244,14 +379,13 @@ function compute_viscous_drag_duct( s_sep_upper = -1.0 end - cdc_lower, usol_lower, stepsol_lower, s_sep_lower = compute_single_side_drag_coefficient_head( + cdc_lower, usol_lower, stepsol_lower, s_sep_lower = compute_single_side_drag_coefficient( + boundary_layer_options, lower_steps, - exit_radius, - operating_point, + rotor_tip_radius, + reference_parameters.Vref, boundary_layer_functions_lower, (; - boundary_layer_options.rk, - boundary_layer_options.separation_criteria, separation_allowance=boundary_layer_options.separation_allowance_lower, separation_penalty=boundary_layer_options.separation_penalty_lower, ); @@ -260,7 +394,11 @@ function compute_viscous_drag_duct( # Calculate total viscous drag total_drag = total_viscous_drag_duct( - cdc_upper, cdc_lower, exit_radius, operating_point.Vinf[], operating_point.rhoinf[] + cdc_upper, + cdc_lower, + rotor_tip_radius, + reference_parameters.Vref[], + operating_point.rhoinf[], ) # Return @@ -276,6 +414,9 @@ function compute_viscous_drag_duct( split_ratio, separation_point_ratio_upper=s_sep_upper, separation_point_ratio_lower=s_sep_lower, + cdc_upper=cdc_upper, + cdc_lower=cdc_lower, + vtdotpv=dots, ) end diff --git a/src/preprocess/geometry/body_geometry.jl b/src/preprocess/geometry/body_geometry.jl index 80c9a1ab..60b5f91c 100644 --- a/src/preprocess/geometry/body_geometry.jl +++ b/src/preprocess/geometry/body_geometry.jl @@ -5,8 +5,9 @@ duct_coordinates, centerbody_coordinates, zwake, + duct_le_coordinates, ncenterbody_inlet, - nduct_inlet; + nduct_inlet, finterp=FLOWMath.akima, ) @@ -17,9 +18,10 @@ Reinterpolate duct and centerbody coordinates in order to make them compatible w - `rp_centerbody_coordinates::Matrix{Float}` : the re-paneled centerbody coordinates - `duct_coordinates::Matrix{Float}` : the input duct coordinates - `centerbody_coordinates::Matrix{Float}` : the input centerbody coordinates -- `zwake::Matrix{Float}` : the wake sheet panel node axial positions -- `ncenterbody_inlet::Matrix{Float}` : the number of panels to use for the centerbody inlet -- `nduct_inlet::Matrix{Float}` : the number of panels to use for the duct inlet +- `zwake::Vector{Float}` : the wake sheet panel node axial positions +- `duct_le_coordinates::Matrix{Float}` : [z r] coordinates of duct leading edge +- `ncenterbody_inlet::Int` : the number of panels to use for the centerbody inlet +- `nduct_inlet,::Int` : the number of panels to use for the duct inlet # Keyword Arguments - `finterp::Function=FLOWMath.akima` : interpolation method @@ -30,6 +32,7 @@ function reinterpolate_bodies!( duct_coordinates, centerbody_coordinates, zwake, + duct_le_coordinates, ncenterbody_inlet, nduct_inlet; finterp=FLOWMath.akima, @@ -37,61 +40,96 @@ function reinterpolate_bodies!( # - separate inner and outer duct coordinates - # dle, dleidx = findmin(view(duct_coordinates, :, 1)) - duct_inner_coordinates = view(duct_coordinates, dleidx:-1:1, :) - duct_outer_coordinates = view(duct_coordinates, dleidx:size(duct_coordinates, 1), :) + duct_inlet_length = zwake[1] - duct_le_coordinates[1] + + if duct_le_coordinates[2] > duct_coordinates[dleidx, 2] + duct_inner_coordinates = [ + duct_le_coordinates + view(duct_coordinates, dleidx:-1:1, :) + ] + duct_outer_coordinates = [ + duct_le_coordinates + view(duct_coordinates, (dleidx + 1):size(duct_coordinates, 1), :) + ] + elseif duct_le_coordinates[2] < duct_coordinates[dleidx, 2] + duct_inner_coordinates = [ + duct_le_coordinates + view(duct_coordinates, (dleidx - 1):-1:1, :) + ] + duct_outer_coordinates = [ + duct_le_coordinates + view(duct_coordinates, dleidx:size(duct_coordinates, 1), :) + ] + else + # coordinate is the zero + duct_inner_coordinates = view(duct_coordinates, dleidx:-1:1, :) + duct_outer_coordinates = view(duct_coordinates, dleidx:size(duct_coordinates, 1), :) + end # - interpolate inner duct geometry to provided grid locations - # - zduct_inner = view(duct_inner_coordinates, :, 1) - rduct_inner = view(duct_inner_coordinates, :, 2) - duct_trailing_index = min(length(zwake), searchsortedfirst(zwake, zduct_inner[end])) - new_zduct_grid = @view(zwake[1:duct_trailing_index]) - new_rduct_grid = finterp(zduct_inner, rduct_inner, new_zduct_grid) - # - interpolate outer duct geometry to provided grid locations - # - zduct_outer = view(duct_outer_coordinates, :, 1) - rduct_outer = view(duct_outer_coordinates, :, 2) + # rename for convenience + casing_z = view(duct_inner_coordinates, :, 1) + casing_r = view(duct_inner_coordinates, :, 2) - new_rduct_grid_outer = finterp(zduct_outer, rduct_outer, new_zduct_grid) + # index of casing trailing edge in wake discretization + casing_te_id = min(length(zwake), searchsortedfirst(zwake, casing_z[end])) - # - interpolate centerbody geometry to provided grid locations - # - zcenterbody = view(centerbody_coordinates, :, 1) - rcenterbody = view(centerbody_coordinates, :, 2) - centerbody_trailing_index = min( - length(zwake), searchsortedfirst(zwake, zcenterbody[end]) + # new points for casing aft of rotor + casing_in_wake_z = @view(zwake[1:casing_te_id]) + casing_in_wake_r = finterp(casing_z, casing_r, casing_in_wake_z) + + # repaneled casing inlet nodes + casing_inlet_z = scaled_cosine_spacing( + nduct_inlet + 1, 2 * duct_inlet_length, duct_le_coordinates[1]; mypi=pi / 2.0 ) - new_zcenterbody_grid = @view(zwake[1:centerbody_trailing_index]) - new_rcenterbody_grid = finterp(zcenterbody, rcenterbody, new_zcenterbody_grid) + casing_inlet_r = finterp(casing_z, casing_r, casing_inlet_z) + + # - interpolate outer duct geometry to provided grid locations - # - scale = new_zcenterbody_grid[1] - zcenterbody[1] - transform = zcenterbody[1] - new_zcenterbody = scaled_cosine_spacing( - ncenterbody_inlet + 1, 2 * scale, transform; mypi=pi / 2 + # rename for convenience + nacelle_z = view(duct_outer_coordinates, :, 1) + nacelle_r = view(duct_outer_coordinates, :, 2) + + nacelle_in_wake_z = collect( + range(zwake[1], nacelle_z[end]; length=length(casing_in_wake_z)) ) - new_rcenterbody = finterp(zcenterbody, rcenterbody, new_zcenterbody) + nacelle_in_wake_r = finterp(nacelle_z, nacelle_r, nacelle_in_wake_z) - scale = new_zduct_grid[1] - duct_coordinates[dleidx, 1] - transform = duct_coordinates[dleidx, 1] - new_zduct_inner = scaled_cosine_spacing( - nduct_inlet + 1, 2 * scale, transform; mypi=pi / 2 + # repaneled nacelle inlet nodes + nacelle_inlet_z = scaled_cosine_spacing( + nduct_inlet + 1, 2 * duct_inlet_length, duct_le_coordinates[1]; mypi=pi / 2.0 ) - new_rduct_inner = finterp(zduct_inner, rduct_inner, new_zduct_inner) + nacelle_inlet_r = finterp(nacelle_z, nacelle_r, nacelle_inlet_z) + + # - interpolate centerbody geometry to provided grid locations - # - # update outer duct geometry - # in order to make sure that the trailing edge panels are (nearly) the same length it's easiest just to use the same x-spacing for the outer side of the duct as well. - new_rduct_outer = finterp(zduct_outer, rduct_outer, new_zduct_inner) + # rename for convenience + centerbody_z = view(centerbody_coordinates, :, 1) + centerbody_r = view(centerbody_coordinates, :, 2) + + centerbody_te_id = min(length(zwake), searchsortedfirst(zwake, centerbody_z[end])) + centerbody_in_wake_z = @view(zwake[1:centerbody_te_id]) + centerbody_in_wake_r = finterp(centerbody_z, centerbody_r, centerbody_in_wake_z) + + centerbody_inlet_length = centerbody_in_wake_z[1] - centerbody_z[1] + centerbody_inlet_z = scaled_cosine_spacing( + ncenterbody_inlet + 1, 2 * centerbody_inlet_length, centerbody_z[1]; mypi=pi / 2.0 + ) + centerbody_inlet_r = finterp(centerbody_z, centerbody_r, centerbody_inlet_z) # assemble new duct coordinates rp_duct_coordinates .= hcat( - [reverse(new_zduct_grid)'; reverse(new_rduct_grid)'], - [reverse(new_zduct_inner)[2:end]'; reverse(new_rduct_inner)[2:end]'], - [new_zduct_inner[2:(end - 1)]'; new_rduct_outer[2:(end - 1)]'], - [new_zduct_grid'; new_rduct_grid_outer'], + [reverse(casing_in_wake_z)'; reverse(casing_in_wake_r)'], + [reverse(casing_inlet_z)[2:end]'; reverse(casing_inlet_r)[2:end]'], + [nacelle_inlet_z[2:end]'; nacelle_inlet_r[2:end]'], + [nacelle_in_wake_z[2:end]'; nacelle_in_wake_r[2:end]'], ) # assemble new centerbody coordinates rp_centerbody_coordinates .= hcat( - [new_zcenterbody[1:(end - 1)]'; new_rcenterbody[1:(end - 1)]'], - [new_zcenterbody_grid'; new_rcenterbody_grid'], + [centerbody_inlet_z[1:(end - 1)]'; centerbody_inlet_r[1:(end - 1)]'], + [centerbody_in_wake_z'; centerbody_in_wake_r'], ) # check that the splining didn't put any of the center body radial coordinates in the negative. diff --git a/src/preprocess/geometry/wake_geometry.jl b/src/preprocess/geometry/wake_geometry.jl index 97b0ff75..707f61dd 100644 --- a/src/preprocess/geometry/wake_geometry.jl +++ b/src/preprocess/geometry/wake_geometry.jl @@ -25,15 +25,35 @@ function discretize_wake( wake_length, npanels, dte_minus_cbte; + le_bracket=5, + fitscale=1e2, ) - # extract duct leading and trailing edge locations - duct_lez = minimum(duct_coordinates[:, 1]) - duct_tez = maximum(duct_coordinates[:, 1]) + # try a root finder to get the true leading edge + duct_lez, duct_minz = findmin(@view(duct_coordinates[:, 1])) + + duct_lesp = FLOWMath.Akima( + @view(duct_coordinates[(duct_minz - le_bracket):(duct_minz + le_bracket), 2]), + @view(duct_coordinates[(duct_minz - le_bracket):(duct_minz + le_bracket), 1]), + ) + + duct_ler = Roots.find_zero( + x -> FLOWMath.derivative(duct_lesp, x) * fitscale, + ( + duct_coordinates[duct_minz - le_bracket, 2], + duct_coordinates[duct_minz + le_bracket, 2], + ), + Roots.Brent(); + atol=eps(), + ) + duct_lez = duct_lesp(duct_ler) + + # trailing edge shouldn't need anything fancy + duct_tez = duct_coordinates[1, 1] # extract hub leading and trailing edge location - cb_lez = minimum(centerbody_coordinates[:, 1]) - cb_tez = maximum(centerbody_coordinates[:, 1]) + cb_lez = centerbody_coordinates[1, 1] + cb_tez = centerbody_coordinates[end, 1] # calculate duct chord duct_chord = max(duct_tez, cb_tez) - min(duct_lez, cb_lez) @@ -79,7 +99,7 @@ function discretize_wake( ridx = findall(x -> x in rotorzloc, zwake) # return dimensionalized wake x-coordinates - return zwake, ridx + return zwake, ridx, [duct_lez duct_ler] end """ @@ -97,7 +117,6 @@ Initialize the wake grid. - `wake_grid::Array{Float,3}` : 3D Array of axial and radial wake_grid points """ function initialize_wake_grid(rp_duct_coordinates, rp_centerbody_coordinates, zwake, rwake) - TF = promote_type( eltype(rp_duct_coordinates), eltype(rp_centerbody_coordinates), @@ -198,7 +217,6 @@ function initialize_wake_grid!( return wake_grid end - """ generate_wake_grid( problem_dimensions, @@ -320,7 +338,6 @@ function generate_wake_grid!( return wake_grid end - """ relax_grid!( grid_solver_options::GridSolverOptionsType, @@ -372,7 +389,10 @@ function relax_grid!( grid_solver_options.converged[1] = false if verbose - println(tabchar^ntab * "Solving Elliptic Grid System using $(grid_solver_options.algorithm) Method") + println( + tabchar^ntab * + "Solving Elliptic Grid System using $(grid_solver_options.algorithm) Method", + ) end # # precondition diff --git a/src/preprocess/preprocess.jl b/src/preprocess/preprocess.jl index 0cd93924..b1b6b725 100644 --- a/src/preprocess/preprocess.jl +++ b/src/preprocess/preprocess.jl @@ -129,6 +129,7 @@ function reinterpolate_geometry!( finterp=(x,y,xp)->FLOWMath.akima(x,y,xp,2.0*eps(),eps()), verbose=false, silence_warnings=true, + le_bracket=1, ) ##### ----- Extract Tuples ----- ##### @@ -147,13 +148,14 @@ function reinterpolate_geometry!( # - Discretize Wake z-coordinates - # # also returns indices of rotor locations and duct and center body trailng edges in the wake - zwake, rotor_indices_in_wake[:] = discretize_wake( + zwake, rotor_indices_in_wake[:], duct_le_coordinates = discretize_wake( duct_coordinates, centerbody_coordinates, rotorzloc, # rotor axial locations wake_length, npanels, - dte_minus_cbte, + dte_minus_cbte; + le_bracket=le_bracket ) # - Re-interpolate Bodies - # @@ -163,6 +165,7 @@ function reinterpolate_geometry!( duct_coordinates, centerbody_coordinates, zwake, + duct_le_coordinates, ncenterbody_inlet, nduct_inlet; finterp=finterp, diff --git a/src/utilities/bookkeeping.jl b/src/utilities/bookkeeping.jl index 21919a4a..3650dfab 100644 --- a/src/utilities/bookkeeping.jl +++ b/src/utilities/bookkeeping.jl @@ -4,20 +4,20 @@ Struct containing dimensions of the problem used throughout the analysis. - `nrotor` : number of rotors -- nwn` : number of wake nodes -- nwp` : number of wake panels -- ncp` : number of casing panels -- ndn` : number of duct nodes -- ncbn` : number of centerbody nodes -- nbn` : number of body nodes -- nbp` : number of body panels -- nws` : number of wake sheets (also rotor nodes) -- nbe` : number of blade elements (also rotor panels) -- nwsn` : number of nodes in each wake sheet -- nwsp` : number of panels in each wake sheet -- ndwin` : number of duct-wake interfacing nodes -- ncbwin` : number of centerbody-wake interfacing nodes -- nbodies=2` : number of bodies (currently hardcoded to 2) +- `nwn` : number of wake nodes +- `nwp` : number of wake panels +- `ncp` : number of casing panels +- `ndn` : number of duct nodes +- `ncbn` : number of centerbody nodes +- `nbn` : number of body nodes +- `nbp` : number of body panels +- `nws` : number of wake sheets (also rotor nodes) +- `nbe` : number of blade elements (also rotor panels) +- `nwsn` : number of nodes in each wake sheet +- `nwsp` : number of panels in each wake sheet +- `ndwin` : number of duct-wake interfacing nodes +- `ncbwin` : number of centerbody-wake interfacing nodes +- `nbodies=2` : number of bodies (currently hardcoded to 2) """ @kwdef struct ProblemDimensions{TI} nrotor::TI # number of rotors @@ -52,8 +52,9 @@ Determine all relevant dimensions to the problem based either on the paneling_co function get_problem_dimensions(paneling_constants::PanelingConstants) # - Extract Paneling Constants - # - (; npanels, ncenterbody_inlet, nduct_inlet, wake_length, nwake_sheets, dte_minus_cbte) = - paneling_constants + (; + npanels, ncenterbody_inlet, nduct_inlet, wake_length, nwake_sheets, dte_minus_cbte + ) = paneling_constants # number of rotors is one less than the length of npanels if the duct and hub trailing edges line up, and is two less if they don't nrotor = iszero(dte_minus_cbte) ? length(npanels) - 1 : length(npanels) - 2 @@ -83,7 +84,7 @@ function get_problem_dimensions(paneling_constants::PanelingConstants) ncbp += npanels[end - 1] end - # duct panels are 2x the number of nacelle panels + # duct panels are casing + nacelle panels ndp = 2 * ncp # duct and center body nodes are 1 more than number of panels @@ -119,21 +120,21 @@ function get_problem_dimensions(paneling_constants::PanelingConstants) end return ProblemDimensions(; - nrotor, # number of rotors - nwn, # number of wake nodes - nwp, # number of wake panels - ncp, # number of casing panels - ndn, # number of duct nodes - ncbn, # number of centerbody nodes - nbn, # number of body nodes - nbp, # number of body panels - nws, # number of wake sheets (also rotor nodes) - nbe, # number of blade elements (also rotor panels) - nwsn, # number of nodes in each wake sheet - nwsp, # number of panels in each wake sheet - ndwin, # number of duct-wake interfacing nodes - ncbwin, # number of centerbody-wake interfacing nodes - nbodies=2, #hard code this for now. + nrotor, # number of rotors + nwn, # number of wake nodes + nwp, # number of wake panels + ncp, # number of casing panels + ndn, # number of duct nodes + ncbn, # number of centerbody nodes + nbn, # number of body nodes + nbp, # number of body panels + nws, # number of wake sheets (also rotor nodes) + nbe, # number of blade elements (also rotor panels) + nwsn, # number of nodes in each wake sheet + nwsp, # number of panels in each wake sheet + ndwin, # number of duct-wake interfacing nodes + ncbwin, # number of centerbody-wake interfacing nodes + nbodies=2, # hard code this for now. ) end diff --git a/src/utilities/inputs.jl b/src/utilities/inputs.jl index 96a919a4..f3173fda 100644 --- a/src/utilities/inputs.jl +++ b/src/utilities/inputs.jl @@ -203,7 +203,7 @@ Note that unlike other input structures, this one, in general, does not define f # Arguments -- `nduct_inlet::Int` : The number of panels to use for the duct inlet (this number is used for both the casing and nacelle re-paneling) +- `nduct_inlet::Int` : The number of panels to use for the casing inlet. - `ncenterbody_inlet::Int` : The number of panels to use for the centerbody inlet. - `npanels::AbstractVector{Int}` : A vector containing the number of panels between discrete locations inside the wake. Specifically, the number of panels between the rotors, between the last rotor and the first body trailing edge, between the body trailing edges (if different), and between the last body trailing edge and the end of the wake. The length of this vector should be N+1 (where N is the number of rotors) if the duct and centerbody trailing edges are aligned, and N+2 if not. - `dte_minus_cbte::Float` : An indicator concerning the hub and duct trailing edge relative locations. Should be set to -1 if the duct trailing edge axial position minus the centerbody trailing edge axial position is negative, +1 if positive (though any positive or negative number will suffice), and zero if the trailing edges are aligned. diff --git a/src/utilities/misc.jl b/src/utilities/misc.jl index ceda89c4..3e491bdd 100644 --- a/src/utilities/misc.jl +++ b/src/utilities/misc.jl @@ -38,16 +38,15 @@ Formatted printing when you have lots of debugging to do and want things to line """ function printdebug(variable_name, variable, nspaces=4; colw=16) if length(variable) == 1 - s = @sprintf "%*s %-*s %2.18f" nspaces " " colw variable_name variable[1] + s = @sprintf "%*s %-*s %2.18f" nspaces " " colw variable_name ForwardDiff.value(variable[1]) println(rstrip(s, ['0', ' '])) else @printf "%*s %-*s\n" nspaces " " colw variable_name for v in variable - s = @sprintf "%*s %-*s %2.18f" (nspaces) " " colw "" v + s = @sprintf "%*s %-*s %2.18f" (nspaces) " " colw "" ForwardDiff.value(v) println(rstrip(s, ['0', ' '])) end end - # println("\t$(variable_name)" * v) return nothing end diff --git a/src/utilities/ode_solvers.jl b/src/utilities/ode_solvers.jl index 5b7c25d4..5ea98dcd 100644 --- a/src/utilities/ode_solvers.jl +++ b/src/utilities/ode_solvers.jl @@ -1,3 +1,6 @@ +#---------------------------------# +# Runge-Kutta # +#---------------------------------# """ RK2(f, y, s, ds, parameters) @@ -12,11 +15,11 @@ """ function RK2(f, y, s, ds, parameters) parameters.verbose && println(" 1st call:") - k1, _ = f(y, s, parameters) + k1 = f(y, parameters, s) parameters.verbose && println(" 2nd call:") - k2, aux = f(y .+ (ds / 2) .* k1, s + (ds / 2), parameters) + k2 = f(y .+ (ds / 2) .* k1, parameters, s + (ds / 2)) unext = @. y + k2 * ds - return unext, aux + return unext end """ @@ -32,12 +35,10 @@ end - `parameters::NamedTuple` : BoundaryLayerOptions and other various parameters """ function RK4(f, y, s, ds, parameters) - k1, aux1 = f(y, s, parameters) - k2, aux2 = f(y .+ (ds / 2) .* k1, s + (ds / 2), parameters) - k3, aux3 = f(y .+ (ds / 2) .* k2, s + (ds / 2), parameters) - k4, aux4 = f(y .+ ds .* k3, s + ds, parameters) - uaux = @. y + (k1 + k2 * 2 + k3 * 2 + k4) * ds / 6 - aux = [(aux1[i] + aux2[i] * 2 + aux3[i] * 2 + aux4[i]) / 6 for i in length(aux1)] - return unext, aux + k1 = f(y, parameters, s) + k2 = f(y .+ (ds / 2) .* k1, parameters, s + (ds / 2)) + k3 = f(y .+ (ds / 2) .* k2, parameters, s + (ds / 2)) + k4 = f(y .+ ds .* k3, parameters, s + ds) + unext = @. y + (k1 + k2 * 2 + k3 * 2 + k4) * ds / 6 + return unext end - diff --git a/src/utilities/options.jl b/src/utilities/options.jl index 7048355d..467b41bd 100644 --- a/src/utilities/options.jl +++ b/src/utilities/options.jl @@ -69,6 +69,14 @@ Used in boundary layer method dispatch """ abstract type BoundaryLayerOptions end +""" + abstract type RK + +Used for selecting boundary layer ODE solver +""" +struct RK end +struct DiffEq end + #---------------------------------# # QUADRATURE TYPES # #---------------------------------# @@ -663,19 +671,23 @@ end - `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) +- `solver_type::Type = DiffEq()` : type of ODE solver (RK() or DiffEq()) +- `ode::Function = RadauIIA5` : solver to use for boundary layer integration (RadauIIA5, RK4, or RK2 available) - `separation_criteria::Float=3.0` : value of H12 after which separation should happen. - `separation_allowance_upper::Int=10` : upper side allowance for how many steps ahead of the trailing edge we'll allow separation without penalty - `separation_allowance_lower::Int=10` : lower side allowance for how many steps ahead of the trailing edge we'll allow separation without penalty - `separation_penalty_upper::Float=0.2` : upper side maximum penalty value for separation (at leading edge) - `separation_penalty_lower::Float=0.2` : lower side maximum penalty value for separation (at leading edge) """ -@kwdef struct HeadsBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp,Ts} <: BoundaryLayerOptions +@kwdef struct HeadsBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp,Ts,Tsol,Tssl,Tssu} <: BoundaryLayerOptions model_drag::Tb=false n_steps::Ti = Int(5e2) first_step_size::Tf = 1e-6 + upper_step_size::Tssu=nothing + lower_step_size::Tssl=nothing offset::To = 1e-3 - rk::Tfun = RK2 + solver_type::Tsol = DiffEq() + ode::Tfun = RadauIIA5 separation_criteria::Ts=3.0 separation_allowance_upper::Ti=10 separation_allowance_lower::Ti=10 @@ -700,13 +712,14 @@ Known Bugs: - `n_steps::Int = Int(2e2)` : number of steps to use in boundary layer integration - `first_step_size::Float = 1e-3` : size of first step in boundary layer integration - `offset::Float = 1e-2` : size of offset for (where to initialize) boundary layer integration -- `rk::Function = RK4` : solver to use for boundary layer integration (RK4 or RK2 available) +- `solver_type::Type = DiffEq()` : type of ODE solver (RK() or DiffEq()) +- `ode::Function = RadauIIA5` : solver to use for boundary layer integration (RadauIIA5, RK4, or RK2 available) - `separation_allowance_upper::Int=3` : upper side allowance for how many steps ahead of the trailing edge we'll allow separation without penalty - `separation_allowance_lower::Int=3` : lower side allowance for how many steps ahead of the trailing edge we'll allow separation without penalty - `separation_penalty_upper::Float=0.2` : upper side maximum penalty value for separation (at leading edge) - `separation_penalty_lower::Float=0.2` : lower side maximum penalty value for separation (at leading edge) """ -@kwdef struct GreensBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp} <: BoundaryLayerOptions +@kwdef struct GreensBoundaryLayerOptions{Tb,Tf,Tfun,Ti,To,Tp,Ts,Tsol,Tssl,Tssu} <: BoundaryLayerOptions model_drag::Tb=true lambda::Tb = false longitudinal_curvature::Tb = true @@ -714,8 +727,12 @@ Known Bugs: dilation::Tb = true n_steps::Ti = Int(2e2) first_step_size::Tf = 1e-6 + upper_step_size::Tssu=nothing + lower_step_size::Tssl=nothing offset::To = 1e-3 - rk::Tfun = RK2 + solver_type::Tsol = DiffEq() + ode::Tfun = RadauIIA5 + separation_criteria::Ts=0.0 separation_allowance_upper::Ti=25 separation_allowance_lower::Ti=25 separation_penalty_upper::Tp=0.2 diff --git a/src/visualization/calculate_velocity_field.jl b/src/visualization/calculate_velocity_field.jl new file mode 100644 index 00000000..2474195b --- /dev/null +++ b/src/visualization/calculate_velocity_field.jl @@ -0,0 +1,108 @@ +function get_induced_velocities!( + velocities, + points, + v_bs, + v_ws, + v_rs, + body_vortex_panels, + body_vortex_strengths, + wake_vortex_panels, + wake_vortex_strengths, + rotor_source_panels, + rotor_source_strengths, + Vinf, + integration_options, +) + v_bs .= 0.0 + v_ws .= 0.0 + v_rs .= 0.0 + + # body-induced velocity + DuctAPE.induced_velocities_from_vortex_panels_on_points!( + v_bs, + points, + body_vortex_panels.node, + body_vortex_panels.nodemap, + body_vortex_panels.influence_length, + integration_options, + ) + + velocities[:, 1] .+= v_bs[:, :, 1] * body_vortex_strengths + velocities[:, 2] .+= v_bs[:, :, 2] * body_vortex_strengths + + # wake-induced velocity + DuctAPE.induced_velocities_from_vortex_panels_on_points!( + v_ws, + @view(points[:, :]), + wake_vortex_panels.node, + wake_vortex_panels.nodemap, + wake_vortex_panels.influence_length, + integration_options, + ) + + velocities[:, 1] .+= v_ws[:, :, 1] * wake_vortex_strengths + velocities[:, 2] .+= v_ws[:, :, 2] * wake_vortex_strengths + + # rotor-induced velocity + DuctAPE.induced_velocities_from_source_panels_on_points!( + v_rs, + @view(points[:, :]), + rotor_source_panels.node, + rotor_source_panels.nodemap, + rotor_source_panels.influence_length, + integration_options, + ) + + velocities[:, 1] .+= v_rs[:, :, 1] * rotor_source_strengths + velocities[:, 2] .+= v_rs[:, :, 2] * rotor_source_strengths + + return nothing +end + +function calculate_velocity_field( + field_points, + body_vortex_panels, + body_vortex_strengths, + wake_vortex_panels, + wake_vortex_strengths, + rotor_source_panels, + rotor_source_strengths, + Vinf; + integration_options=IntegrationOptions(), +) + + TF = promote_type( + eltype(body_vortex_panels.node), + eltype(wake_vortex_panels.node), + eltype(rotor_source_panels.node), + ) + + # Initialize + velocities = zeros(TF, length(field_points[1, :]), 2) + + # So we don't have to re-allocate all the memory: + v_bs = zeros(TF, length(field_points[1,:]), Int(body_vortex_panels.totnode[]), 2) + v_ws = zeros(TF, length(field_points[1,:]), Int(wake_vortex_panels.totnode[]), 2) + v_rs = zeros(TF, length(field_points[1,:]), Int(rotor_source_panels.totnode[]), 2) + + velocities[:, 1] .= Vinf + + # - Calculate Velocity at Current Point - # + get_induced_velocities!( + velocities, + field_points, + v_bs, + v_ws, + v_rs, + body_vortex_panels, + body_vortex_strengths, + wake_vortex_panels, + wake_vortex_strengths, + rotor_source_panels, + rotor_source_strengths, + Vinf, + integration_options, + ) + + return velocities[:, 1], velocities[:, 2] +end diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl index 6d0668f9..b4aab4de 100644 --- a/src/visualization/plot_recipes.jl +++ b/src/visualization/plot_recipes.jl @@ -18,6 +18,7 @@ struct plotMomentum end # Streamlines struct plotStreamlines end +# struct plotVelocityField end #---------------------------------# # PLOT Duct GEOMETRY # @@ -554,7 +555,7 @@ end quinary=RGB(155 / 255, 82 / 255, 162 / 255), #purple plotsgray=RGB(128 / 255, 128 / 255, 128 / 255), #gray ), - scale_thickness=50.0, + scale_thickness=5.0, bl_ylim=nothing, ) color_palette --> [ @@ -746,3 +747,65 @@ end return nothing end + +##---------------------------------# +## PLOT VELOCITY FIELD # +##---------------------------------# + +#@recipe function plot_velocity_field( +# ::plotVelocityField, +# points, +# bvp, +# bvs, +# wvp, +# wvs, +# rsp, +# rss, +# vinf; +# integration_options=IntegrationOptions(), +# default_colors=(; +# primary=RGB(1 / 255, 149 / 255, 226 / 255), #blue +# secondary=RGB(189 / 255, 10 / 255, 53 / 255), #red +# tertiary=RGB(76 / 255, 173 / 255, 59 / 255), #green +# quaternary=RGB(238 / 255, 167 / 255, 46 / 255), #orange +# quinary=RGB(155 / 255, 82 / 255, 162 / 255), #purple +# plotsgray=RGB(128 / 255, 128 / 255, 128 / 255), #gray +# ), +# scale_thickness=50.0, +#) +# color_palette --> [ +# default_colors.primary, +# default_colors.secondary, +# default_colors.tertiary, +# default_colors.quaternary, +# default_colors.quinary, +# default_colors.plotsgray, +# ] + +# # Plot specific values +# aspect_ratio --> 1 + +# u, v = calculate_velocity_field( +# points, +# bvp, +# bvs, +# wvp, +# wvs, +# rsp, +# rss, +# vinf; +# integration_options=integration_options, +# ) + +# for (z, r) in zip(eachcol(points[1, :]), eachcol(points[2, :])) +# @series begin +# color --> default_colors.plotsgray +# label --> "" +# linewidth --> 0.5 +# quiver --> (u, v) +# return z, r +# end +# end + +# return nothing +#end diff --git a/test/pre_processing_tests.jl b/test/pre_processing_tests.jl index f50ff20f..8897f90f 100644 --- a/test/pre_processing_tests.jl +++ b/test/pre_processing_tests.jl @@ -10,13 +10,14 @@ println("\nPRECOMPUTED ROTOR & WAKE INPUTS") # get input data include("data/basic_two_rotor_for_test_NEW.jl") - zwake, rotor_indices_in_wake = dt.discretize_wake( + zwake, rotor_indices_in_wake, duct_le_coordinates = dt.discretize_wake( duct_coordinates, centerbody_coordinates, rotor.rotorzloc, # rotor axial locations paneling_constants.wake_length, paneling_constants.npanels, paneling_constants.dte_minus_cbte; + le_bracket=1, ) @test zwake == [0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0] @@ -45,6 +46,7 @@ println("\nPRECOMPUTED ROTOR & WAKE INPUTS") duct_coordinates, centerbody_coordinates, zwake, + duct_le_coordinates, paneling_constants.ncenterbody_inlet, paneling_constants.nduct_inlet; finterp=fm.linear, @@ -55,10 +57,13 @@ println("\nPRECOMPUTED ROTOR & WAKE INPUTS") @test size(rp_duct_coordinates, 2) == 2 * size(duct_coordinates, 1) - 1 @test size(rp_centerbody_coordinates, 2) == 2 * size(centerbody_coordinates, 1) - 1 - @test rp_duct_coordinates == [ - 1.0 0.75 0.5 0.25 0.0 0.25 0.5 0.75 1.0 - 2.0 1.75 1.5 1.75 2.0 2.25 2.5 2.25 2.0 - ] + @test isapprox( + rp_duct_coordinates, + [ + 1.0 0.75 0.5 0.25 0.0 0.25 0.5 0.75 1.0 + 2.0 1.75 1.5 1.75 2.0 2.25 2.5 2.25 2.0 + ], + ) @test rp_centerbody_coordinates == [ 0.0 0.25 0.5 0.75 1.0 0.0 0.25 0.5 0.25 0.0 @@ -554,7 +559,7 @@ end # copy over operating point for f in fieldnames(typeof(operating_point)) if f != :units - solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) end end