diff --git a/Project.toml b/Project.toml index 9217b1b6..82a10455 100644 --- a/Project.toml +++ b/Project.toml @@ -1,23 +1,28 @@ name = "DuctAPE" uuid = "ad8e49fd-fab7-444e-af4a-0daba3b8bf11" authors = ["Judd Mehr"] -version = "0.5.1" +version = "0.6.0" [deps] BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" FLOWMath = "6cb5d3fb-0fe8-4cc2-bd89-9fe0b19a99d3" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" FixedPoint = "3325f569-5a18-4e7d-8356-246b69339eea" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" ImplicitAD = "e7cbb90b-9b31-4eb2-a8c8-45099c074ee1" +LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f" LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9" NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56" PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc" +RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" +Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4" SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" @@ -33,7 +38,6 @@ FastGaussQuadrature = "1" FiniteDiff = "2" FixedPoint = "1" ForwardDiff = "0.10" -ImplicitAD = "0.3" LineSearches = "7" MINPACK = "1" NLsolve = "4" diff --git a/docs/.gitignore b/docs/.gitignore index a303fff2..96db14da 100644 --- a/docs/.gitignore +++ b/docs/.gitignore @@ -1,2 +1,3 @@ build/ site/ +gr-temp/ diff --git a/docs/make.jl b/docs/make.jl index f55876cf..664cb24c 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -29,6 +29,7 @@ makedocs(; "Outputs" => "DuctAPE/advanced_usage/outputs.md", # "Manual Geometry" => "DuctAPE/advanced_usage/manual_repaneling.md", ], + "Visualization" => "DuctAPE/visualization.md", "API" => [ "Public API Reference" => "DuctAPE/api/public_api.md", "Private API Reference" => [ diff --git a/docs/src/DuctAPE/advanced_usage/option.md b/docs/src/DuctAPE/advanced_usage/option.md index 5705ef8b..9f5a11f4 100644 --- a/docs/src/DuctAPE/advanced_usage/option.md +++ b/docs/src/DuctAPE/advanced_usage/option.md @@ -145,6 +145,37 @@ DuctAPE.set_options(; !!! note "Iteration Counters" The `iterations` field (not to be confused with the `iterations_limit` field) in the solver options should generally not be changed. They automatically save (in-place) the number of iterations the solver performs and can be accessed after the analysis is run. +## Boundary Layer Solvers + +If desired, a one-way turbulent boundary layer can be modeled, from which an approximate viscous drag can be determined. +Currently, only the Head's method for turbulent boundary layer computation is working. +In the case of separation, a separation drag penalty is applied based on values selected in the options. + +```@docs; canonical=false +DuctAPE.HeadsBoundaryLayerOptions +``` + +Here is an example of a possible boundary layer option setting: + +```julia +# Define Boundary Layer Settings +boundary_layer_options = DuctAPE.HeadsBoundaryLayerOptions(; + model_drag=true, + separation_penalty_upper=0.1, + separation_penalty_lower=0.1, + separation_allowance_upper=3, + separation_allowance_lower=25, +) + +# set all the options +DuctAPE.set_options(; + integration_options=integration_options, + grid_solver_options=wake_solve_options, + solver_options=aero_solver_options, + boundary_layer_options=boundary_layer_options, +) +``` + ## Advanced Options for Multi-point analyses For using advanced options in multi-point analyses, there are various changes that need to be made to avoid run-time errors. diff --git a/docs/src/DuctAPE/api/private_postprocess.md b/docs/src/DuctAPE/api/private_postprocess.md index 395eb941..0806cb0a 100644 --- a/docs/src/DuctAPE/api/private_postprocess.md +++ b/docs/src/DuctAPE/api/private_postprocess.md @@ -1,3 +1,7 @@ +```@contents +Pages = ["private_postprocess.md"] +Depth = 5 +``` ```@docs DuctAPE.post_process @@ -51,4 +55,54 @@ DuctAPE.get_blade_loads DuctAPE.get_blade_loads! ``` +## Boundary Layer + +### Thermodynamics +```@docs +DuctAPE.sa1 +DuctAPE.sa2 +DuctAPE.standard_atmosphere +DuctAPE.ideal_gas_rho +DuctAPE.sutherlands_law +DuctAPE.speed_of_sound +DuctAPE.calculate_mach +DuctAPE.total_temperature +DuctAPE.total_pressure +DuctAPE.static_temperature +DuctAPE.static_pressure +DuctAPE.static_density +DuctAPE.convert_temperature_to_kelvin +DuctAPE.convert_viscosity +``` + +### General Boundary Layer Functions + +```@docs +DuctAPE.arc_lengths_from_panel_lengths +DuctAPE.split_at_stagnation_point +DuctAPE.bl_step_fun +DuctAPE.set_boundary_layer_steps +DuctAPE.RK2 +DuctAPE.RK4 +``` +### Head's Method Specific Functions + +```@docs +DuctAPE.setup_boundary_layer_functions_head +DuctAPE.calculate_H +DuctAPE.calculate_cf +DuctAPE.boundary_layer_residual_head +DuctAPE.boundary_layer_residual_head! +DuctAPE.solve_head_boundary_layer! +``` + +### Viscous Drag + +```@docs +DuctAPE.squire_young +DuctAPE.total_viscous_drag_duct +DuctAPE.compute_viscous_drag_duct +DuctAPE.compute_single_side_drag_coefficient_head +DuctAPE.compute_viscous_drag_duct_schlichting +``` diff --git a/docs/src/DuctAPE/api/private_preprocess.md b/docs/src/DuctAPE/api/private_preprocess.md index c4ccde46..4dca3a6f 100644 --- a/docs/src/DuctAPE/api/private_preprocess.md +++ b/docs/src/DuctAPE/api/private_preprocess.md @@ -1,4 +1,3 @@ - ```@contents Pages = ["private_preprocess.md"] Depth = 5 diff --git a/docs/src/DuctAPE/api/public_api.md b/docs/src/DuctAPE/api/public_api.md index aacee094..c5f0018d 100644 --- a/docs/src/DuctAPE/api/public_api.md +++ b/docs/src/DuctAPE/api/public_api.md @@ -69,6 +69,14 @@ DuctAPE.setup_analysis DuctAPE.analyze ``` +## Postprocess + +```@docs +DuctAPE.BoundaryLayerOptions +DuctAPE.HeadsBoundaryLayerOptions +DuctAPE.GreensBoundaryLayerOptions +``` + ## Miscellaneous ### Airfoil/Geometry Manipulation diff --git a/docs/src/DuctAPE/api/visualiztion_api.md b/docs/src/DuctAPE/api/visualiztion_api.md new file mode 100644 index 00000000..29efc9b9 --- /dev/null +++ b/docs/src/DuctAPE/api/visualiztion_api.md @@ -0,0 +1,3 @@ +```@docs +DuctAPE.generate_plots +``` diff --git a/docs/src/DuctAPE/theory_latex/build.sh b/docs/src/DuctAPE/theory_latex/build.sh old mode 100755 new mode 100644 diff --git a/docs/src/DuctAPE/theory_latex/clean.sh b/docs/src/DuctAPE/theory_latex/clean.sh old mode 100755 new mode 100644 diff --git a/docs/src/DuctAPE/tutorial.md b/docs/src/DuctAPE/tutorial.md index a343cce6..383560ee 100644 --- a/docs/src/DuctAPE/tutorial.md +++ b/docs/src/DuctAPE/tutorial.md @@ -346,7 +346,7 @@ DuctAPE.OperatingPoint ```@example tutorial # Freestream -Vinf = 0.0 # hover condition +Vinf = 30.0 rhoinf = 1.226 asound = 340.0 muinf = 1.78e-5 @@ -356,7 +356,7 @@ RPM = 8000.0 Omega = RPM * pi / 30 # if using RPM, be sure to convert to rad/s # utilizing the constructor function to put things in vector types -operating_point = DuctAPE.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) +operating_point = DuctAPE.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) nothing # hide ``` @@ -436,7 +436,7 @@ Running a multi-point analysis on the example geometry given there, it might loo ```@example tutorial # - Advance Ratio Range - # -Js = range(0.0, 2.0; step=0.01) +Js = range(0.1, 2.0; step=0.01) # - Calculate Vinfs - # D = 2.0 * rotor.Rtip[1] # rotor diameter @@ -523,8 +523,6 @@ nothing #hide And then we can plot the data to compare DFDC and DuctAPE. ```@example tutorial -using Plots - # set up efficiency plot pe = plot(; xlabel="Advance Ratio", ylabel="Efficiency") diff --git a/docs/src/DuctAPE/visualization.md b/docs/src/DuctAPE/visualization.md new file mode 100644 index 00000000..b8e96c16 --- /dev/null +++ b/docs/src/DuctAPE/visualization.md @@ -0,0 +1,84 @@ +# Visualization + +There are several convenience plotting methods implemented in DuctAPE based on [RecipesBase](https://juliaplots.org/RecipesBase.jl/stable/). +In addition a general function for plotting the suite of available plots or animations is provided in the `generate_plots` function. + +```@docs; canonical=false +DuctAPE.generate_plots +``` + +```@setup visualize +using DuctAPE +include("../assets/plots_default.jl") +gr() +``` + +The following generates animations across the given advance ratios. + +!!! warning "Plotting Streamlines" + Currently, plotting streamlines, especially animations, takes an exceptionally long time. + +```@julia +using Plots + +# - Advance Ratio Range - # +advance_ratios = range(0.1, 2.0; step=0.01) + +# - Calculate Vinfs - # +D = 2.0 * rotor.Rtip[1] # rotor diameter +n = RPM / 60.0 # rotation rate in revolutions per second +Vinfs = advance_ratios * n * D + +# - Set Operating Points - # +operating_points = [deepcopy(operating_point) for i in 1:length(Vinfs)] +for (iv, v) in enumerate(Vinfs) + operating_points[iv].Vinf[] = v +end + +# - Run Multi-point Analysis - # +outs, ins, success_flags = DuctAPE.analyze( + ducted_rotor, + operating_points, + reference_parameters, + DuctAPE.set_options( + operating_points; + boundary_layer_options=DuctAPE.HeadsBoundaryLayerOptions(; + model_drag=true, n_steps=1000, separation_criteria=3.0 + ), + ); + return_inputs=true, +) + +DuctAPE.generate_plots( + DuctAPE.animatedPlots(), + Plots, # Pass in the Plots namespace + ins, + outs; + plot_pressure=true, + plot_velocity=true, + plot_boundary_layer=true, + plot_streamlines=true, + save_path="../assets/", + static_file_type=".png", + (; + custom_defaults..., + # size=(600,400), # causes misalignment issues + cp_ylim=(-3, 3), # keyword argument to set ylim for cp plots + vtan_ylim=(0, 3), # keyword argument to set ylim for vtan plots + bl_ylim=(0.1, 0.25), # keyword argument to set ylim for boundary layer plots + )..., +) +nothing # hide +``` + +!!! note "Custom Defaults" + Additional arguments splatted into `generate_plots` are passed into `Plots.plot` directly as keyword arguments. In this case, `custom_defaults` happens to be the defaults associated with the plot formatting used in these docs. Any arguments passed in this way will override any options set in the plots recipes for all the plots. In the plots shown here, we have overridden the color palatte, but nothing else. + +!!! warning "Custom Sizes" + For some reason, specifying figure size causes misalignment with pressure and velocity distributions and the underlayed geometry. + +![](../assets/geometry.png) +![](../assets/surface_pressure.gif) +![](../assets/surface_velocity.gif) +![](../assets/boundary_layer.gif) +![](../assets/streamlines.gif) diff --git a/docs/src/assets/boundary_layer.gif b/docs/src/assets/boundary_layer.gif new file mode 100644 index 00000000..bbebda29 Binary files /dev/null and b/docs/src/assets/boundary_layer.gif differ diff --git a/docs/src/assets/define_propulsor.jl b/docs/src/assets/define_propulsor.jl index 307f783a..166ee84c 100644 --- a/docs/src/assets/define_propulsor.jl +++ b/docs/src/assets/define_propulsor.jl @@ -40,7 +40,7 @@ RPM = 8000.0 Omega = RPM * pi / 30 # if using RPM, be sure to convert to rad/s # utilizing the constructor function to put things in vector types -operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) +operating_point = dt.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) nduct_inlet = 50 ncenterbody_inlet = 30 diff --git a/docs/src/assets/gen_trim_logo.sh b/docs/src/assets/gen_trim_logo.sh old mode 100755 new mode 100644 diff --git a/docs/src/assets/geometry.png b/docs/src/assets/geometry.png new file mode 100644 index 00000000..edaedf1a Binary files /dev/null and b/docs/src/assets/geometry.png differ diff --git a/docs/src/assets/plots_default.jl b/docs/src/assets/plots_default.jl index 00f12c5a..20593747 100644 --- a/docs/src/assets/plots_default.jl +++ b/docs/src/assets/plots_default.jl @@ -2,7 +2,7 @@ Default Plots Settings for creating tikz native figures for dissertation =# -using Plots; +using Plots gr() using LaTeXStrings using Measures @@ -13,11 +13,6 @@ byublue = RGB(0.0, 46.0 / 255, 93.0 / 255) #BYU Navy Blue darkblue = RGB(0 / 255, 25 / 255, 50 / 255) byured = RGB(155.0 / 255, 0, 0) #"BYU" Red darkred = RGB(50 / 255, 0 / 255, 25 / 255) -# middlegray = RGB(128.0 / 255, 128.0 / 255, 128.0 / 255) #Middle Gray -# myblue = RGB(0.0 / 255, 92.0 / 255, 171.0 / 255) # royal blue -# myred = RGB(192.0 / 255, 83.0 / 255, 103.0 / 255) # royal red -# mygreen = RGB(143.0 / 255, 166.0 / 255, 81.0 / 255) # royal green -# mygray = RGB(130.0 / 255, 130.0 / 255, 130.0 / 255) # royal gray primary = RGB(0, 92 / 255, 171 / 255) #blue secondary = RGB(190 / 255, 76 / 255, 77 / 255) #red @@ -27,113 +22,14 @@ quinary = RGB(190 / 255, 147 / 255, 61 / 255) #yellow plotsgray = RGB(128 / 255, 128 / 255, 128 / 255) #gray function plots_default(; - # #:Plot - # background_color=RGBA(1, 1, 1, 0), background_color=nothing, - # background_color_outside = nothing, - # display_type, - dpi=300, - # extra_kwargs, - # extra_plot_kwargs, - fontfamily="Palatino_Roman", - # fontfamily="Palatino Roman", + fontfamily="cmunrm", foreground_color=RGB(128 / 255, 128 / 255, 128 / 255), #gray, plotsgray, - # html_output_format, - # inset_subplots, - # layout, - # link, - # overwrite_figure, - # plot_title, - # plot_title_location, - # plot_titlefontcolor, - # plot_titlefontfamily, - # plot_titlefonthalign, - # plot_titlefontrotation, - # plot_titlefontsize, - # plot_titlefontvalign, - # pos, - # show, - size=(400, 300) .* 5.0 ./ 4.0, #it appears that 100 ≈ 1inch in LaTeX - # size=(600, 450), #it appears that 100 ≈ 1inch in LaTeX - # size=(800, 600), #it appears that 100 ≈ 1inch in LaTeX - # tex_output_standalone, - # thickness_scaling, - # warn_on_unsupported, - # window_title, - - ####################### - # :Series - ####################### - # arrow, - # bar_edges, - # bar_position, - # bar_width, - # bins, - # colorbar_entry, - # connections, - # contour_labels, - # contours, - # extra_kwargs, - # fill_z, - fillalpha=0.125, - fillcolor=RGB(128 / 255, 128 / 255, 128 / 255), - # fillrange, - # group, - # hover, - # label, - # levels, - # line_z, - # linealpha, - # linecolor, - # linestyle, - linewidth=1.0, - # marker_z, - # markeralpha, - # markercolor, - # markershape, - markersize=1, - markerstrokealpha=0, - # markerstrokecolor, - # markerstrokestyle, - # markerstrokewidth, - # normalize, - # orientation, - # primary, - # quiver, - # ribbon, - # series_annotations, - # seriesalpha, - # seriescolor, - # seriestype, - # show_empty_bins, - # smooth, - # stride, - # subplot, - # weights, - # x, - # xerror, - # y, - # yerror, - # z, - # zerror - - ####################### - # :Subplot - ####################### - # annotationcolor, - annotationfontfamily="Palatino_Roman", - annotationfontsize=10, - # annotationhalign, - # annotationrotation, - # annotations, - # annotationvalign, - # aspect_ratio, + size=(400, 300), + annotationfontfamily="cmunrm", background_color_inside=nothing, background_color_legend=nothing, background_color_subplot=nothing, - # bottom_margin, - # camera, - # clims, color_palette=[ RGB(0, 92 / 255, 171 / 255), #blue RGB(190 / 255, 76 / 255, 77 / 255), @@ -142,336 +38,29 @@ function plots_default(; RGB(190 / 255, 147 / 255, 61 / 255), RGB(128 / 255, 128 / 255, 128 / 255), #gray ], - #color_palette=[ - # RGB(0.0, 46.0 / 255.0, 93.0 / 255.0), #BYU Blue - # RGB(155.0 / 255.0, 0.0, 0.0), #"BYU" Red - # RGB(128.0 / 255.0, 128.0 / 255.0, 128.0 / 255.0), #Middle Gray - # RGB(162.0 / 255.0, 227.0 / 255.0, 162.0 / 255.0), #Light Green - # RGB(243.0 / 255.0, 209.0 / 255.0, 243.0 / 255.0), #Pink - # RGB(205.0 / 255.0, 179.0 / 255.0, 0.0), #Yellow - # RGB(161.0 / 255.0, 161.0 / 255.0, 226.0 / 255.0), #Purple - #], - # colorbar, - # colorbar_continuous_values, - # colorbar_discrete_values, - # colorbar_fontfamily, - # colorbar_formatter, - # colorbar_scale, - # colorbar_tickfontcolor, - # colorbar_tickfontfamily, - # colorbar_tickfonthalign, - # colorbar_tickfontrotation, - # colorbar_tickfontsize, - # colorbar_tickfontvalign, - # colorbar_ticks, - # colorbar_title, - # colorbar_title_location, - # colorbar_titlefontcolor, - # colorbar_titlefontfamily, - # colorbar_titlefonthalign, - # colorbar_titlefontrotation, - # colorbar_titlefontsize, - # colorbar_titlefontvalign, - # extra_kwargs, - # fontfamily_subplot, foreground_color_legend=nothing, - # foreground_color_subplot, - # foreground_color_title, - # framestyle = :zerolines, - # left_margin, - # legend=false, # include legend true/false - # legendfontcolor, - # legendfontfamily, - # legendfonthalign, - # legendfontrotation, - # legendfontsize, - # legendfontvalign, - # legendtitle, - # legendtitlefontcolor, - # legendtitlefontfamily, - # legendtitlefonthalign, - # legendtitlefontrotation, - # legendtitlefontsize, - # legendtitlefontvalign, - margin=0mm, - # projection, - # right_margin, - # subplot_index, - # title, - # titlefontcolor, - # titlefontfamily, - # titlefonthalign, - # titlefontrotation, - titlefontsize=14, - # titlefontvalign, - # titlelocation, - # top_margin - - ##################### - # :Axis - ##################### - # discrete_values, - # draw_arrow, - # flip, - # foreground_color_axis, - # foreground_color_border, - # foreground_color_grid, - # foreground_color_guide, - # foreground_color_minor_grid, - # foreground_color_text, - # formatter, + margin=5mm, grid=false, # background grid true/false - # gridalpha, - # gridlinewidth=0.5, - # gridstyle, - # guide, - # guide_position, - # guidefontcolor, - # guidefontfamily, - # guidefonthalign, - # guidefontrotation, - guidefontsize=11, - # guidefontvalign, - # ylims=(0, 3), - # xlims=(0, 2), - # link, - # minorgrid, - # minorgridalpha, - # minorgridlinewidth, - # minorgridstyle, - # minorticks, - # mirror, - # rotation, - # scale, - # showaxis = false, #turns off spines and tick labels, but not ticks - # tick_direction, - # tickfontcolor, - # tickfontfamily, - # tickfonthalign, - # tickfontrotation, - tickfontsize=8, - # tickfontvalign, - # ticks=false, #turns off tick marks - # widen, + yguidefontrotation=-90, ) Plots.default() Plots.default(; - # #:Plot background_color=background_color, - # display_type, - dpi=dpi, - # extra_kwargs=extra_kwargs, - # extra_plot_kwargs=extra_plot_kwargs, fontfamily=fontfamily, foreground_color=foreground_color, - # html_output_format, - # inset_subplots, - # layout, - # link, - # overwrite_figure, - # plot_title, - # plot_title_location, - # plot_titlefontcolor, - # plot_titlefontfamily, - # plot_titlefonthalign, - # plot_titlefontrotation, - # plot_titlefontsize, - # plot_titlefontvalign, - # pos, - # show, size=size, - # tex_output_standalone, - # thickness_scaling, - # warn_on_unsupported, - # window_title, - - ####################### - # :Series - ####################### - # arrow, - # bar_edges, - # bar_position, - # bar_width, - # bins, - # colorbar_entry, - # connections, - # contour_labels, - # contours, - # extra_kwargs, - # fill_z, - # fillalpha=fillalpha, - # fillcolor=fillcolor, - # fillrange, - # group, - # hover, - # label, - # levels, - # line_z, - # linealpha, - # linecolor, - # linestyle, - linewidth=linewidth, - # marker_z, - # markeralpha, - # markercolor, - # markershape, - markersize=markersize, - markerstrokealpha=markerstrokealpha, - # markerstrokecolor, - # markerstrokestyle, - # markerstrokewidth, - # normalize, - # orientation, - # primary, - # quiver, - # ribbon, - # series_annotations, - # seriesalpha, - # seriescolor, - # seriestype, - # show_empty_bins, - # smooth, - # stride, - # subplot, - # weights, - # x, - # xerror, - # y, - # yerror, - # z, - # zerror - - ####################### - # :Subplot - ####################### - # annotationcolor, annotationfontfamily=annotationfontfamily, - annotationfontsize=annotationfontsize, - # annotationhalign, - # annotationrotation, - # annotations, - # annotationvalign, - # aspect_ratio, background_color_inside=background_color_inside, background_color_legend=background_color_legend, background_color_subplot=background_color_subplot, - # bottom_margin, - # camera, - # clims, color_palette=color_palette, - # colorbar, - # colorbar_continuous_values, - # colorbar_discrete_values, - # colorbar_fontfamily, - # colorbar_formatter, - # colorbar_scale, - # colorbar_tickfontcolor, - # colorbar_tickfontfamily, - # colorbar_tickfonthalign, - # colorbar_tickfontrotation, - # colorbar_tickfontsize, - # colorbar_tickfontvalign, - # colorbar_ticks, - # colorbar_title, - # colorbar_title_location, - # colorbar_titlefontcolor, - # colorbar_titlefontfamily, - # colorbar_titlefonthalign, - # colorbar_titlefontrotation, - # colorbar_titlefontsize, - # colorbar_titlefontvalign, - # extra_kwargs, - # fontfamily_subplot, foreground_color_legend=foreground_color_legend, - # foreground_color_subplot, - # foreground_color_title, - # framestyle = :zerolines, - # left_margin, - # legend=false, # include legend true/false - # legendfontcolor, - # legendfontfamily, - # legendfonthalign, - # legendfontrotation, - # legendfontsize, - # legendfontvalign, - # legendtitle, - # legendtitlefontcolor, - # legendtitlefontfamily, - # legendtitlefonthalign, - # legendtitlefontrotation, - # legendtitlefontsize, - # legendtitlefontvalign, margin=margin, - # projection, - # right_margin, - # subplot_index, - # title, - # titlefontcolor, - # titlefontfamily, - # titlefonthalign, - # titlefontrotation, - titlefontsize=titlefontsize, - # titlefontvalign, - # titlelocation, - # top_margin - - ##################### - # :Axis - ##################### - # discrete_values, - # draw_arrow, - # flip, - # foreground_color_axis, - # foreground_color_border, - # foreground_color_grid, - # foreground_color_guide, - # foreground_color_minor_grid, - # foreground_color_text, - # formatter, grid=grid, - # gridalpha, - # gridlinewidth=0.5, - # gridstyle, - # guide, - # guide_position, - # guidefontcolor, - # guidefontfamily, - # guidefonthalign, - # guidefontrotation, - guidefontsize=guidefontsize, - # guidefontvalign, - # ylims=(0, 3), - # xlims=(0, 2), - # link, - # minorgrid, - # minorgridalpha, - # minorgridlinewidth, - # minorgridstyle, - # minorticks, - # mirror, - # rotation, - # scale, - # showaxis = false, #turns off spines and tick labels, but not ticks - # tick_direction, - # tickfontcolor, - # tickfontfamily, - # tickfonthalign, - # tickfontrotation, - tickfontsize=tickfontsize, - # tickfontvalign, - # ticks=false, #turns off tick marks - # widen, + yguidefontrotation=yguidefontrotation, ) - return nothing + return (; color_palette) end -plots_default() - -function latexsize(x, y) - return (x * 96, y * 96) -end - -# - Rotate y-axis - # -extra_kwargs = Dict(:subplot => Dict("ylabel style" => "{align=center, rotate=-90}")) - +custom_defaults = plots_default() diff --git a/docs/src/assets/streamlines.gif b/docs/src/assets/streamlines.gif new file mode 100644 index 00000000..b286fcfe Binary files /dev/null and b/docs/src/assets/streamlines.gif differ diff --git a/docs/src/assets/surface_pressure.gif b/docs/src/assets/surface_pressure.gif new file mode 100644 index 00000000..42f3b07b Binary files /dev/null and b/docs/src/assets/surface_pressure.gif differ diff --git a/docs/src/assets/surface_velocity.gif b/docs/src/assets/surface_velocity.gif new file mode 100644 index 00000000..e03141a4 Binary files /dev/null and b/docs/src/assets/surface_velocity.gif differ diff --git a/src/C4Blade/dfdc_polar_fitting.jl b/src/C4Blade/dfdc_polar_fitting.jl new file mode 100644 index 00000000..6e337707 --- /dev/null +++ b/src/C4Blade/dfdc_polar_fitting.jl @@ -0,0 +1,125 @@ +project_dir = dirname(dirname(@__FILE__)) +if project_dir == "" + project_dir = "." +end + +# TODO: change datapath to test directory or validation data directory +datapath = project_dir * "/validation/data/" +savepath = project_dir * "/validation/data/" +# savepath = project_dir * "/dfdc_airfoil_approximations/" +afname = "dfdc_naca4412" + +using LsqFit +using FLOWMath +const fm = FLOWMath +using Plots +using DuctAPE.C4Blade +const c4b = C4Blade + +function linfit(aoas, cl0, dclda) + return dclda .* aoas .+ cl0 +end + +function quadfit(cls, cdmin, clcdmin, dcddcl2) + return dcddcl2 * (cls .- clcdmin) .^ 2 .+ cdmin +end + +# - Read in Polar - # +info, Re_ref, mach, aoa, cl, cd = c4b.parsefile(datapath * "naca4412.dat", true) +clsp = fm.Akima(aoa, cl) +cdsp = fm.Akima(aoa, cd) + +# - Plot Polar - # +plot(aoa * 180 / pi, cl) + +# - Chop Polar Down for Fits and replot - # +aoachop = range(-5, 10, 100) * pi / 180 +clchop = clsp(aoachop) +cdchop = cdsp(aoachop) +pcl = plot(aoachop, clchop) +pcd = plot(cdchop, clchop) + +# - Get cdmin and cl at cdmin as well as zero lift angle of attack- # +cdmin, cdminid = findmin(cdchop) +clcdmin = clchop[cdminid] + +_, zeroclid = findmin(abs.(clchop)) +alpha0 = fm.linear( + [clchop[zeroclid - 1], clchop[zeroclid]], + [aoachop[zeroclid - 1], aoachop[zeroclid]], + 0.0, +) + +_, zeroalphaid = findmin(abs.(aoachop)) +cl0 = fm.linear( + [aoachop[zeroalphaid - 1], aoachop[zeroalphaid]], + [clchop[zeroalphaid - 1], clchop[zeroalphaid]], + 0.0, +) + +# - Fit drag curve - # +p0 = [0.05] +quadfitwrap(cls, p) = quadfit(clchop, cdmin, clcdmin, p[1]) +cdfit = LsqFit.curve_fit(quadfitwrap, clchop, cdchop, p0) + +dcddcl2 = cdfit.param[1] + +# - Plot fitted curve - # +plot!(pcd, quadfit(clchop, cdmin, clcdmin, cdfit.param[1]), clchop) + +# - Fit Lift Curve - # +p0 = [2.0 * pi] +linfitwrap(aoas, p) = linfit(aoas, cl0, p[1]) +clfit = LsqFit.curve_fit(linfitwrap, aoachop, clchop, p0) + +dclda = clfit.param[1] + +# - Plot fitted curve - # +plot!(pcl, aoachop, linfit(aoachop, cl0, clfit.param[1])) + +# - Select clmin, clmax - # +clmin, clminid = findmin(clchop) +clmaxid = findlast(x -> x < 0.1, clfit.resid) +clmax = clchop[clmaxid] + +# - Manually set various paramters - # +Re_exp = -0.5 +mcrit = 0.7 +cmcon = 0.0 +dclda_stall = 0.1 +dcl_stall = 0.1 + +# - `alpha0::Float` : zero lift angle of attack +# - `clmax::Float` : maximum cl +# - `clmin::Float` : minimum cl +# - `dclda::Float` : lift curve slope (1/radians) +# - `dclda_stall::Float` : lift curve slope post-stall (1/radians) +# - `dcl_stall::Float` : cl increment from initial to total stall. +# - `cdmin::Float` : minimum cd +# - `cldmin::Float` : cl at cdmin +# - `dcdcl2::Float` : quadratic curve factor for cd curve (d(cd)/d(cl^2)) +# - `cmcon::Float` : pitching moment constant +# - `Re_ref::Float` : reference Reynolds number at which cd values apply +# - `Re_exp::Float` : Reynolds number exponent scaling (cd = cd*(Re/Re_ref)^Re_exp) +# - `mcrit::Float` : critical Mach number + +f = open(savepath * afname * ".jl", "w") + +write( + f, + "alpha0 = $alpha0 + clmax = $clmax + clmin = $clmin + dclda = $dclda + dclda_stall = $dclda_stall + dcl_stall = $dcl_stall + cdmin = $cdmin + clcdmin = $clcdmin + dcddcl2 = $dcddcl2 + cmcon = $cmcon + Re_ref = $Re_ref + Re_exp = $Re_exp + mcrit = $mcrit", +) + +close(f) diff --git a/src/DuctAPE.jl b/src/DuctAPE.jl index 1778631f..9dc42366 100644 --- a/src/DuctAPE.jl +++ b/src/DuctAPE.jl @@ -42,9 +42,18 @@ using NLsolve #Includes Anderson Solver using LineSearches # used in newton solver using ForwardDiff # used for jacobian for newton solver +# For boundary layer stuff +using Roots + # - Utility Packages - # using FLOWMath # used for various items, mostly interpolation using Printf # used when verbose option is selected +using RecipesBase # for plotting + +# - Visualization Packages - # +import Colors.RGB +using LaTeXStrings +using ProgressMeter #---------------------------------# # EXPORTS # @@ -56,8 +65,7 @@ export c4b # - Types - # # Inputs -export DuctedRotor, - Rotor, OperatingPoint, PanelingConstants, ReferenceParameters +export DuctedRotor, Rotor, OperatingPoint, PanelingConstants, ReferenceParameters # - Preallocations - # export allocate_prepost_container_cache, @@ -76,6 +84,7 @@ export ChainSolverOptions, SpeedMappingOptions, FixedPointOptions, CSORSolverOptions +export BoundaryLayerOptions # - Preprocess - # export setup_analysis @@ -83,6 +92,20 @@ export setup_analysis # - Analyses - # export analyze +# - Visualization - # +export generate_plots +export plotGeometry, + plotDuctGeometry, + plotBodyGeometry, + underlayGeometry, + plotCP, + plotVtan, + plotStagnation, + plotMomentum, + plotStreamlines, + staticPlots, + animatedPlots + #---------------------------------# # INCLUDES # #---------------------------------# @@ -99,6 +122,8 @@ include("utilities/caching/allocate_caches.jl") include("utilities/caching/reshape_caches.jl") include("utilities/caching/integration_caches.jl") include("utilities/caching/elliptic_grid_parameter_packaging.jl") +include("utilities/thermodynamics.jl") +include("utilities/ode_solvers.jl") # Airfoil utility functions include("utilities/airfoils/airfoil_utilities.jl") @@ -156,6 +181,16 @@ include("postprocess/velocities.jl") include("postprocess/pressures.jl") include("postprocess/rotor_performance.jl") include("postprocess/utils.jl") +include("postprocess/boundary_layer_utils.jl") +include("postprocess/boundary_layer_green.jl") +include("postprocess/boundary_layer_head.jl") +include("postprocess/viscous_drag.jl") + +##### ----- VISUALIZATION ----- ##### +# include("visualization/plot_recipe_defaults.jl") +include("visualization/plot_recipes.jl") +include("visualization/calculate_streamlines.jl") +include("visualization/convenience_plots.jl") ##### ----- DEBUGGING ----- ##### include("../test/test_utils.jl") diff --git a/src/analysis/analyses.jl b/src/analysis/analyses.jl index a15c9e67..5c5322b1 100644 --- a/src/analysis/analyses.jl +++ b/src/analysis/analyses.jl @@ -185,6 +185,7 @@ function analyze( solve_parameter_cache_dims, operating_point, reference_parameters, + options.boundary_layer_options, A_bb_LU, airfoils, idmaps, @@ -457,7 +458,9 @@ function analyze_multipoint( if options.verbose println("\n Operating Point:") for fn in fieldnames(typeof(operating_point)) - println(@sprintf " %6s = %5.3e" fn getfield(operating_point, fn)[]) + if fn!=:units + println(@sprintf " %6s = %5.3e" fn getfield(operating_point, fn)[]) + end end end @@ -468,7 +471,9 @@ function analyze_multipoint( # - copy over operating point - # for f in fieldnames(typeof(operating_point)) + if f != :units solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + end end # - Set up Body Linear System RHS - # diff --git a/src/analysis/setup.jl b/src/analysis/setup.jl index fcfeda1c..56dc4e1f 100644 --- a/src/analysis/setup.jl +++ b/src/analysis/setup.jl @@ -118,8 +118,10 @@ function setup_analysis( # copy over operating point for f in fieldnames(typeof(operating_point)) + if f!=:units solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) end +end ##### ----- PERFORM PREPROCESSING COMPUTATIONS ----- ##### if options.verbose diff --git a/src/postprocess/boundary_layer.jl b/src/postprocess/boundary_layer.jl deleted file mode 100644 index 41352b7d..00000000 --- a/src/postprocess/boundary_layer.jl +++ /dev/null @@ -1,664 +0,0 @@ -using FLOWMath -using Plots -using DifferentialEquations - -#---------------------------------# -# Pre-solve Functions # -#---------------------------------# - -""" -""" -function sa1(alt; hardness=5) - # Sigmoid blend of fits for <11000 and <25000 - T1(alt) = 15.04 - 0.00649 * alt + 273.1 #units: K - P1(alt) = 101.29 * (T1(alt) / 288.08)^5.256 #units: kPa - T2(alt) = -56.46 + 273.1 #units: K - P2(alt) = 22.65 * exp(1.73 - 0.000157 * alt) #units: kPa - - return FLOWMath.sigmoid_blend(T1(alt), T2(alt), alt, 11000, hardness), - FLOWMath.sigmoid_blend(P1(alt), P2(alt), alt, 11000, hardness) -end - -""" -""" -function sa2(alt; hardness=5) - # Sigmoid blend of <25000 and >25000 - T2(alt) = -56.46 + 273.1 #units: K - P2(alt) = 22.65 * exp(1.73 - 0.000157 * alt) #units: kPa - T3(alt) = -131.21 + 0.00299 * alt + 273.1 #units: K - P3(alt) = 2.488 / ((T3(alt) / 216.6)^11.388) #units: kPa - - return FLOWMath.sigmoid_blend(T2(alt), T3(alt), alt, 25000, hardness), - FLOWMath.sigmoid_blend(P2(alt), P3(alt), alt, 25000, hardness) -end - -""" -""" -function ideal_gas_rho(P, T) - return @. P / (0.2869 * T) -end - -""" - standard_atmosphere(alt; hardness=25) - -Smoothed fits to the Standard Atmosphere model. - -Assumes calorically imperfect gas. - -# Arguments -- `alt::Float` : Altitude - -# Keyword Arguments: -- `hardness::float` : hardness factor for sigmoid blend - -# Returns -- `static_temperature::Float` : Static temperature -- `static_pressure::Float` : Static pressure -- `static_density::Float` : Static density -""" -function standard_atmosphere(alt; hardness=25) - #Get temperature (T) and pressure (P) from table fits - if alt < (11000 + 25000) / 2.0 - T, P = sa1(alt; hardness=hardness) - else - T, P = sa2(alt; hardness=hardness) - end - - # return T, P, rho - return T, P * 1000, ideal_gas_rho(P, T) -end - -""" -""" -function speed_of_sound(static_pressure, static_density; gamma=1.4) - return sqrt(gamma * static_pressure / static_density) -end - -""" -""" -function calc_edge_mach(edge_velocity, speed_of_sound) - return edge_velocity / speed_of_sound -end - -""" -""" -function total_pressure(static_pressure, M; gamma=1.4) - return static_pressure * (1.0 + (gamma - 1.0) * M^2 / 2.0)^(gamma / (gamma - 1.0)) -end - -""" -""" -function total_temperature(static_temperature, M; gamma=1.4) - return static_temperature * (1.0 + (gamma - 1.0) * M^2 / 2.0) -end - -# """ -# """ -# function t0_ts(minf; gamma=1.4) -# return 1.0 + minf^2 * (gamma - 1.0) / 2.0 -# end - -# """ -# assumes calorically perfect gas -# """ -# function static_temperature(M, minf; gamma=1.4) -# return t0_ts(minf; gamma=gamma) / t0_ts(M; gamma=gamma) -# end - -""" -""" -function static_pressure(total_pressure, M; gamma=1.4) - return total_pressure / ((1.0 + (gamma - 1.0) * M^2 / 2.0)^(gamma / (gamma - 1.0))) -end - -""" -""" -function static_temperature(total_temperature, M; gamma=1.4) - return total_temperature / (1.0 + (gamma - 1.0) * M^2 / 2.0) -end - -""" -""" -function static_density(static_pressure, speed_of_sound; gamma=1.4) - return gamma * static_pressure / speed_of_sound^2 -end - -""" -""" -function calc_edge_viscosity( - static_temperure, mu_sea_level=0.0000181206, T_sea_level=288.15, S=110.4 -) - return mu_sea_level * (static_temperure / T_sea_level)^(3.0 / 2.0) * (T_sea_level + S) / - (static_temperure + S) -end - -#---------------------------------# -# solve functions # -#---------------------------------# -#= -note: function variables are in the following order: -- state variables -- variables dependent on state variables -- precomputed parameters -=# - -""" -""" -function Fc(M) - return sqrt(1.0 + 0.2 * M^2) -end - -""" -""" -function FR(M) - return 1.0 + 0.056 * M^2 -end - -""" -""" -function calc_Re(rhoe, Ue, d2, mue) - return rhoe * Ue * d2 / mue -end - -""" -""" -function calc_Cf0(Red2, M; hardness=5) - logblend = FLOWMath.sigmoid_blend( - 1.05, log(10, FR(M) * Red2), FR(M) * Red2, 10^1.02, hardness - ) - - return (0.01013 / (logblend - 1.02) - 0.00075) / Fc(M) -end - -""" -TODO: is there ever any way this could end up with a zero in the denominator? -""" -function calc_H12bar0(Cf0, M) - return 1.0 / (1.0 - 6.55 * sqrt(Cf0 / 2 * (1.0 + 0.04 * M^2))) -end - -""" -TODO: need to make sure that H12bar0 is never initialzied to zero -""" -function calc_Cf(H12bar, H12bar0, Cf0) - return Cf0 * (0.9 / (H12bar / H12bar0 - 0.4) - 0.5) -end - -""" -""" -function calc_H12(H12bar, M, Pr=1.0) - return (H12bar + 1.0) * (1.0 + Pr^(1.0 / 3.0) * M^2 / 5) - 1.0 -end - -""" -""" -function calc_Ctau(CE, Cf0, M) - return (0.024 * CE + 1.2 * CE^2 + 0.32Cf0) * (1.0 + 0.1 * M^2) -end - -""" -TODO: CE cannot be -0.01 or we'll have a divide by zero. -Green talks about handling this. -""" -function calc_F(CE, Cf0) - return (0.02 * CE + CE^2 + 0.8 * Cf0 / 3) / (0.01 + CE) -end - -""" -TODO: H12bar can't be 1, or we divide by zero -can H12bar ever get that low? -""" -function calc_H1(H12bar) - return 3.15 + 1.72 / (H12bar - 1.0) - 0.01 * (H12bar - 1.0)^2 -end - -""" -TODO: at what value(s) of H12bar does the denominator go to zero? -""" -function calc_dH12bardH1(H12bar) - return -(H12bar - 1.0)^2 / (1.72 + 0.02 * (H12bar - 1.0)^3) -end - -""" -TODO: can't evaluate this when r=0 or we'll divide by zero -""" -function calc_richardson_number(H12bar, d2, H12, H1, R) - return 2.0 * d2 / (3.0 * R) * (H12 + H1) * (H1 / H12bar + 0.3) -end - -""" - 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 beta(Ri, hardness=100.0) - return flowmath.sigmoid_blend(4.5, 7.0, Ri, 0.0, hardness) -end - -""" -""" -function longitudinal_curvature_influence(M, Ri) - return 1 + beta(Ri) * (1.0 + M^2 / 5) * Ri -end - -""" -""" -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 - -""" -TODO: Ue can't be zero or we divide by zero -""" -function dilation_influence(H12bar, d2, H12, H1, M, Ue, dUedx) - return 1.0 + 7.0 / 3.0 * M^2 * (H1 / H12bar + 1.0) * (H12 + H1) * d2 * dUedx / Ue -end - -""" -""" -function calc_lambda(args...) - return *(args...) -end - -""" -TODO; make sure H12 can't be zero -""" -function calc_d2dUedxUeeq0(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 - -""" -""" -function calc_CEeq0(H1, H12, Cf, d2dUedxUeeq0) - return H1 * (Cf / 2.0 - (H12 + 1.0) * d2dUedxUeeq0) -end - -""" -""" -function calc_Ctaueq0(CEeq0, Cf0, M) - return (0.24 * CEeq0 + 1.2 * CEeq0^2 + 0.32 * Cf0) * (1.0 + 0.1 * M^2) -end - -""" -TODO: c can be negative here. -""" -function calc_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 - -""" -TODO: make sure H12 can't be -1 -""" -function calc_d2dUedxUeeq(H1, H12, Cf, CEeq) - return (Cf / 2.0 - CEeq / H1) / (H12 + 1.0) -end - - -#---------------------------------# -# Viscous Drag # -#---------------------------------# - -""" -""" -function squire_young(d2,Ue, Uinf, H12, c) - return 2.0*d2/c *(Ue/Uinf)^((5.0+H12)/2.0) -end - -#---------------------------------# -# state initilization functions # -#---------------------------------# - -""" -initialize momentum thickness state with flat plate schlichting solution given a dx from the stagnation point to the next panel edge. -""" -function d2_init(dx, Rex) - return 0.036 * dx / Rex^0.2 -end - -""" -todo: initialize with H12bar0 -""" -function H12bar_init(Cf0, M) - return calc_H12bar0(Cf0, M) -end - -""" -todo: initialize with CEeq -""" -function CE_init(Ctaueq0, M, Cf0) - return calc_CEeq(Ctaueq0, M, 1, Cf0) -end - -""" -""" -function f(u, params, s) - du = zeros(length(u)) - return f!(du, u, params, s) -end - -""" -""" -function f!(du, u, params, s) - - # - unpack variables - # - rd2, H12bar, CE = u - (; - r_coords, - edge_velocity, - edge_density, - edge_viscosity, - edge_mach, - edge_acceleration, - surface_curvature, - surface_derivative, - ) = params - - # - rename for convenience - # - r = r_coords(s) - Ue = edge_velocity(s) - rhoe = edge_density(s) - mue = edge_viscosity(s) - M = edge_mach(s) - dUedx = edge_acceleration(s) - R = surface_curvature(s) - drds = surface_derivative(s) - - # - calculate intermediate variables - # - d2 = rd2 / r - Red2 = calc_Re(d2, Ue, rhoe, mue) - Cf0 = calc_Cf0(Red2, M) - H12bar0 = calc_H12bar0(Cf0, M) - Cf = calc_Cf(H12bar, H12bar0, Cf0) - H12 = calc_H12(H12bar, M) - - dH12bardH1 = calc_dH12bardH1(H12bar) #yes, this should be negative according to example - H1 = calc_H1(H12bar) - - F = calc_F(CE, Cf0) - - d2dUedxUeeq0 = calc_d2dUedxUeeq0(H12bar, H12, Cf, M) - CEeq0 = calc_CEeq0(H1, H12, Cf, d2dUedxUeeq0) - Ctaueq0 = calc_Ctaueq0(CEeq0, Cf0, M) - - Ctau = calc_Ctau(CE, Cf0, M) - - if params.lambda - if params.longitudinal_curvature - Ri = calc_richardson_number(H12bar, d2, H12, H1, r) - l1 = longitudinal_curvature_influence(M, Ri) - else - l1 = 1 - end - - if params.lateral_strain - l2 = lateral_strain_influence(H12bar, d2, H12, H1, r, drds) - else - l2 = 1 - end - - if params.dilation - l3 = dilation_influence(H12bar, d2, H12, H1, M, Ue, dUedx) - else - l3 = 1 - end - - lambda = calc_lambda(l1, l2, l3) - else - lambda = 1 - end - CEeq = calc_CEeq(Ctaueq0, M, lambda, Cf0) - d2dUedxUeeq = calc_d2dUedxUeeq(H1, H12, Cf, CEeq) - - # - system of equations - # - - # momentum - du[1] = r * Cf / 2.0 - (H12 + 2.0 - M^2) * rd2 * dUedx / Ue - - # entrainment - du[2] = dH12bardH1 * (CE - H1 * (Cf / 2.0 - (H12 + 1.0) * d2 * dUedx / Ue)) / d2 - - # lag - du[3] = - F / d2 * ( - 2.8 / (H12 + H1) * (sqrt(Ctaueq0) - lambda * sqrt(Ctau)) + d2dUedxUeeq - - d2 * dUedx / Ue * (1.0 + 0.075 * M^2 * ((1.0 + 0.2 * M^2) / (1.0 + 0.1 * M^2))) - ) - - return du, Cf, H12 -end - -""" -""" -function initialize_turbulent_boundary_layer_states(r_init, rhoe, Ue, dx, mue, M) - - # - Initialize States - # - # initialize momentum thickness (d2) using flat plate schlichting model - Rex = calc_Re(rhoe, Ue, dx, mue) - d2 = d2_init(dx, Rex) - - # initialize H12bar (compressible shape factor) using _0 equations - Red2 = calc_Re(rhoe, Ue, d2, mue) - Cf0 = calc_Cf0(Red2, M) - H12bar0 = H12bar_init(Cf0, M) - - # initialize the entrainment coefficient (CE) using equilibrium equations - H1 = calc_H1(H12bar0) - H12 = calc_H12(H12bar0, M) - Cf = calc_Cf(H12bar0, H12bar0, Cf0) - d2dUedxUeeq0 = calc_d2dUedxUeeq0(H12bar0, H12, Cf, M) - CEeq0 = calc_CEeq0(H1, H12, Cf, d2dUedxUeeq0) - Ctaueq0 = calc_Ctaueq0(CEeq0, Cf0, M) - CE = CE_init(Ctaueq0, M, Cf0) - - # initial states - return [r_init * d2, H12bar0, CE], Cf, H12 -end - -""" -""" -function RK2(f, u, params, s, ds) - k1, _, _ = f(u, params, s) - k2, Cf, H12 = f(u .+ (ds / 2) .* k1, params, s + (ds / 2)) - unext =@. u + k2 * ds - return unext, Cf, H12 -end - -""" -""" -function RK4(f, u, params, s, ds) - k1, Cf1, H12_1 = f(u, params, s) - k2, Cf2, H12_2 = f(u .+ (ds / 2) .* k1, params, s + (ds / 2)) - k3, Cf3, H12_3 = f(u .+ (ds / 2) .* k2, params, s + (ds / 2)) - k4, Cf4, H12_4 = f(u .+ ds .* k3, params, s + ds) - unext = @. u + (k1 + k2 * 2 + k3 * 2 + k4) * ds / 6 - Cfnext = (Cf1 + Cf2 * 2 + Cf3 * 2 + Cf4) / 6 - H12next = (H12_1 + H12_2 * 2 + H12_3 * 2 + H12_4) / 6 - return unext, Cfnext, H12next -end - -""" -""" -function solve_turbulent_boundary_layer_rk4!(f, init, params, svec; s_init_default=0.01) - u0, Cf0, H12_0 = init - s = svec[1] > s_init_default ? [svec[1]] : [s_init_default] - sep = [false] - sepid = [1] - us = zeros(eltype(u0), length(u0), length(svec)) - Cfs = zeros(eltype(u0), length(svec)) - H12s = zeros(eltype(u0), length(svec)) - us[:, 1] = u0 - Cfs[1] = Cf0 - H12s[1] = H12_0 - - for i in 1:(length(svec) - 1) - - # take step - us[:, i + 1], Cfs[i + 1], H12s[i+1] = RK4(f, us[:, i], params, svec[i], svec[i + 1] - svec[i]) - - if Cfs[i + 1] <= 0.0 - sep[1] = true - sepid[1] = i - break - end - end - - if sep[1] == true - u1sep = FLOWMath.akima(Cfs[(sepid - 1):sepid], us[1, (sepid - 1):sepid], 0.0) - u2sep = FLOWMath.akima(Cfs[(sepid - 1):sepid], us[2, (sepid - 1):sepid], 0.0) - u3sep = FLOWMath.akima(Cfs[(sepid - 1):sepid], us[3, (sepid - 1):sepid], 0.0) - usep = [u1sep; u2sep; u3sep] - H12sep = FLOWMath.akima(Cfs[(sepid - 1):sepid], H12s[(sepid - 1):sepid], 0.0) - ssep = FLOWMath.akima(Cfs[(sepid - 1):sepid], svec[(sepid - 1):sepid], 0.0) - else - usep = us[:, end] - H12sep = H12s[end] - ssep = svec[end] - end - - # return states at separate, and separation state, and separation postition - return us, usep, Cfs, H12sep, ssep, sepid[1] -end - -""" -""" -function solve_turbulent_boundary_layer(f!, u0, xspan, params) - prob = DifferentialEquations.ODEProblem(f!, u0, xspan, params) - # alg = DifferentialEquations.radau() - alg = DifferentialEquations.Tsit5() - sol = DifferentialEquations.solve(prob, alg) - return sol -end - -""" -""" -function test_flat_plate(s0, L; ds=1e-2) - - # - Parameters - # - lambda = false - longitudinal_curvature = false - lateral_strain = false - dilation = false - s = range(s0 + ds, s0 + L; step=ds) - N = length(s) - r_coords = ones(N) - drds = zeros(N) - Uinf = 10.0 - edge_velocity = Uinf * ones(N) - edge_acceleration = zeros(N) - surface_curvature = zeros(N) - surface_derivative = zeros(N) - alt = 0.0 - - Tinf, Pinf, rhoinf = standard_atmosphere(alt; hardness=25) - - a = speed_of_sound(Pinf, rhoinf) - - edge_mach = calc_edge_mach(edge_velocity, a) - Minf = calc_edge_mach(Uinf, a) - - P0 = total_pressure(Pinf, Minf) - - T0 = total_temperature(Tinf, Minf) - - Pe = static_pressure.(Ref(P0), edge_mach) - - Te = static_temperature.(Ref(T0), edge_mach) - - edge_density = static_density(Pe, a) - - edge_viscosity = calc_edge_viscosity.(Te) - - # - Spline things in terms of curve length - # - params = (; - lambda, - longitudinal_curvature, - lateral_strain, - dilation, - edge_velocity=FLOWMath.Akima(s, edge_velocity), - edge_mach=FLOWMath.Akima(s, edge_mach), - edge_density=FLOWMath.Akima(s, edge_density), - edge_viscosity=FLOWMath.Akima(s, edge_viscosity), - edge_acceleration=FLOWMath.Akima(s, edge_acceleration), - surface_curvature=FLOWMath.Akima(s, surface_curvature), - surface_derivative=FLOWMath.Akima(s, surface_derivative), - r_coords=FLOWMath.Akima(s, r_coords), - Uinf, - ) - - u0, Cf0, H12_0 = initialize_turbulent_boundary_layer_states( - r_coords[2], - edge_density[1], - edge_velocity[1], - s[1], - edge_viscosity[1], - edge_mach[1], - ) - - # return solve_turbulent_boundary_layer(f!, u0, (s[1], s[end]), params) - return u0, Cf0, H12_0, params, s -end - -ds = 1e-4 -s0 = 0.01 -L = 1.0 -u0, Cf0, H12_0, params, svec = test_flat_plate(s0, L; ds=ds) - -""" -""" -function schlichting(x, Ue, rhoe, mue) - Rex = rhoe * Ue * x / mue - d1 = 0.046 * x / Rex^0.2 - d2 = 0.036 * x / Rex^0.2 - cf = 0.0592 / Rex^0.2 - return d1, d2, cf -end - -schs = svec -sch = - schlichting.( - schs, - params.edge_velocity.(schs), - params.edge_density.(schs), - params.edge_viscosity.(schs), - ) - -sd1 = getindex.(sch, 1) -sd2 = getindex.(sch, 2) -scf = getindex.(sch, 3) - -plotstep = Int(1 / ds * L * 5 * 1e-3) -p1 = plot(; xlabel="s", ylabel="momentum thickness") -plot!( - p1, - schs[1:plotstep:end], - sd2[1:plotstep:end]; - label="Schlicting", - linewidth=3, - linestyle=:dot, -) - -us, usep, Cfs, H12sep, ssep, separation_id = solve_turbulent_boundary_layer_rk4!( - f, [u0, Cf0, H12_0], params, svec; s_init_default=0.01 -) - -cd = squire_young(usep[1],params.edge_velocity(ssep), params.Uinf, H12sep, L) - -plot!(p1, svec[1:plotstep:end], us[1, 1:plotstep:end]; label="Green RK4") -savefig("mom_thick_schlichting_green.png") - -# pcf = plot(xlabel="s", ylabel="Cf") -# plot!(pcf, svec[1:plotstep:end], Cfs[1:plotstep:end]; label="Green RK4") diff --git a/src/postprocess/boundary_layer_green.jl b/src/postprocess/boundary_layer_green.jl new file mode 100644 index 00000000..e5535776 --- /dev/null +++ b/src/postprocess/boundary_layer_green.jl @@ -0,0 +1,660 @@ +""" + 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) + # 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) + 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 FLOWMath.ksmax([21.0; 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) + # 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) +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) + # apply smooth-max to make sure we stay non-negative. + return FLOWMath.ksmax([0.0; 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 = FLOWMath.ksmax([0.0; 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, FLOWMath.ksmax([1.0; H12sep]), s_sep, sepid +end diff --git a/src/postprocess/boundary_layer_head.jl b/src/postprocess/boundary_layer_head.jl new file mode 100644 index 00000000..19810d2d --- /dev/null +++ b/src/postprocess/boundary_layer_head.jl @@ -0,0 +1,244 @@ +""" + setup_boundary_layer_functions_head( + 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_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 +""" +function setup_boundary_layer_functions_head( + 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) + + # local density + edge_mach(ss) = calculate_mach(edge_velocity(ss), operating_point.asound[]) + 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 + + # # r's + # r_coords = Akima_smooth(s, duct_control_points[2, :]) + + return (; edge_velocity, edge_acceleration, edge_density, edge_viscosity, verbose) +end + +""" + calculate_H(H1) + +Calculate the value of the shape factor used in Head's method. +""" +function calculate_H(H1) + + # get each side of the piecewise equation + hgeq = 0.86 * (H1 - 3.3)^(-0.777) + 1.1 + hlt = 1.1538 * (H1 - 3.3)^(-0.326) + 0.6778 + + # blend the pieces smoothly + return FLOWMath.sigmoid_blend(hlt, hgeq, H1, 5.3) +end + +""" + calculate_cf(H, Red2) + +Calculate the skin friction coefficient used in Head's method +""" +function calculate_cf(H, Red2) + return 0.246 * 10^(-0.678 * H) * Red2^(-0.268) +end + +""" + boundary_layer_residual_head(y, s, parameters) + +Out of place residual function for Head's method. +""" +function boundary_layer_residual_head(y, s, parameters) + dy = similar(y) .= 0 + return boundary_layer_residual_head!(dy, y, s, parameters) +end + +""" + boundary_layer_residual_head!(dy, y, s, parameters) + +Calculate dy give the current states, y, the input position, s, and various parameters. +""" +function boundary_layer_residual_head!(dy, y, s, parameters) + + # - unpack parameters - # + (; verbose) = parameters + + # - unpack variables - # + d2, H1 = y + verbose && printdebug("d2: ", d2) + verbose && printdebug("H1: ", H1) + + # limit H1 to be greater than 3.3 + H1lim = FLOWMath.ksmax([H1; 3.3 + 1e-2]) + verbose && printdebug("H1lim: ", H1lim) + + # - unpack variables - # + (; edge_velocity, edge_acceleration, edge_density, edge_viscosity) = parameters + + # - Intermediate Calculations - # + # determine dUedx + dUedx = edge_acceleration(s) + verbose && printdebug("dUedx: ", dUedx) + + # determine H + H = calculate_H(H1lim) + verbose && printdebug("H: ", H) + + # determine local edge velocity + Ue = edge_velocity(s) + verbose && printdebug("Ue: ", Ue) + + # determine momentum thickness reynolds number + Red2 = calculate_Re(edge_density(s), Ue, d2, edge_viscosity(s)) + verbose && printdebug("Red2: ", Red2) + + # determine cf + cf = calculate_cf(H, Red2) + verbose && printdebug("cf: ", cf) + + # - "Residuals" - # + # dδ₂/dx + dy[1] = cf / 2.0 - dUedx * d2 / Ue * (H + 2.0) + verbose && printdebug("dy[1]: ", dy[1]) + + # dH₁/dx + 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 +end + +""" + solve_head_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_head_boundary_layer!(f, rk, initial_states, steps, parameters; verbose=false) + + # Unpack States and variables for viscous drag + u0, H0 = 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)) + Hs = zeros(eltype(u0), length(steps)) + + # Initialize intermediate states and outputs + us[:, 1] = u0 + Hs[1] = H0 + + # 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], Hs[i + 1] = rk( + f, us[:, i], steps[i], abs(steps[i + 1] - steps[i]), parameters + ) + + sepid[1] = i + 1 + if Hs[i + 1] >= parameters.separation_criteria + sep[1] = true + break + end + end + + if sep[1] == true + + # - Interpolate to find actual s_sep - # + + usep = [ + FLOWMath.linear( + Hs[(sepid[] - 1):sepid[]], + us[1, (sepid[] - 1):sepid[]], + parameters.separation_criteria, + ) + FLOWMath.linear( + Hs[(sepid[] - 1):sepid[]], + us[2, (sepid[] - 1):sepid[]], + parameters.separation_criteria, + ) + ] + + Hsep = FLOWMath.linear( + Hs[(sepid[] - 1):sepid[]], + Hs[(sepid[] - 1):sepid[]], + parameters.separation_criteria, + ) + + s_sep = FLOWMath.linear( + Hs[(sepid[] - 1):sepid[]], + steps[(sepid[] - 1):sepid[]], + parameters.separation_criteria, + ) + stepsol = steps[1:sepid[]] + usol = us[:, 1:sepid[]] + else + usep = us[:, end] + Hsep = Hs[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, Hsep, s_sep, sepid, usol, stepsol +end diff --git a/src/postprocess/boundary_layer_utils.jl b/src/postprocess/boundary_layer_utils.jl new file mode 100644 index 00000000..0e80318f --- /dev/null +++ b/src/postprocess/boundary_layer_utils.jl @@ -0,0 +1,126 @@ +#---------------------------------# +# Setup functions # +#---------------------------------# + +""" + arc_lengths_from_panel_lengths(duct_panel_lengths, bl_ids) + +Cumulative sum of panel lengths for the given section of surface associated with the upper or lower boundary layer. + +# Arguments: +- `duct_panel_lengths::Vector{Float}` : vector of panel lengths (called influence_length in body_vortex_panels) associated with the duct (casing + nacelle). + +# Returns: +- `s::Vector{Float}` : cumulative sum of panel lengths between control points in the given index range, starting from zero. +""" +function arc_lengths_from_panel_lengths(duct_panel_lengths) + return cumsum( + [ + 0.0 + duct_panel_lengths[1] # partial panel length + [ + 0.5 * (duct_panel_lengths[i] + duct_panel_lengths[i - 1]) for + i in 3:length(duct_panel_lengths) + ] + ], + ) +end + +""" + split_at_stagnation_point(duct_panel_lengths, n_panels_casing) + +Split the duct body surface at the leading edge of the duct. + +# Arguments: +- `duct_panel_lengths::Vector{Float}` : Vector of panel lengths for the duct from casing trailing edge clockwise to nacelle trailing edge. +- `n_panels_casing::Int` : number of panels comprising the casing side of the duct + +# Returns: +- `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 + for i in 2: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 + end + dp[1] = dp[2] + end + + # interpolate the lengths between control points + sum_length = 0.5 * sum(duct_panel_lengths[stag_ids]) + stag_interp = FLOWMath.linear(dp, [0.0, sum_length], 0.0) + + partial_panel_lengths = [stag_interp, sum_length - stag_interp] + + s_upper = arc_lengths_from_panel_lengths( + [abs(partial_panel_lengths[2]); duct_panel_lengths[stag_ids[2]:end]] + ) + + s_lower = arc_lengths_from_panel_lengths( + [abs(partial_panel_lengths[1]); duct_panel_lengths[stag_ids[1]:-1:1]] + ) + + return s_upper, s_lower, stag_ids, stag_interp / sum_length +end + +""" + bl_step_fun(n, m, p) + +Function used in determining step sizes for boundary layer calculation. f(n) = m*n^p + +Given a number of steps, n ∈ [1:N], provides the cumulative step lengths according to the power, p, and the multiplicative factor, m; where p determined from the `set_boundary_layer_steps` functions. +""" +function bl_step_fun(n, m, p) + return m .* n .^ p +end + +""" + set_boundary_layer_steps(N::Int, first_step_size, total_length) + +Sets boundary layer steps based on desired number of steps (must be an Integer), an initial step size, and the total cumulative length of the steps. + +# Arguments: +- `N::Int` : Number of steps to take +- `first_step_size::Float` : size of first step (which is `m` in `bl_step_fun`) +- `total_length::Float` : total surface length to divide up. + +# Returns: +- `steps::Vector{Float}` : steps along surface length satisfying the equation: f(n) = m*n^p with the condition that `m` is the first step size and f(N) = `total_length` +""" +function set_boundary_layer_steps(N::Int, first_step_size, total_length) + + # residual + function res(y, x, p) + total_length = x[1] + first_step_size = x[2] + N = x[3] + return total_length - first_step_size * N^y + end + + # solver (uses Roots.jl, for which the docs indicate implicit differentiation is required) + function solvewrap(x, p) + reswrap(y) = res(y, x, p) + return Roots.find_zero(reswrap, 1.0) + end + + # solve for power coefficient + p = ImplicitAD.implicit(solvewrap, res, [total_length; first_step_size; N]) + + # return steps + return bl_step_fun(1:N, first_step_size, p) +end + +function set_boundary_layer_steps(N::Int, total_length) + return range(0.0, total_length; length=N) +end diff --git a/src/postprocess/postprocess.jl b/src/postprocess/postprocess.jl index 9e28be0b..a6267946 100644 --- a/src/postprocess/postprocess.jl +++ b/src/postprocess/postprocess.jl @@ -8,6 +8,7 @@ solve_parameter_cache_dims, operating_point, reference_parameters, + boundary_layer_options, A_bb_LU, airfoils, idmaps, @@ -31,6 +32,7 @@ Post-process a converged nonlinear solve solution. - `solve_parameter_cache_dims::NamedTuple` : the dimensions of the solver parameters - `operating_point::OperatingPoint` : the operating point being analyzed - `reference_parameters::ReferenceParameters` : a ReferenceParameters object +- `BoundaryLayerOptions::BoundaryLayerOptions` : a BoundaryLayerOptions object - `A_bb_LU::LinearAlgebra.LU` : LinearAlgebra LU factorization of the LHS matrix - `airfoils::Vector{AFType}` : A matrix of airfoil types associated with each of the blade elements - `idmaps::NamedTuple` : A named tuple containing index mapping used in bookkeeping throughout solve and post-process @@ -77,6 +79,17 @@ Post-process a converged nonlinear solve solution. - `vtan_nacelle_out` - `vtan_centerbody_in` - `vtan_centerbody_out` + - `boundary_layers` + - `stagnation_indices` + - `upper_solved_states` + - `upper_solved_steps` + - `lower_solved_states` + - `lower_solved_steps` + - `surface_length_upper` + - `surface_length_lower` + - `split_ratio` + - `separation_point_ratio_upper` + - `separation_point_ratio_lower` - `rotors` - `circulation` - `panel_strengths` @@ -143,6 +156,7 @@ function post_process( solve_parameter_cache_dims, operating_point, reference_parameters, + boundary_layer_options, A_bb_LU, airfoils, idmaps, @@ -164,12 +178,12 @@ function post_process( (; Vref, Rref) = reference_parameters # - Extract PrePost Cache - # - reset_containers!(prepost_containers; exception_keys=(:panels,:ivb)) + reset_containers!(prepost_containers; exception_keys=(:panels, :ivb)) (; - # stuff from pre-process + # Stuff from Pre-process panels, ivb, - # rotor stuff + # Rotor Stuff rotor_inviscid_thrust, rotor_inviscid_thrust_dist, rotor_viscous_thrust, @@ -193,14 +207,13 @@ function post_process( blade_normal_force_per_unit_span, blade_tangential_force_per_unit_span, blade_loading_intermediate_containers, - # body stuff + # Body Stuff zpts, vtan_tuple, cp_tuple, body_thrust, body_force_coefficient, - # cp_tuple, - # totals stuff + # Totals Stuff total_thrust, total_torque, total_power, @@ -209,6 +222,9 @@ function post_process( total_CT, total_CQ, total_CP, + # Boundary Layer Stuff #TODO: add these to caches at some point (requires re-work of boundary layer implementation likey) + # duct_viscous_drag, + # boundary_layer_outputs, ) = prepost_containers # - Extract Panels - # @@ -255,7 +271,7 @@ function post_process( blade_elements, wakeK, idmaps, - multipoint_index + multipoint_index, ) (; @@ -487,9 +503,32 @@ function post_process( Vref=Vref[1], ) + # - Duct Viscous Drag - # + + if boundary_layer_options.model_drag + + #TODO; make this in place + # compute_viscous_drag_duct!(duct_viscous_drag, boundary_layer_outputs, + duct_viscous_drag, boundary_layer_outputs = compute_viscous_drag_duct( + boundary_layer_options, + Vtan_out[1:Int(body_vortex_panels.npanel[1])], + Vtot_out[:, 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=verbose, + ) + + body_thrust[1] -= duct_viscous_drag + else + boundary_layer_outputs = nothing + end + ### --- TOTAL OUTPUTS --- ### # - Total Thrust - # + total_thrust[] = sum([rotor_inviscid_thrust'; rotor_viscous_thrust']) total_thrust[] = sum([rotor_inviscid_thrust'; rotor_viscous_thrust'; body_thrust]) # - Total Torque - # @@ -563,6 +602,8 @@ function post_process( vtan_nacelle_out, vtan_centerbody_in, vtan_centerbody_out, + # boundary layers + boundary_layers = boundary_layer_outputs, ), # - Rotor Values - # rotors=(; @@ -626,10 +667,7 @@ function post_process( res_vals.vr_wake, res_vals.Cm_avg, ), - reference_values=(; - Vinf = operating_point.Vinf[], - Vref = reference_parameters.Vref[], - ), + reference_values=(; Vinf=operating_point.Vinf[], Vref=reference_parameters.Vref[]), ) if write_outputs @@ -695,47 +733,11 @@ Run through the residual function post-convergence to save needed intermediate v - `res_vals::NamedTuple` : A named tuple containing the state variables and populated solve containers. """ function run_residual!( -solver_options::TS, -converged_states, -state_dims, -solve_container_cache, -solve_container_cache_dims, -operating_point, -ivr, -ivw, -linsys, -blade_elements, -wakeK, -idmaps, -multipoint_index -) where {TS<:ExternalSolverOptions} - -#= - NOTE: we want to get all the intermediate values available to user if desired. - The solve_containers cache will contain all the intermediate values after running the estimate states function. -=# -# - Separate out the state variables - # -vz_rotor, vtheta_rotor, Cm_wake = extract_state_variables( - solver_options, converged_states, state_dims -) - -# - Extract and Reset Cache - # -# get cache vector of correct types -solve_container_cache_vec = @views PreallocationTools.get_tmp( - solve_container_cache, converged_states -) -solve_containers = withdraw_solve_container_cache( - solver_options, solve_container_cache_vec, solve_container_cache_dims -) -reset_containers!(solve_containers) #note: also zeros out state estimates - -# - Estimate New States - # -# currently has 280 allocations -estimate_states!( - solve_containers, - vz_rotor, - vtheta_rotor, - Cm_wake, + solver_options::TS, + converged_states, + state_dims, + solve_container_cache, + solve_container_cache_dims, operating_point, ivr, ivw, @@ -743,54 +745,55 @@ estimate_states!( blade_elements, wakeK, idmaps, -) + multipoint_index, +) where {TS<:ExternalSolverOptions} + + #= + NOTE: we want to get all the intermediate values available to user if desired. + The solve_containers cache will contain all the intermediate values after running the estimate states function. + =# + # - Separate out the state variables - # + vz_rotor, vtheta_rotor, Cm_wake = extract_state_variables( + solver_options, converged_states, state_dims + ) + + # - Extract and Reset Cache - # + # get cache vector of correct types + solve_container_cache_vec = @views PreallocationTools.get_tmp( + solve_container_cache, converged_states + ) + solve_containers = withdraw_solve_container_cache( + solver_options, solve_container_cache_vec, solve_container_cache_dims + ) + reset_containers!(solve_containers) #note: also zeros out state estimates + + # - Estimate New States - # + # currently has 280 allocations + estimate_states!( + solve_containers, + vz_rotor, + vtheta_rotor, + Cm_wake, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps, + ) -return (; vz_rotor, vtheta_rotor, Cm_wake, solve_containers...) + return (; vz_rotor, vtheta_rotor, Cm_wake, solve_containers...) end """ """ function run_residual!( -solver_options::CSORSolverOptions, -converged_states, -state_dims, -solve_container_cache, -solve_container_cache_dims, -operating_point, -ivr, -ivw, -linsys, -blade_elements, -wakeK, -idmaps, -multipoint_index -) - -#= - NOTE: we want to get all the intermediate values available to user if desired. - The solve_containers cache will contain all the intermediate values after running the insides of the residual function -=# -# - Separate out the state variables - # -Gamr, sigr, gamw = extract_state_variables(solver_options, converged_states, state_dims) - -# - Extract and Reset Cache - # -# get cache vector of correct types -solve_container_cache_vec = @views PreallocationTools.get_tmp( - solve_container_cache, converged_states -) -solve_container_cache_vec .= 0 -solve_containers = withdraw_solve_container_cache( - solver_options, solve_container_cache_vec, solve_container_cache_dims -) - -# - Run Residual - # -compute_CSOR_residual!( - zeros(eltype(converged_states), 2), - solver_options, - solve_containers, - Gamr, - sigr, - gamw, + solver_options::CSORSolverOptions, + converged_states, + state_dims, + solve_container_cache, + solve_container_cache_dims, operating_point, ivr, ivw, @@ -798,8 +801,43 @@ compute_CSOR_residual!( blade_elements, wakeK, idmaps, - multipoint_index; - verbose=false, + multipoint_index, +) + + #= + NOTE: we want to get all the intermediate values available to user if desired. + The solve_containers cache will contain all the intermediate values after running the insides of the residual function + =# + # - Separate out the state variables - # + Gamr, sigr, gamw = extract_state_variables(solver_options, converged_states, state_dims) + + # - Extract and Reset Cache - # + # get cache vector of correct types + solve_container_cache_vec = @views PreallocationTools.get_tmp( + solve_container_cache, converged_states + ) + solve_container_cache_vec .= 0 + solve_containers = withdraw_solve_container_cache( + solver_options, solve_container_cache_vec, solve_container_cache_dims + ) + + # - Run Residual - # + compute_CSOR_residual!( + zeros(eltype(converged_states), 2), + solver_options, + solve_containers, + Gamr, + sigr, + gamw, + operating_point, + ivr, + ivw, + linsys, + blade_elements, + wakeK, + idmaps, + multipoint_index; + verbose=false, ) return (; Gamr, sigr, gamw, solve_containers...) diff --git a/src/postprocess/viscous_drag.jl b/src/postprocess/viscous_drag.jl new file mode 100644 index 00000000..b790e522 --- /dev/null +++ b/src/postprocess/viscous_drag.jl @@ -0,0 +1,299 @@ +#---------------------------------# +# Viscous Drag # +#---------------------------------# + +""" + squire_young(d2, Ue, Uinf, H12) + +Squire-Young formula for the viscous drag coeffiecient of one side of a 2D body. + +# Arguments: +- `d2::Float` : Momentum thickness at separation extrapolated back to the trailing edge (d2 = d2sep+rsep-rTE) +- `Ue::Float` : Edge velocity at separation point +- `Uinf::Float` : Freestream velocity +- `H12::Float` : Boundary layer shape factor at separation point + +Note: if no separation occurs, the inputs are simply the final values for the boundary layer. + +# Returns: +- `cdc::Float` : viscous drag coefficient times reference chord +""" +function squire_young(d2, Ue, Uinf, H12) + # note: formula also divides by chord, but we're going to multiply by chord later. + return 2.0 * d2 * (Ue / Uinf)^((5.0 + H12) / 2.0) +end + +""" + total_viscous_drag_duct(cd_upper, cd_lower, exit_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 +- `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) + # 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 +end + +#---------------------------------# +# HEAD'S METHOD # +#---------------------------------# + +""" + compute_single_side_drag_coefficient_head( + steps, + exit_radius, + operating_point, + boundary_layer_options; + verbose=false, + ) + +Solve integral boundary layer and obtain viscous drag coefficient from Squire-Young formula for one side of the duct (side being defined as portion of surface on once side of the stagnation point) + +# 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 +- `boundary_layer_options::BoundaryLayerOptions` : BoundaryLayerOptions object + +# Returns: +- `cd::Float` : viscous drag coefficient +""" +function compute_single_side_drag_coefficient_head( + steps, + exit_radius, + operating_point, + boundary_layer_functions, + boundary_layer_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 + + usep, Hsep, s_sep, sepid, usol, stepsol = solve_head_boundary_layer!( + boundary_layer_residual_head, + boundary_layer_options.rk, + [initial_states, H0], + steps, + (; boundary_layer_functions..., boundary_layer_options.separation_criteria); + verbose=false, + ) + + # verbose && println(s_sep / steps[end]) + + cdsq = squire_young( + usep[1], boundary_layer_functions.edge_velocity(s_sep), operating_point.Vinf[], Hsep + ) + + cd = FLOWMath.ksmin([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], + s_sep, + ) + ], + 100, + ) + + cd += cdadd + + return cd, usol, stepsol, s_sep / steps[end] +end + +""" + compute_viscous_drag_duct( + boundary_layer_options::BoundaryLayerOptions, + Vtan_duct, + cp_duct, + duct_panel_lengths, + exit_radius, + operating_point, + ) + +Determine total, dimensional viscous drag on the duct. + +# Arguments: +- `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 + +# Returns: +- `duct_viscous_drag::Float` : total viscous drag of duct +- `boundary_layer_outputs::NamedTuple` : named tuple of various boundary layer related outputs: + - `stagnation_indices` + - `upper_solved_states` + - `upper_solved_steps` + - `lower_solved_states` + - `lower_solved_steps` + - `surface_length_upper` + - `surface_length_lower` + - `split_ratio` + - `separation_point_ratio_upper` + - `separation_point_ratio_lower` +""" +function compute_viscous_drag_duct( + boundary_layer_options::HeadsBoundaryLayerOptions, + Vtan_duct, + Vtot_duct, + duct_panel_lengths, + duct_panel_tangents, + exit_radius, + operating_point; + 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 + ) + + # set up boundary layer solve parameters + boundary_layer_functions_upper = setup_boundary_layer_functions_head( + s_upper, + [0.0; Vtan_duct[stag_ids[2]:end]], + operating_point, + boundary_layer_options; + verbose=verbose, + ) + + # set up boundary layer solve parameters + boundary_layer_functions_lower = setup_boundary_layer_functions_head( + s_lower, + [0.0; Vtan_duct[stag_ids[1]:-1:1]], + operating_point, + boundary_layer_options; + verbose=verbose, + ) + + # - Set integration steps - # + # upper side + upper_steps = + set_boundary_layer_steps( + boundary_layer_options.n_steps, + boundary_layer_options.first_step_size, + s_upper[end] - boundary_layer_options.offset, + ) .+ boundary_layer_options.offset + + # 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 + + # - Get drag coeffients - # + + cdc_upper, usol_upper, stepsol_upper, s_sep_upper = compute_single_side_drag_coefficient_head( + upper_steps, + exit_radius, + operating_point, + boundary_layer_functions_upper, + (; + boundary_layer_options.rk, + boundary_layer_options.separation_criteria, + separation_allowance=boundary_layer_options.separation_allowance_upper, + separation_penalty=boundary_layer_options.separation_penalty_upper, + ); + verbose=verbose, + ) + + cdc_lower, usol_lower, stepsol_lower, s_sep_lower = compute_single_side_drag_coefficient_head( + lower_steps, + exit_radius, + operating_point, + 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, + ); + verbose=verbose, + ) + + # Calculate total viscous drag + total_drag = total_viscous_drag_duct( + cdc_upper, cdc_lower, exit_radius, operating_point.Vinf[], operating_point.rhoinf[] + ) + + # Return + return total_drag, + (; + stagnation_indices=stag_ids, + upper_solved_states=usol_upper, + upper_solved_steps=stepsol_upper, + lower_solved_states=usol_lower, + lower_solved_steps=stepsol_lower, + surface_length_upper=s_upper, + surface_length_lower=s_lower, + split_ratio, + separation_point_ratio_upper=s_sep_upper, + separation_point_ratio_lower=s_sep_lower, + ) +end + +#---------------------------------# +# SCHLICHTING APPROXIMATION # +#---------------------------------# + +""" + compute_viscous_drag_duct_schlichting( + vtan_duct_TE, duct_chord, TE_radius, operating_point + ) + +Computes Schlichting approximation of skin friction drag dimensionalized to drag per unit length using duct chord, and using the trailing edge circuference as the total length. + +# Arguments: +- `vtan_duct_TE::Float` : tangential velocity at the duct trailing edge +- `duct_chord::Float` : length of duct +- `TE_radius::Float` : radius of the trailing edge point +- `operating_point::OperatingPoint` : OperatingPoint object +""" +function compute_viscous_drag_duct_schlichting( + vtan_duct_TE, duct_chord, TE_radius, operating_point +) + (; rhoinf, muinf) = operating_point + q = 0.5 * rhoinf[] * vtan_duct_TE .^ 2 + + Cf = + 0.074 ./ + calculate_Re.(Ref(rhoinf[]), vtan_duct_TE, Ref(duct_chord), Ref(muinf[])) .^ 0.2 + + return (Cf[1] * q[1] + Cf[2] * q[2]) * duct_chord * 2.0 * pi * TE_radius, Cf +end diff --git a/src/preprocess/geometry/body_geometry.jl b/src/preprocess/geometry/body_geometry.jl index 59bcb284..80c9a1ab 100644 --- a/src/preprocess/geometry/body_geometry.jl +++ b/src/preprocess/geometry/body_geometry.jl @@ -96,7 +96,7 @@ function reinterpolate_bodies!( # check that the splining didn't put any of the center body radial coordinates in the negative. for rpcb in eachcol(rp_centerbody_coordinates) - if rpcb[2] < 0.0 && rpcb[2] > -2.0*eps() + if rpcb[2] < 0.0 && rpcb[2] > -2.0 * eps() rpcb[2] = 0.0 end end diff --git a/src/preprocess/geometry/panel_geometry.jl b/src/preprocess/geometry/panel_geometry.jl index 996cddc3..08c99722 100644 --- a/src/preprocess/geometry/panel_geometry.jl +++ b/src/preprocess/geometry/panel_geometry.jl @@ -41,15 +41,15 @@ function generate_panels( # - Initialize Outputs - # # control points - controlpoint = zeros(TF, 2, totpanel[]) # x-r, panel + controlpoint = zeros(TF, 2, totpanel[]) # z-r, panel # nodes - node = zeros(TF, 2, totnode[]) # x-r, node + node = zeros(TF, 2, totnode[]) # z-r, node # node map - nodemap = zeros(Int, 2, totpanel[]) # node, x-r + nodemap = zeros(Int, 2, totpanel[]) # node, z-r # first and last nodes of bodies, rotor, or wakes #TODO: check where this is used and update order of indices if needed - endnodes = zeros(TF, nbodies[], 2, 2) # TE, upper-lower, x-r + endnodes = zeros(TF, nbodies[], 2, 2) # TE, upper-lower, z-r # indices of endnodess endnodeidxs = ones(Int, 2, nbodies[]) # lower idx, upper idx, lower or upper # indices of endpanels @@ -68,8 +68,8 @@ function generate_panels( ittangent = zeros(TF, 2, nbodies[]) # trailing edge gap - tenode = zeros(TF, nbodies[], 2, 2) # body, node1-2, x-r - tenormal = zeros(TF, 2, nbodies[]) #body, x-r + tenode = zeros(TF, nbodies[], 2, 2) # body, node1-2, z-r + tenormal = zeros(TF, 2, nbodies[]) #body, z-r teinfluence_length = zeros(TF, nbodies[]) teadjnodeidxs = similar(endnodeidxs) .= 1 #can't be the same as endpoints, because we may have repeated values for non-duct bodies tendotn = zeros(TF, 2, nbodies[]) #bodies, node1,2 @@ -330,7 +330,7 @@ function def_te_panel!( endnodeidxs, ) for ib in 1:Int(nbodies[]) - # check if signs of x-tangents are the same + # check if signs of z-tangents are the same if sign(tangent[1, Int(endpanelidxs[1, ib])]) != sign(tangent[1, Int(endpanelidxs[2, ib])]) # if not: it's a duct diff --git a/src/utilities/caching/allocate_caches.jl b/src/utilities/caching/allocate_caches.jl index 9439dc79..cc296702 100644 --- a/src/utilities/caching/allocate_caches.jl +++ b/src/utilities/caching/allocate_caches.jl @@ -722,12 +722,18 @@ function allocate_solve_parameter_cache( l = lfs(s) Vinf = cache_dims!(total_length, l, s) + Minf = cache_dims!(total_length, l, s) + rhoinf = cache_dims!(total_length, l, s) muinf = cache_dims!(total_length, l, s) asound = cache_dims!(total_length, l, s) + Ptot = cache_dims!(total_length, l, s) + + Ttot = cache_dims!(total_length, l, s) + s = (nrotor,) l = lfs(s) Omega = cache_dims!(total_length, l, s) @@ -827,7 +833,7 @@ function allocate_solve_parameter_cache( Gamr, sigr, gamw, - operating_point=(; Vinf, rhoinf, muinf, asound, Omega), + operating_point=(; Vinf, Minf, rhoinf, muinf, asound, Ptot, Ttot, Omega), ivr=(; v_rb, v_rr, v_rw), ivw=(; v_wb, v_wr, v_ww), linsys=(; A_bb, b_bf, A_bw, A_pw, A_br, A_pr), @@ -899,12 +905,18 @@ function allocate_solve_parameter_cache( l = lfs(s) Vinf = cache_dims!(total_length, l, s) + Minf = cache_dims!(total_length, l, s) + rhoinf = cache_dims!(total_length, l, s) muinf = cache_dims!(total_length, l, s) asound = cache_dims!(total_length, l, s) + Ptot = cache_dims!(total_length, l, s) + + Ttot = cache_dims!(total_length, l, s) + s = (nrotor,) l = lfs(s) Omega = cache_dims!(total_length, l, s) @@ -1004,7 +1016,7 @@ function allocate_solve_parameter_cache( vz_rotor, vtheta_rotor, Cm_wake, - operating_point=(; Vinf, rhoinf, muinf, asound, Omega), + operating_point=(; Vinf, Minf, rhoinf, muinf, asound, Ptot, Ttot, Omega), ivr=(; v_rb, v_rr, v_rw), ivw=(; v_wb, v_wr, v_ww), linsys=(; A_bb, b_bf, A_bw, A_pw, A_br, A_pr), diff --git a/src/utilities/caching/reshape_caches.jl b/src/utilities/caching/reshape_caches.jl index d43fc51b..e5e7f267 100644 --- a/src/utilities/caching/reshape_caches.jl +++ b/src/utilities/caching/reshape_caches.jl @@ -623,6 +623,9 @@ function withdraw_solve_parameter_cache(solver_options::TS, vec, dims) where {TS Vinf=reshape( @view(vec[dims.operating_point.Vinf.index]), dims.operating_point.Vinf.shape ), + Minf=reshape( + @view(vec[dims.operating_point.Minf.index]), dims.operating_point.Minf.shape + ), rhoinf=reshape( @view(vec[dims.operating_point.rhoinf.index]), dims.operating_point.rhoinf.shape ), @@ -632,6 +635,12 @@ function withdraw_solve_parameter_cache(solver_options::TS, vec, dims) where {TS asound=reshape( @view(vec[dims.operating_point.asound.index]), dims.operating_point.asound.shape ), + Ptot=reshape( + @view(vec[dims.operating_point.Ptot.index]), dims.operating_point.Ptot.shape + ), + Ttot=reshape( + @view(vec[dims.operating_point.Ttot.index]), dims.operating_point.Ttot.shape + ), Omega=reshape( @view(vec[dims.operating_point.Omega.index]), dims.operating_point.Omega.shape ), @@ -759,6 +768,9 @@ function withdraw_solve_parameter_cache( Vinf=reshape( @view(vec[dims.operating_point.Vinf.index]), dims.operating_point.Vinf.shape ), + Minf=reshape( + @view(vec[dims.operating_point.Minf.index]), dims.operating_point.Minf.shape + ), rhoinf=reshape( @view(vec[dims.operating_point.rhoinf.index]), dims.operating_point.rhoinf.shape ), @@ -768,6 +780,12 @@ function withdraw_solve_parameter_cache( asound=reshape( @view(vec[dims.operating_point.asound.index]), dims.operating_point.asound.shape ), + Ptot=reshape( + @view(vec[dims.operating_point.Ptot.index]), dims.operating_point.Ptot.shape + ), + Ttot=reshape( + @view(vec[dims.operating_point.Ttot.index]), dims.operating_point.Ttot.shape + ), Omega=reshape( @view(vec[dims.operating_point.Omega.index]), dims.operating_point.Omega.shape ), diff --git a/src/utilities/inputs.jl b/src/utilities/inputs.jl index 8bc83b98..ec82614b 100644 --- a/src/utilities/inputs.jl +++ b/src/utilities/inputs.jl @@ -1,36 +1,165 @@ """ - OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) + struct Imperial + +Disptach type used for setting Operating Point Units to Imperial + +# Units: feet, pounds, seconds +- Velocity: feet per second +- Temperature: Fahrenheit +- Pressure: pounds per square foot +- Density: slugs per cubic foot +- Dynamic Viscosity: slugs per square foot +""" +struct Imperial end + +""" + struct SI + +Disptach type used for setting Operating Point Units to SI. +Note that the default is SI, and is really only used in the backend. + +# Units: meters, kilograms, seconds +- Velocity: meters per second +- Temperature: Kelvin +- Pressure: Pascals +- Density: kilograms per meter cubed +- Dynamic Viscosity: Pascal seconds +""" +struct SI end + +""" + OperatingPoint(Vinf, Minf, rhoinf, muinf, asound, Ptot, Ttot, Omega) + OperatingPoint( + Vinf, Omega, rhoinf=nothing, muinf=nothing, asound=nothing; altitude=0.0 + ) + OperatingPoint( + ::Imperial, Vinf, Omega, rhoinf=nothing, muinf=nothing, asound=nothing; altitude=0.0 + ) DuctedRotor operating point information. -# Arguments +Functions that take in `altitude` will populate undefined thermodynamic properties of the freestream using a standard_atmosphere model, ideal gas law, and Sutherland's law; defaulting to SI units. +If the `::Imperial` dispatch type is input, then the thermodynamic properties will be converted to Imperial units. + +# Fields/Arguments - `Vinf::AbstractVector{Float}` : Freestream velocity magnitude (which is only in the axial direction). +- `Minf::AbstractVector{Float}` : Freestream Mach number - `rhoinf::AbstractVector{Float}` : Freestream density - `muinf::AbstractVector{Float}` : Freestream viscosity - `asound::AbstractVector{Float}` : Freestream speed of sound +- `Ptot::AbstractVector{Float}` : Freestream total pressure +- `Ttot::AbstractVector{Float}` : Freestream total temperature - `Omega::AbstractVector{Float}` : Rotor rototation rate(s) """ struct OperatingPoint{ - Tv<:AbstractVector, - Tr<:AbstractVector, - Tm<:AbstractVector, Ta<:AbstractVector, + TM<:AbstractVector, + Tm<:AbstractVector, To<:AbstractVector, + Tp<:AbstractVector, + Tr<:AbstractVector, + Tt<:AbstractVector, + Tv<:AbstractVector, + Tu, } + units::Tu Vinf::Tv + Minf::TM rhoinf::Tr muinf::Tm asound::Ta + Ptot::Tp + Ttot::Tt Omega::To end -function OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) +function OperatingPoint( + units::Imperial, + Vinf, + Omega, + rhoinf=nothing, + muinf=nothing, + asound=nothing, + altitude=0.0, +) + + # Get thermodynamic properties + Tinf, Pinf, rho_inf, mu_inf = standard_atmosphere(units, altitude) + + # freestream density + if isnothing(rhoinf) + rhoinf = rho_inf + end + + # freestream dynamic viscosity + if isnothing(muinf) + muinf = mu_inf + end + + # freestream speed of sound + if isnothing(asound) + asound = speed_of_sound(Pinf, rhoinf) + end + + # freestream Mach + Minf = calculate_mach.(Vinf, asound) + + # freestream total pressure and temperature + Ptot = total_pressure.(Pinf, Minf) + Ttot = total_temperature.(Tinf, Minf) + + return OperatingPoint( + Imperial(), + isscalar(Vinf) ? [Vinf] : Vinf, + isscalar(Minf) ? [Minf] : Minf, + isscalar(rhoinf) ? [rhoinf] : rhoinf, + isscalar(muinf) ? [muinf] : muinf, + isscalar(asound) ? [asound] : asound, + isscalar(Ptot) ? [Ptot] : Ptot, + isscalar(Ttot) ? [Ttot] : Ttot, + isscalar(Omega) ? [Omega] : Omega, + ) +end + +function OperatingPoint( + Vinf, Omega, rhoinf=nothing, muinf=nothing, asound=nothing; altitude=0.0 +) + + # Get thermodynamic properties + Tinf, Pinf, rho_inf, mu_inf = standard_atmosphere(altitude) + + # freestream density + if isnothing(rhoinf) + rhoinf = rho_inf + end + + # freestream dynamic viscosity + if isnothing(muinf) + muinf = mu_inf + end + + # freestream speed of sound + if isnothing(asound) + asound = speed_of_sound(Pinf, rhoinf) + end + + # freestream Mach + Minf = calculate_mach.(Vinf, asound) + + # freestream total pressure and temperature + Ptot = total_pressure.(Pinf, Minf) + Ttot = total_temperature.(Tinf, Minf) + return OperatingPoint( + SI(), isscalar(Vinf) ? [Vinf] : Vinf, + isscalar(Minf) ? [Minf] : Minf, isscalar(rhoinf) ? [rhoinf] : rhoinf, isscalar(muinf) ? [muinf] : muinf, isscalar(asound) ? [asound] : asound, + isscalar(Ptot) ? [Ptot] : Ptot, + isscalar(Ttot) ? [Ttot] : Ttot, isscalar(Omega) ? [Omega] : Omega, ) end @@ -148,7 +277,6 @@ function Rotor( fliplift; i_know_what_im_doing=false, ) - if length(findall(t -> t > pi / 2, twists)) > 2 @warn "It looks like your input twist angles may be in degrees. Note that the required units for twist are radians. Converting to radians for you (set the `i_know_what_im_doing` keyword argument to true to disable automatic conversion)." twists .*= pi / 180.0 @@ -179,8 +307,8 @@ end - `duct_coordinates::AbstractMatrix` : The [z, r] coordinates of the duct geometry beginning at the inner (casing) side trailing edge and proceeding clockwise. Note that the duct geometry absolute radial position does not need to be included here if the `autoshiftduct` option is selected. - `centerbody_coordinates::AbstractMatrix` : The [z, r] coordinates of the centerbody beginning at the leading edge and ending at the trailing edge. Note that the leading edge is assumed to be placed at a radial distance of 0.0 from the axis of rotation. -- `paneling_constants::PanelingConstants` : Constants used in re-paneling the geometry. - `rotor::Rotor` : Rotor (and possibly stator) geometric paramters. +- `paneling_constants::PanelingConstants` : Constants used in re-paneling the geometry. """ struct DuctedRotor{Td<:AbstractMatrix,Tcb<:AbstractMatrix,Trp<:Rotor,Tpc<:PanelingConstants} duct_coordinates::Td diff --git a/src/utilities/misc.jl b/src/utilities/misc.jl index cf456053..00e6b0a3 100644 --- a/src/utilities/misc.jl +++ b/src/utilities/misc.jl @@ -31,6 +31,16 @@ function printval(text, val) return nothing end +""" + printdebug(variable_name, variable, nspaces=4) + +Formatted printing when you have lots of debugging to do and want things to line up nicely. +""" +function printdebug(variable_name, variable, nspaces=4) + @printf "%*s %-14s %2.12f\n" nspaces " " variable_name variable + # println("\t$(variable_name)" * v) + return nothing +end """ dot(A, B) = sum(a * b for (a, b) in zip(A, B)) @@ -189,3 +199,21 @@ function cache_dims!(total_length, l, s) total_length[] += l return dims end + +""" + akima_smooth(x, y, xpt; delta=2.0 * eps(), eps=eps()) + +Wrapper for FLOWMath.akima with different optional argument values. +""" +function akima_smooth(x, y, xpt; delta=2.0 * eps(), eps=eps()) + return FLOWMath.akima(x, y, xpt, delta, eps) +end + +""" + Akima_smooth(x, y; delta=2.0 * eps(), eps=eps()) + +Wrapper for FLOWMath.Akima with different optional argument values. +""" +function Akima_smooth(x, y; delta=2.0 * eps(), eps=eps()) + return FLOWMath.Akima(x, y, delta, eps) +end diff --git a/src/utilities/ode_solvers.jl b/src/utilities/ode_solvers.jl new file mode 100644 index 00000000..5b7c25d4 --- /dev/null +++ b/src/utilities/ode_solvers.jl @@ -0,0 +1,43 @@ +""" + RK2(f, y, s, ds, parameters) + +2nd Order Runge-Kutta integration scheme. + +# Arguments: +- `f::function_handle` : residual function for integration +- `y::Vector{Float}` : current states +- `s::Float` : current position +- `ds::Float` : step size +- `parameters::NamedTuple` : BoundaryLayerOptions and other various parameters +""" +function RK2(f, y, s, ds, parameters) + parameters.verbose && println(" 1st call:") + k1, _ = f(y, s, parameters) + parameters.verbose && println(" 2nd call:") + k2, aux = f(y .+ (ds / 2) .* k1, s + (ds / 2), parameters) + unext = @. y + k2 * ds + return unext, aux +end + +""" + RK4(f, y, s, ds, parameters) + +4th Order Runge-Kutta integration scheme. + +# Arguments: +- `f::function_handle` : residual function for integration +- `y::Vector{Float}` : current states +- `s::Float` : current position +- `ds::Float` : step size +- `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 +end + diff --git a/src/utilities/options.jl b/src/utilities/options.jl index ee97a60b..fa209747 100644 --- a/src/utilities/options.jl +++ b/src/utilities/options.jl @@ -62,6 +62,13 @@ Used in integration method dispatch """ abstract type IntegrationMethod end +""" + abstract type BoundaryLayerOptions + +Used in boundary layer method dispatch +""" +abstract type BoundaryLayerOptions end + #---------------------------------# # QUADRATURE TYPES # #---------------------------------# @@ -550,7 +557,8 @@ Options for Chain Solvers (try one solver, if it doesn't converge, try another) - `converged::AbstractArray{Bool} = [false]` : flag to track if convergence took place. - `iterations::AbstractArray{Int} = [0]` : iteration counter """ -@kwdef struct ChainSolverOptions{TB,TI,TS<:ExternalSolverOptions} <:ExternalPolyAlgorithmOptions +@kwdef struct ChainSolverOptions{TB,TI,TS<:ExternalSolverOptions} <: + ExternalPolyAlgorithmOptions solvers::AbstractArray{TS} = [ NLsolveOptions(; algorithm=:anderson, atol=1e-10, iteration_limit=200), MinpackOptions(; atol=1e-10, iteration_limit=100), @@ -644,6 +652,76 @@ Options for Newton elliptic grid solver. iterations::AbstractArray{TI} = [0] end +#---------------------------------# +# BOUNDARY LAYER TYPES # +#---------------------------------# +""" + struct HeadsBoundaryLayerOptions + +# Fields: +- `model_drag::Tb=false` : flag to turn on viscous drag approximation +- `n_steps::Int = Int(2e2)` : number of steps to use in boundary layer integration +- `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) +- `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 + model_drag::Tb=false + n_steps::Ti = Int(2e2) + first_step_size::Tf = 1e-6 + offset::To = 1e-3 + rk::Tfun = RK2 + separation_criteria::Ts=3.0 + separation_allowance_upper::Ti=10 + separation_allowance_lower::Ti=10 + separation_penalty_upper::Tp=0.2 + separation_penalty_lower::Tp=0.2 +end + +""" + struct GreensBoundaryLayerOptions + +NOTE: Green's method is mostly implemented, but there are several bugs still, especially when using Imperial units. +Known Bugs: +- Imperial units overestimate momentum thickness. Likely a unit conversion bug. +- In some cases of non-separation, the momentum thickens or shape parameter becomes exceedingly large, vastly overestimating the drag coefficient. Likely the product of one or more of the adjustments to try and make the method more robust. + +# Fields: +- `model_drag::Tb=true` : flag to turn off viscous drag approximation +- `lambda::Bool = true` : flag to add secondary influences into boundary layer residuals +- `longitudinal_curvature::Bool = true` : if `lambda`=true, flag to add longitudinal curvature influence into boundary layer residuals +- `lateral_strain::Bool = true` : if `lambda`=true, flag to add lateral strain influence into boundary layer residuals +- `dilation::Bool = true` : if `lambda`=true, flag to add dilation influence into boundary layer residuals +- `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) +- `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 + model_drag::Tb=true + lambda::Tb = false + longitudinal_curvature::Tb = true + lateral_strain::Tb = true + dilation::Tb = true + n_steps::Ti = Int(2e2) + first_step_size::Tf = 1e-6 + offset::To = 1e-3 + rk::Tfun = RK2 + separation_allowance_upper::Ti=25 + separation_allowance_lower::Ti=25 + separation_penalty_upper::Tp=0.2 + separation_penalty_lower::Tp=0.2 +end + #---------------------------------# # OPTION SET TYPES # #---------------------------------# @@ -670,6 +748,7 @@ Type containing (nearly) all the available user options. ## Integration Options - `integration_options::IntegrationOptions type = IntegrationOptions()` : integration options ## Post-processing Options +- `boundary_layer_options::BoundaryLayerOptions` : BoundaryLayerOptions object - `write_outputs::AbstractArray{Bool} = [false]` : Bool for whether to write the outputs of the analysis to an external file (slow) - `outfile::AbstractArray{String} = ["outputs.jl"]` : External output file name (including path information) for files to write - `checkoutfileexists::Bool = false` : Flag for whether to check if file exists before overwriting @@ -680,6 +759,7 @@ Type containing (nearly) all the available user options. """ @kwdef struct Options{ TB, + TBL, TBwo, TF, TI, @@ -705,6 +785,7 @@ Type containing (nearly) all the available user options. # - Integration Options - # integration_options::TIo = IntegrationOptions() # - Post-processing Options - # + boundary_layer_options::TBL = HeadsBoundaryLayerOptions() write_outputs::TBwo = [false] outfile::TSf = ["outputs.jl"] checkoutfileexists::TB = false diff --git a/src/utilities/thermodynamics.jl b/src/utilities/thermodynamics.jl new file mode 100644 index 00000000..af0e1a64 --- /dev/null +++ b/src/utilities/thermodynamics.jl @@ -0,0 +1,216 @@ +""" + sa1(altitude; hardness=50) + +Standard atmosphere temperature and pressure in SI units blended between the first linear portion and the constant portion. +""" +function sa1(altitude; hardness=50) + # Sigmoid blend of fits for <11000 and <25000 + T1(altitude) = 15.04 - 0.00649 * altitude + 273.1 #units: K + P1(altitude) = 101.29 * (T1(altitude) / 288.08)^5.256 #units: kPa + T2(altitude) = -56.46 + 273.1 #units: K + P2(altitude) = 22.65 * exp(1.73 - 0.000157 * altitude) #units: kPa + + return FLOWMath.sigmoid_blend(T1(altitude), T2(altitude), altitude, 11000, hardness), + FLOWMath.sigmoid_blend(P1(altitude), P2(altitude), altitude, 11000, hardness) +end + +""" + sa2(altitude; hardness=50) + +Standard atmosphere temperature and pressure in SI units blended between the the constant portion and the second linear portion. +""" +function sa2(altitude; hardness=50) + # Sigmoid blend of <25000 and >25000 + T2(altitude) = -56.46 + 273.1 #units: K + P2(altitude) = 22.65 * exp(1.73 - 0.000157 * altitude) #units: kPa + T3(altitude) = -131.21 + 0.00299 * altitude + 273.1 #units: K + P3(altitude) = 2.488 / ((T3(altitude) / 216.6)^11.388) #units: kPa + + return FLOWMath.sigmoid_blend(T2(altitude), T3(altitude), altitude, 25000, hardness), + FLOWMath.sigmoid_blend(P2(altitude), P3(altitude), altitude, 25000, hardness) +end + +""" + ideal_gas_rho(P, T) + +Ideal gas law for calculating density in SI units. +""" +function ideal_gas_rho(P, T) + return @. P / (0.2869 * T) +end + +""" + sutherlands_law( + static_temperure, mu_sea_level=1.789e-5, T_sea_level=288.15, S=110.4 + ) + +Sutherland's law in SI units for calculating air viscosity relative to sea level. +""" +function sutherlands_law( + static_temperure, mu_sea_level=1.789e-5, T_sea_level=288.15, S=110.4 +) + return mu_sea_level * (static_temperure / T_sea_level)^(3.0 / 2.0) * (T_sea_level + S) / + (static_temperure + S) +end + +""" + standard_atmosphere(altitude; hardness=25) + standard_atmosphere(::Imperial, altitude; hardness=25) + +Smoothed fits to the Standard Atmosphere model. + +Assumes calorically imperfect gas. + +# Arguments +- `altitude::Float` : Altitude in meters for SI units, or feet for Imperial units + +# Keyword Arguments: +- `hardness::float` : hardness factor for sigmoid blend + +# Returns +- `static_temperature::Float` : Static temperature +- `static_pressure::Float` : Static pressure +- `static_density::Float` : Static density +- `static_dynamic_viscosity::Float` : Static dynamic Viscosity +""" +function standard_atmosphere(altitude; hardness=25) + + #Get temperature (T) and pressure (P) from table fits + if altitude < (11000 + 25000) / 2.0 + T, P = sa1(altitude; hardness=hardness) + else + T, P = sa2(altitude; hardness=hardness) + end + + # return T, P, rho, mu + return T, P * 1000, ideal_gas_rho(P, T), sutherlands_law(T) +end + +function standard_atmosphere(imperial_units, altitude; hardness=25) + + # convert from feet to meters + altitude *= 0.3048 + + #Get temperature (T) and pressure (P) from table fits in SI units + if altitude < (11000 + 25000) / 2.0 + T, P = sa1(altitude; hardness=hardness) + else + T, P = sa2(altitude; hardness=hardness) + end + + # - Convert to Imperial Units - # + + # convert from kg/m^3 to slugs/ft^3 + rho = ideal_gas_rho(P, T) * 0.00194032 + + # convert from Pa-s to slugs/ft-s + mu = sutherlands_law(T) * 0.0208854342 + + # convert from Kelvin to Fahrenheit + T -= 273.15 + T *= 9.0 / 5.0 + T += 32.0 + + # convert from kilo pascals to slugs/ft^2 + P *= 20.885434273039 + + return T, P, rho, mu +end + +""" + speed_of_sound(static_pressure, static_density; gamma=1.4) + +Speed of sound from isentropic relations +""" +function speed_of_sound(static_pressure, static_density; gamma=1.4) + return sqrt(gamma * static_pressure / static_density) +end + +""" + calculate_mach(edge_velocity, speed_of_sound) + +Mach number from velocity and speed of sound +""" +function calculate_mach(edge_velocity, speed_of_sound) + return edge_velocity / speed_of_sound +end + +""" + total_pressure(static_pressure, M; gamma=1.4) + +Total pressure from isentropic relations +""" +function total_pressure(static_pressure, M; gamma=1.4) + return static_pressure * (1.0 + (gamma - 1.0) * M^2 / 2.0)^(gamma / (gamma - 1.0)) +end + +""" + total_temperature(static_temperature, M; gamma=1.4) + +Total temperature from isentropic relations +""" +function total_temperature(static_temperature, M; gamma=1.4) + return static_temperature * (1.0 + (gamma - 1.0) * M^2 / 2.0) +end + +""" + static_pressure(total_pressure, M; gamma=1.4) + +Static pressure from isentropic relations +""" +function static_pressure(total_pressure, M; gamma=1.4) + return total_pressure / ((1.0 + (gamma - 1.0) * M^2 / 2.0)^(gamma / (gamma - 1.0))) +end + +""" + static_temperature(total_temperature, M; gamma=1.4) + +Static temperature from isentropic relations +""" +function static_temperature(total_temperature, M; gamma=1.4) + return total_temperature / (1.0 + (gamma - 1.0) * M^2 / 2.0) +end + +""" + static_density(static_pressure, speed_of_sound; gamma=1.4) + +Static density from isentropic relations +""" +function static_density(static_pressure, speed_of_sound; gamma=1.4) + return gamma * static_pressure / speed_of_sound^2 +end + +""" + convert_temperature_to_kelvin(::Units, T) + +Convert from Fahrenheit to Kelvin or Return temperature if already in SI units +""" +function convert_temperature_to_kelvin(::SI, T) + return T +end + +function convert_temperature_to_kelvin(::Imperial, T) + + # convert to celsius + T -= 32.0 + T *= 5.0 / 9.0 + + # convert to kelvin + T += 273.15 + + return T +end + +""" + convert_viscosity(::SI, mu) + +Convert viscosity from Imperial units to SI or Return input if already SI units. +""" +function convert_viscosity(::SI, mu) + return mu +end + +function convert_viscosity(::Imperial, mu) + return mu * 0.0208854342 +end + diff --git a/src/visualization/calculate_streamlines.jl b/src/visualization/calculate_streamlines.jl new file mode 100644 index 00000000..a3b5ce92 --- /dev/null +++ b/src/visualization/calculate_streamlines.jl @@ -0,0 +1,187 @@ +function wake_streamlines() + + # set body streamlines + # determine centerbody wake streamline + # determine inner rotor streamlines + # determine duct wake streamline + + return nothing +end + +function get_streamline_velocities!( + velocities, + points, + z, + 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, + @view(points[:, z, :]), + body_vortex_panels.node, + body_vortex_panels.nodemap, + body_vortex_panels.influence_length, + integration_options, + ) + + velocities[z, :, 1] .+= v_bs[:, :, 1] * body_vortex_strengths + velocities[z, :, 2] .+= v_bs[:, :, 2] * body_vortex_strengths + + # wake-induced velocity + DuctAPE.induced_velocities_from_vortex_panels_on_points!( + v_ws, + @view(points[:, z, :]), + wake_vortex_panels.node, + wake_vortex_panels.nodemap, + wake_vortex_panels.influence_length, + integration_options, + ) + + velocities[z, :, 1] .+= v_ws[:, :, 1] * wake_vortex_strengths + velocities[z, :, 2] .+= v_ws[:, :, 2] * wake_vortex_strengths + + # rotor-induced velocity + DuctAPE.induced_velocities_from_source_panels_on_points!( + v_rs, + @view(points[:, z, :]), + rotor_source_panels.node, + rotor_source_panels.nodemap, + rotor_source_panels.influence_length, + integration_options, + ) + + velocities[z, :, 1] .+= v_rs[:, :, 1] * rotor_source_strengths + velocities[z, :, 2] .+= v_rs[:, :, 2] * rotor_source_strengths + + return nothing +end + +function calculate_streamlines( + body_vortex_panels, + body_vortex_strengths, + wake_vortex_panels, + wake_vortex_strengths, + rotor_source_panels, + rotor_source_strengths, + Vinf; + starting_radial_points=range(0.001, 1.0; length=20), + axial_range=[0.5, 1.0], + nominal_step_size=1e-2, + step_limit=Int(1e2), + integration_options=IntegrationOptions(), + dot_tol=0.999, + stag_tol=0.4, +) + + # TODO; + # figure out how to ignore things after a stagnation point is hit, maybe set it to NaN? + + TF = promote_type( + eltype(body_vortex_panels.node), + eltype(wake_vortex_panels.node), + eltype(rotor_source_panels.node), + ) + + # Initialize + points = zeros(TF, 2, step_limit, length(starting_radial_points)) + points[1, 1, :] .= axial_range[1] + points[2, 1, :] = starting_radial_points + velocities = zeros(TF, step_limit, length(starting_radial_points), 2) + + # So we don't have to allocate all the memory: + v_bs = zeros(TF, length(starting_radial_points), Int(body_vortex_panels.totnode[]), 2) + v_ws = zeros(TF, length(starting_radial_points), Int(wake_vortex_panels.totnode[]), 2) + v_rs = zeros(TF, length(starting_radial_points), Int(rotor_source_panels.totnode[]), 2) + + # March along each streamline + for z in 1:step_limit + velocities[z, :, 1] .= Vinf + + # - Calculate Velocity at Current Point - # + get_streamline_velocities!( + velocities, + points, + z, + 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, + ) + + if z > 1 + dp = [ + dot( + velocities[z, r, :] / norm(velocities[z, r, :]), + velocities[z - 1, r, :] / norm(velocities[z - 1, r, :]), + ) for r in eachindex(starting_radial_points) + ] + while all(dp .< dot_tol) + points[:, z, :] .= 0.5 .* (points[:, z, :] .+ points[:, z - 1, :]) + velocities[z, :, 1] .= Vinf + + get_streamline_velocities!( + velocities, + points, + z, + 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, + ) + + dp .= [ + dot( + velocities[z, r, :] / norm(velocities[z, r, :]), + velocities[z - 1, r, :] / norm(velocities[z - 1, r, :]), + ) for r in eachindex(starting_radial_points) + ] + end + end + + if z != step_limit + # - Take Step and Set Next Point - # + vmag = sqrt.(velocities[z, :, 2] .^ 2 .+ velocities[z, :, 1] .^ 2) + points[1, z + 1, :] .= + points[1, z, :] .+ velocities[z, :, 1] .* nominal_step_size ./ vmag + points[2, z + 1, :] .= + points[2, z, :] .+ velocities[z, :, 2] .* nominal_step_size ./ vmag + if z > 1 + for r in eachindex(starting_radial_points) + if vmag[r] < stag_tol + points[:, z + 1, r] .= NaN + end + end + end + end + end + + return points, velocities +end diff --git a/src/visualization/convenience_plots.jl b/src/visualization/convenience_plots.jl new file mode 100644 index 00000000..5e1e5eff --- /dev/null +++ b/src/visualization/convenience_plots.jl @@ -0,0 +1,389 @@ +#---------------------------------# +# DISPATCH TYPES # +#---------------------------------# +struct staticPlots end +struct animatedPlots end + +""" + generate_plots( + ::staticPlots, (or ::animatedPlots) + Plots, + ins, + outs; + save_path="", + static_file_type=".pdf", + suffix=nothing, + plot_geometry=true, + plot_pressure=false, + plot_velocity=false, + plot_boundary_layer=false, + plot_streamlines=false, + verbose=false, + kwargs..., + ) + +Generate standard suite of plots or animations from input and output objects. + +# Arguments: +- `::staticPlots (or ::animatedPlots)` : +- `Plots::` : the Plots package namespace +- `ins::NamedTuple` : returned inputs from `analyze` function +- `outs::Vector{NamedTuple}` : retured outputs from `analyze` function + +# Keyword Arguments: +- `save_path=""` : custom save path +- `static_file_type=".pdf"` : file type for static files (must be compatible with the desired backend) +- `suffix=nothing` : custom suffixes, if unused plots files will be numbered starting from 1. +- `plot_geometry=true` : flag to generate geometry plot +- `plot_panels=false` : flag to include markers indicating panel edges in geometry plot +- `plot_pressure=false` : flag to generate surface pressures plot +- `plot_velocity=false` : flag to generate surface velocities plot +- `plot_boundary_layer=false` : flag to generate boundary layer plot +- `plot_streamlines=false` : flag to generate streamlines plot +- `verbose=false` : print verbose statements +- `kwargs...` : arguments passed into the plot functions (Plots keyword arguments/defaults to be used in every plot) +""" +function generate_plots( + ::staticPlots, + Plots, + ins, + outs; + save_path="", + static_file_type=".pdf", + suffix=nothing, + plot_geometry=true, + plot_panels=false, + plot_pressure=false, + plot_velocity=false, + plot_boundary_layer=false, + plot_streamlines=false, + verbose=false, + kwargs..., +) + if isnothing(suffix) + suffix = ["_$(o)" for o in eachindex(outs)] + lpad.(suffix, ndigits(length(outs)) + 1, "0") + end + + # Extract useful items + (; body_vortex_panels, wake_vortex_panels, rotor_source_panels) = ins.panels + + # - Geometry - # + if plot_geometry + verbose && println("Plotting Geometry") + Plots.plot( + plotGeometry(), + body_vortex_panels, + rotor_source_panels, + wake_vortex_panels; + plot_panels=plot_panels, + kwargs..., + ) + Plots.savefig(save_path * "geometry" * static_file_type) + end + + # - Surface Pressure - # + if plot_pressure + verbose && println("Plotting Surface Pressure") + for (i, out) in zip(suffix, outs) + + # underlay geometry first + plt = Plots.plot( + underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs... + ) + + # then Plots.plot pressure distributions + Plots.plot!( + plotCP(), + body_vortex_panels, + out.bodies, + rotor_source_panels; + subplot=2, + inset=(1, Plots.bbox(0, 0, 1, 1)), + kwargs..., + ) + + Plots.savefig(save_path * "surface_pressure" * i * static_file_type) + end + end + + # - Tangential Velocity - # + if plot_velocity + verbose && println("Plotting Surface Velocity") + for (i, out) in zip(suffix, outs) + # underlay geometry first + plt = Plots.plot( + underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs... + ) + + # then Plots.plot velocity distributions + Plots.plot!( + plotVtan(), + body_vortex_panels, + out.bodies, + out.reference_values.Vref[], + rotor_source_panels; + subplot=2, + inset=(1, Plots.bbox(0, 0, 1, 1)), + kwargs..., + ) + Plots.savefig(save_path * "surface_velocity" * i * static_file_type) + end + end + + # - Boundary Layer Stuff - # + if plot_boundary_layer + verbose && println("Plotting Boundary Layer") + for (i, out) in zip(suffix, outs) + # Plots.plot momentum thicknesses on top + plt = Plots.plot( + plotDuctGeometry(), body_vortex_panels; color=6, linewidth=0.5, kwargs... + ) + + Plots.plot!( + plotMomentum(), + out.bodies.boundary_layers, + body_vortex_panels; + scale_thickness=5.0, + kwargs..., + ) + + # Plots.plot stagnation point on top + Plots.plot!( + plotStagnation(), + out.bodies.boundary_layers, + body_vortex_panels; + markersize=4, + color=3, + markerstrokecolor=3, + kwargs..., + ) + Plots.savefig(save_path * "boundary_layer" * i * static_file_type) + end + end + + # - Streamlines - # + if plot_streamlines + verbose && println("Plotting Streamlines") + for (i, out) in zip(suffix, outs) + plt = Plots.plot(plotBodyGeometry(), body_vortex_panels; kwargs...) + + Plots.plot!( + plotStreamlines(), + body_vortex_panels, + out.bodies.panel_strengths, + wake_vortex_panels, + out.wake.panel_strengths, + rotor_source_panels, + out.rotors.panel_strengths, + out.reference_values.Vinf[]; + starting_radial_points=range(0.01, 0.7; length=50), + axial_range=[-0.15, 0.5], + step_limit=75, + nominal_step_size=1e-2, + integration_options=IntegrationOptions(), + stag_tol=0.0, + kwargs..., + ) + + Plots.plot!( + plotStagnation(), + out.bodies.boundary_layers, + body_vortex_panels; + markersize=4, + color=2, + markerstrokecolor=2, + kwargs..., + ) + + Plots.savefig(save_path * "streamlines" * i * static_file_type) + end + end + + return nothing +end + +function generate_plots( + ::animatedPlots, + Plots, + ins, + outs; + save_path="", + static_file_type=".pdf", + plot_geometry=true, + plot_panels=false, + plot_pressure=false, + plot_velocity=false, + plot_boundary_layer=false, + plot_streamlines=false, + verbose=false, + kwargs..., +) + + # Extract useful items + (; body_vortex_panels, wake_vortex_panels, rotor_source_panels) = ins.panels + + # - Geometry - # + if plot_geometry + verbose && println("Plotting Geometry") + Plots.plot( + plotGeometry(), + body_vortex_panels, + rotor_source_panels, + wake_vortex_panels; + plot_panels=plot_panels, + kwargs..., + ) + Plots.savefig(save_path * "geometry" * static_file_type) + end + + # - Surface Pressure - # + if plot_pressure + verbose && println("Animating Surface Pressure") + anim = Plots.Animation() + if verbose + pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true) + end + for out in outs + + # underlay geometry first + plt = Plots.plot( + underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs... + ) + + # then Plots.plot pressure distributions + Plots.plot!( + plotCP(), + body_vortex_panels, + out.bodies, + rotor_source_panels, + wake_vortex_panels; + subplot=2, + inset=(1, Plots.bbox(0, 0, 1, 1)), + (; kwargs..., background_color=:white)..., + ) + + Plots.frame(anim, plt) + verbose && next!(pbar) + end + verbose && println("Saving $(save_path)surface_pressure.gif") + Plots.gif(anim, save_path * "surface_pressure.gif") + end + + # - Tangential Velocity - # + if plot_velocity + verbose && println("Animating Surface Velocity") + anim = Plots.Animation() + if verbose + pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true) + end + for out in outs + # underlay geometry first + plt = Plots.plot( + underlayGeometry(), body_vortex_panels, rotor_source_panels; kwargs... + ) + + # then Plots.plot velocity distributions + Plots.plot!( + plotVtan(), + body_vortex_panels, + out.bodies, + out.reference_values.Vref[], + rotor_source_panels, + wake_vortex_panels; + subplot=2, + inset=(1, Plots.bbox(0, 0, 1, 1)), + (; kwargs..., background_color=:white)..., + ) + Plots.frame(anim, plt) + verbose && next!(pbar) + end + verbose && println("Saving $(save_path)surface_velocity.gif") + Plots.gif(anim, save_path * "surface_velocity.gif") + end + + # - Boundary Layer Stuff - # + if plot_boundary_layer + verbose && println("Animating Boundary Layer") + anim = Plots.Animation() + if verbose + pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true) + end + for out in outs + # Plots.plot momentum thicknesses on top + plt = Plots.plot( + plotDuctGeometry(), body_vortex_panels; color=6, linewidth=0.5, kwargs... + ) + + Plots.plot!( + plotMomentum(), + out.bodies.boundary_layers, + body_vortex_panels; + scale_thickness=5.0, + kwargs..., + ) + + # Plots.plot stagnation point on top + Plots.plot!( + plotStagnation(), + out.bodies.boundary_layers, + body_vortex_panels; + color=3, + markerstrokecolor=3, + (; kwargs..., background_color=:white)..., + ) + Plots.frame(anim, plt) + verbose && next!(pbar) + end + verbose && println("Saving $(save_path)boundary_layer.gif") + Plots.gif(anim, save_path * "boundary_layer.gif") + end + + # - Streamlines - # + if plot_streamlines + verbose && println("Animating Streamlines") + if verbose + pbar = Progress(length(outs); desc="Generating Frames...", showspeed=true) + end + anim = Plots.Animation() + for out in outs + plt = Plots.plot(plotBodyGeometry(), body_vortex_panels; kwargs...) + + Plots.plot!( + plotStreamlines(), + body_vortex_panels, + out.bodies.panel_strengths, + wake_vortex_panels, + out.wake.panel_strengths, + rotor_source_panels, + out.rotors.panel_strengths, + out.reference_values.Vinf[]; + starting_radial_points=range(0.01, 0.7; length=50), + axial_range=[-0.15, 0.5], + step_limit=75, + nominal_step_size=1e-2, + integration_options=IntegrationOptions(), + stag_tol=0.0, + kwargs..., + ) + + Plots.plot!( + plotStagnation(), + out.bodies.boundary_layers, + body_vortex_panels; + markersize=2, + color=2, + markerstrokecolor=2, + (; kwargs..., background_color=:white)..., + ) + + Plots.frame(anim, plt) + + verbose && next!(pbar) + end + verbose && println("Saving $(save_path)streamlines.gif") + Plots.gif(anim, save_path * "streamlines.gif") + end + + return nothing +end diff --git a/src/visualization/plot_recipes.jl b/src/visualization/plot_recipes.jl new file mode 100644 index 00000000..9d3637f2 --- /dev/null +++ b/src/visualization/plot_recipes.jl @@ -0,0 +1,741 @@ +#---------------------------------# +# DISPATCH TYPES # +#---------------------------------# + +# Geometry +struct plotGeometry end +struct plotDuctGeometry end +struct plotBodyGeometry end +struct underlayGeometry end + +# Surface Distributions +struct plotCP end +struct plotVtan end + +# Boundary Layer +struct plotStagnation end +struct plotMomentum end + +# Streamlines +struct plotStreamlines end + +#---------------------------------# +# PLOT Duct GEOMETRY # +#---------------------------------# +@recipe function plot_duct_geometry( + ::plotDuctGeometry, + bvp; + plot_panels=false, + 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 + ), +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Default Labels + xguide --> L"z" + yguide --> L"r" + + # Aspect Ratio + aspect_ratio --> 1 + xlim --> ( + minimum(bvp.node[1, 1:Int(bvp.npanel[1])]) - maximum(bvp.node[1, :]) * 0.1, + maximum(bvp.node[1, :]) * 1.1, + ) + + # - Plot Body Geometry - # + @series begin + label --> false + seriescolor --> 1 + if plot_panels + marker --> true + markersize --> 0.75 + markerstrokecolor --> 1 + end + return bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])], + bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])] + end + + return nothing +end +# --------------------------------# +# PLOT BODY GEOMETRY # +#---------------------------------# +@recipe function plot_body_geometry( + ::plotBodyGeometry, + bvp; + plot_panels=false, + 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 + ), +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Default Labels + xguide --> L"z" + yguide --> L"r" + + # Aspect Ratio + aspect_ratio --> 1 + + # - Plot Body Geometry - # + for b in 1:Int(bvp.nbodies[]) + @series begin + label --> false + seriescolor --> 1 + if plot_panels + marker --> true + markersize --> 0.75 + markerstrokecolor --> 1 + end + return bvp.node[1, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])], + bvp.node[2, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])] + end + end + + return nothing +end + +#---------------------------------# +# PLOT ALL GEOMETRY # +#---------------------------------# +@recipe function plot_geometry( + ::plotGeometry, + bvp, + rsp, + wvp; + plot_panels=false, + discrete_labels=true, + xticktol=1e-2, + yticktol=1e-2, + tickdigits=2, + 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 + ), +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Default Labels + xguide --> L"z" + yguide --> L"r" + + # Aspect Ratio + aspect_ratio --> 1 + ylim --> (0.0, maximum(bvp.node[2, :]) * 1.05) + + if discrete_labels + xticks --> DuctAPE.determine_geometry_xlabels( + bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits + ) + yticks --> determine_geometry_ylabels(bvp, rsp; tol=yticktol, tickdigits=tickdigits) + xgrid --> true + ygrid --> true + else + grid --> false + end + + # - Plot Rotor Geometry - # + for r in 1:Int(rsp.nbodies[]) + @series begin + label --> false + seriescolor --> 2 + if plot_panels + markershape --> :square + markersize --> 2 + markercolor --> 2 + markerstrokecolor --> 2 + markerstrokealpha --> 0 + linewidth --> 0.5 + else + linewidth --> 2 + end + return rsp.node[1, Int(rsp.endnodeidxs[1, r]):Int(rsp.endnodeidxs[2, r])], + rsp.node[2, Int(rsp.endnodeidxs[1, r]):Int(rsp.endnodeidxs[2, r])] + end + end + + # - Plot Body Geometry - # + for b in 1:Int(bvp.nbodies[]) + @series begin + label --> false + seriescolor --> 1 + if plot_panels + markershape --> :circle + markersize --> 1.5 + markercolor --> 1 + markerstrokecolor --> 1 + markerstrokealpha --> 0 + linewidth --> 0.5 + else + linewidth --> 1.5 + end + + return bvp.node[1, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])], + bvp.node[2, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])] + end + end + + # - Plot Wake Geometry - # + for w in 1:Int(wvp.nbodies[]) + @series begin + label --> false + seriescolor --> 3 + linewidth --> 0.5 + if plot_panels + markershape --> :vline + markersize --> 1 + markercolor --> 3 + markerstrokecolor --> 3 + end + return wvp.node[1, Int(wvp.endnodeidxs[1, w]):Int(wvp.endnodeidxs[2, w])], + wvp.node[2, Int(wvp.endnodeidxs[1, w]):Int(wvp.endnodeidxs[2, w])] + end + end + + return nothing +end + +#---------------------------------# +# PLOT CP DISTRIBUTINS # +#---------------------------------# +@recipe function plot_cp( + ::plotCP, + bvp, + bout, + rsp=nothing, + wvp=nothing; + xticktol=1e-2, + tickdigits=2, + cp_ylim=nothing, + 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 + ), +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Flip y axis + yflip --> true + # Label Axes + yguide --> L"C_P" + xguide --> L"z" + + if !isnothing(rsp) && !isnothing(wvp) + xticks --> + determine_geometry_xlabels(bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits) + xgrid --> true + ygrid --> false + else + grid --> false + end + + if !isnothing(cp_ylim) + ylim --> cp_ylim + end + + # - Plot Body Surface Pressure - # + for b in 1:Int(bvp.nbodies[]) + @series begin + label --> false + seriescolor --> b + linewidth = 2 + + return bvp.controlpoint[ + 1, Int(bvp.endpanelidxs[1, b]):Int(bvp.endpanelidxs[2, b]) + ], + bout.cp_out[Int(bvp.endpanelidxs[1, b]):Int(bvp.endpanelidxs[2, b])] + end + end +end + +#---------------------------------# +# PLOT Vtan DISTRIBUTIONS # +#---------------------------------# +@recipe function plot_vtan( + ::plotVtan, + bvp, + bout, + Vref, + rsp=nothing, + wvp=nothing; + xticktol=1e-2, + tickdigits=2, + vtan_ylim=nothing, + 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 + ), +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Label Axes + yguide --> L"\frac{V_\mathrm{surf}}{V_\infty}" + xguide --> L"z" + if !isnothing(vtan_ylim) + ylim --> vtan_ylim + end + + if !isnothing(rsp) && !isnothing(wvp) + xticks --> + determine_geometry_xlabels(bvp, rsp, wvp; tol=xticktol, tickdigits=tickdigits) + xgrid --> true + ygrid --> false + else + grid --> false + end + + # - Plot Body Surface Pressure - # + for b in 1:Int(bvp.nbodies[]) + @series begin + label --> false + seriescolor --> b + linewidth = 2 + return bvp.controlpoint[ + 1, Int(bvp.endpanelidxs[1, b]):Int(bvp.endpanelidxs[2, b]) + ], + bout.Vtan_out[Int(bvp.endpanelidxs[1, b]):Int(bvp.endpanelidxs[2, b])] ./ Vref + end + end +end + +#---------------------------------# +# SET XTICKS AT LE, TE, ETC. # +#---------------------------------# +function add_ticks(xt, xl, tmp, tol, tickdigits=2) + if !any(abs.(tmp .- xt) .< tol / 10) + if !any(abs.(tmp .- xt) .< tol) + push!(xl, @sprintf "%3.*f" tickdigits tmp) + else + push!(xl, "") + end + push!(xt, tmp) + end + return xt, xl +end + +function determine_geometry_xlabels(bvp, rsp, wvp; tol=1e-2, tickdigits=2) + xt = [] + xl = [] + + # Rotor + for r in 1:Int(rsp.nbodies[]) + push!(xt, rsp.node[1, Int(rsp.endnodeidxs[r, r])]) + push!(xl, @sprintf "%3.*f" tickdigits rsp.node[1, Int(rsp.endnodeidxs[r, r])]) + end + + # Centerbody + for i in 1:2 + tmp = bvp.node[1, Int(bvp.endnodeidxs[i, 2])] + xt, xl = add_ticks(xt, xl, tmp, tol, tickdigits) + end + + # Duct + for f in [maximum, minimum] + tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) + xt, xl = add_ticks(xt, xl, tmp, tol, tickdigits) + end + + # Wake + tmp = maximum(wvp.node[1, Int(wvp.endnodeidxs[1, 1]):Int(wvp.endnodeidxs[2, 1])]) + xt, xl = add_ticks(xt, xl, tmp, tol, tickdigits) + + return (xt, xl) +end + +function determine_geometry_ylabels(bvp, rsp; tol=1e-2, tickdigits=2) + yt = [] + yl = [] + + # Rotor + for r in 1:Int(rsp.nbodies[]) + push!(yt, rsp.node[2, Int(rsp.endnodeidxs[r, r])]) + push!(yl, @sprintf "%3.*f" tickdigits rsp.node[2, Int(rsp.endnodeidxs[r, r])]) + end + + # Centerbody + for i in 1:2 + tmp = bvp.node[2, Int(bvp.endnodeidxs[i, 2])] + yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits) + end + + # Duct front/back + for f in [findmax, findmin] + tmp = f(bvp.node[1, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) + tmpy = bvp.node[2, tmp[2]] + yt, yl = add_ticks(yt, yl, tmpy, tol, tickdigits) + end + + # Duct inner/outer + for f in [maximum, minimum] + tmp = f(bvp.node[2, Int(bvp.endnodeidxs[1, 1]):Int(bvp.endnodeidxs[2, 1])]) + yt, yl = add_ticks(yt, yl, tmp, tol, tickdigits) + end + + return (yt, yl) +end + +#---------------------------------# +#PLOT GEOMETRY UNDER DISTRIBUTIONS# +#---------------------------------# +@recipe function plot_underlayed_geometry( + ::underlayGeometry, + bvp, + rsp; + 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 + ), + discrete_labels=true, +) + color_palette --> [ + default_colors.primary, + default_colors.secondary, + default_colors.tertiary, + default_colors.quaternary, + default_colors.quinary, + default_colors.plotsgray, + ] + + # Plot specific values + axis --> false + ticks --> false + aspect_ratio --> 1 + + # - Plot Body Geometry - # + for b in 1:Int(bvp.nbodies[]) + @series begin + label --> false + seriescolor --> b + linealpha --> 1 / 2 + linewidth --> 0.5 + return bvp.node[1, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])], + bvp.node[2, Int(bvp.endnodeidxs[1, b]):Int(bvp.endnodeidxs[2, b])] + end + end + + # - Close Center Body - # + @series begin + label --> false + seriescolor --> Int(bvp.nbodies[]) + linealpha --> 1 / 2 + linewidth --> 0.5 + return ones(2) * bvp.node[1, Int(bvp.endnodeidxs[2, 2])], + [0.0; bvp.node[2, Int(bvp.endnodeidxs[2, 2])]] + end + + # - Plot Rotor Geometry - # + for r in 1:Int(rsp.nbodies[]) + @series begin + z_order --> 1 + label --> false + seriescolor --> 3 + linewidth --> 1.5 + return rsp.node[1, Int(rsp.endnodeidxs[1, r]):Int(rsp.endnodeidxs[2, r])], + rsp.node[2, Int(rsp.endnodeidxs[1, r]):Int(rsp.endnodeidxs[2, r])] + end + end +end + +#---------------------------------# +# PLOT STAGNATION # +#---------------------------------# +@recipe function plot_stagnation( + ::plotStagnation, + blo, + bvp, + 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 + ), +) + @series begin + seriestype --> :scatter + color --> 2 + markerstrokealpha --> 0 + label --> "" + + stagz = sum( + bvp.controlpoint[1, blo.stagnation_indices] .* + [1.0 - blo.split_ratio; blo.split_ratio], + ) + stagr = sum( + bvp.controlpoint[2, blo.stagnation_indices] .* + [1.0 - blo.split_ratio; blo.split_ratio], + ) + return [stagz], [stagr] + end + + return nothing +end + +#---------------------------------# +# PLOT δ₂ # +#---------------------------------# +@recipe function plot_momentum_thickness( + ::plotMomentum, + blo, + bvp; + 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, + bl_ylim=nothing, +) + 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 + + # stagnation point + stagz = sum( + bvp.controlpoint[1, blo.stagnation_indices] .* + [1.0 - blo.split_ratio; blo.split_ratio], + ) + stagr = sum( + bvp.controlpoint[2, blo.stagnation_indices] .* + [1.0 - blo.split_ratio; blo.split_ratio], + ) + + # stagnation point normal + stagnz = sum( + bvp.normal[1, blo.stagnation_indices] .* [1.0 - blo.split_ratio; blo.split_ratio] + ) + stagnr = sum( + bvp.normal[2, blo.stagnation_indices] .* [1.0 - blo.split_ratio; blo.split_ratio] + ) + + if !isnothing(bl_ylim) + ylim --> bl_ylim + end + + @series begin + color --> 1 + label --> false + # linewidth --> 0.5 + if isapprox(blo.separation_point_ratio_upper, 1.0) + arrow --> :closed + end + + # spline the states vs steps + d2_upper = blo.upper_solved_states[1, :] + + # spline the control points vs surface length + cpzu = akima_smooth( + blo.surface_length_upper, + [stagz; bvp.controlpoint[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]], + blo.upper_solved_steps, + ) + cpru = akima_smooth( + blo.surface_length_upper, + [stagr; bvp.controlpoint[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]], + blo.upper_solved_steps, + ) + + nhatzu = akima_smooth( + blo.surface_length_upper, + [stagnz; bvp.normal[1, blo.stagnation_indices[2]:Int(bvp.npanel[1])]], + blo.upper_solved_steps, + ) + nhatru = akima_smooth( + blo.surface_length_upper, + [stagnr; bvp.normal[2, blo.stagnation_indices[2]:Int(bvp.npanel[1])]], + blo.upper_solved_steps, + ) + + # thickness point is thickness normal to control points + d2zu = cpzu .+ (d2_upper .* nhatzu) .* scale_thickness + d2ru = cpru .+ (d2_upper .* nhatru) .* scale_thickness + + return d2zu, d2ru + end + + @series begin + color --> 2 + label --> false + # linewidth --> 0.5 + if isapprox(blo.separation_point_ratio_lower, 1.0) + arrow --> :closed + end + + # spline the states vs steps + d2_lower = blo.lower_solved_states[1, :] + + # spline the control points vs surface length + cpzl = akima_smooth( + blo.surface_length_lower, + [bvp.controlpoint[1, blo.stagnation_indices[2]:-1:1]; stagz], + blo.lower_solved_steps, + ) + cprl = akima_smooth( + blo.surface_length_lower, + [bvp.controlpoint[2, blo.stagnation_indices[2]:-1:1]; stagr], + blo.lower_solved_steps, + ) + + nhatzl = akima_smooth( + blo.surface_length_lower, + [bvp.normal[1, blo.stagnation_indices[2]:-1:1]; stagnz], + blo.lower_solved_steps, + ) + nhatrl = akima_smooth( + blo.surface_length_lower, + [bvp.normal[2, blo.stagnation_indices[2]:-1:1]; stagnr], + blo.lower_solved_steps, + ) + + # thickness point is thickness normal to control points + d2zl = cpzl .+ (d2_lower .* nhatzl) .* scale_thickness + d2rl = cprl .+ (d2_lower .* nhatrl) .* scale_thickness + + return d2zl, d2rl + end + + return nothing +end + +#---------------------------------# +# PLOT STREAMLINES # +#---------------------------------# +@recipe function plot_streamlines( + ::plotStreamlines, + bvp, + bvs, + wvp, + wvs, + rsp, + rss, + vinf; + starting_radial_points=range(0.001, 1.0; length=20), + axial_range=[-0.2, 1.0], + nominal_step_size=1e-2, + step_limit=Int(1e2), + dot_tol=0.999, + stag_tol=0.4, + 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 + + p, _ = DuctAPE.calculate_streamlines( + bvp, + bvs, + wvp, + wvs, + rsp, + rss, + vinf; + starting_radial_points=starting_radial_points, + axial_range=axial_range, + step_limit=step_limit, + nominal_step_size=nominal_step_size, + integration_options=integration_options, + stag_tol=stag_tol, + ) + + for (z, r) in zip(eachcol(p[1, :, :]), eachcol(p[2, :, :])) + @series begin + color --> default_colors.plotsgray + label --> "" + linewidth --> 0.5 + return z, r + end + end + + return nothing +end diff --git a/test/boundary_layer_tests.jl b/test/boundary_layer_tests.jl new file mode 100644 index 00000000..8e204162 --- /dev/null +++ b/test/boundary_layer_tests.jl @@ -0,0 +1,129 @@ +@testset "Boundary Layer Functions" begin + ##### ----- Head's Method Tests ----- ##### + + # @testset "Step Setter" begin + N = 10 + first_step_size = 1 + total_length = 100 + steps = dt.set_boundary_layer_steps(N, first_step_size, total_length) + + @test steps[1] == first_step_size + @test steps[end] == total_length + @test length(steps) == N + + # @testset "Raw Arc Length Calcuations" begin + panel_lengths = ones(10) + s = dt.arc_lengths_from_panel_lengths(panel_lengths) + + @test s == 0:9 +end + +# these are for the green's method which may get resurrected at some point... +# @testset "Boundary Layer Functions" begin +# # - Stagnation Point - # +# cp_duct = abs.(range(-10, 10; step=1)) .+ 1.1 +# panel_lengths = ones(length(cp_duct)) +# s, s_stagnation, lower_length, upper_length = dt.split_at_stagnation_point( +# panel_lengths, cp_duct +# ) +# @test s == cumsum([0.0; panel_lengths[1:(end - 1)]]) +# @test s_stagnation == 10.0 +# @test lower_length == 10.0 +# @test upper_length == 10.0 + +# # @testset "Radius of Curvature Calculation" begin +# s = range(1, 10) +# controlpoint = [s'; s' .^ 2] +# ss = 4.0 +# R = dt.calculate_radius_of_curvature(s, controlpoint, ss) +# @test isapprox(R, 1.0 / 0.003816453371976) + +# # @testset "Boundary Layer Initialization" begin +# x = 1.0 +# Rex = 1.0 +# @test 0.036 == dt.d2_init(x, Rex) + +# @test dt.calculate_H12bar0(2.0, 0.0) == 1.0 / (1.0 - 6.55) + +# @test dt.calculate_CEeq(1.0, 0.0, 1.0, 1.0) == sqrt((1.0 - 0.32) / 1.2 + 0.0001) - 0.01 + +# s_init = 0.1 +# r_init = 1.0 +# Ue = 1.0 +# M = 0.1 +# rhoe = 1.0 +# mue = 1e-7 + +# states, Cf_init, H12_init = dt.initialize_turbulent_boundary_layer_states( +# s_init, r_init, Ue, M, rhoe, mue +# ) + +# @test isapprox( +# states, [0.00022714464401286951, 1.3836007395668004, 0.10980454099439806] +# ) +# @test isapprox(Cf_init, 0.003581879845103732) +# @test isapprox(H12_init, 1.3883679410459338) + +# # @testset "Schlichting" begin +# @test dt.d2_init(1.0, 1.0) == 0.036 + +# # - H12bar_init - # +# @test dt.H12bar_init(1.0, 1.0) == dt.calculate_H12bar0(1.0, 1.0) + +# # - CE_init - # +# @test dt.CE_init(1.0, 0.1, 0.1) == dt.calculate_CEeq(1.0, 0.1, 1, 0.1) + +# # - Fc - # +# @test dt.Fc(0) == 1.0 + +# # - FR - # +# @test dt.FR(0) == 1.0 + +# # - Reynolds - # +# @test dt.calculate_Re(rhoe, Ue, 1.0, mue) == rhoe * Ue / mue + +# # - H12bar0 - # +# @test dt.calculate_H12bar0(0, 1.0) == 1.0 + +# # - Cf - # +# @test isapprox(dt.calculate_Cf(2.2, 1.0, 1.0), 0.0, atol=eps()) + +# # - H12 - # +# @test dt.calculate_H12(0.0, sqrt(5)) == 1.0 + +# # - Ctau - # +# @test dt.calculate_Ctau(1.0, 0.0, sqrt(10)) == 2.0 * (0.024 + 1.2) + +# # - F - # +# @test dt.calculate_F(1.0, 1.0) == (0.02 + 1.0 + 0.8 * 1.0 / 3) / 1.01 + +# # - H1 - # +# @test dt.calculate_H1(2.0) == 3.15 + 1.72 - 0.01 + +# # - dH12bardH1 - # +# @test dt.calculate_dH12bardH1(2.0) == -1.0 / (1.72 + 0.02) + +# # - Ri - # +# @test dt.calculate_richardson_number(1.0, 1.0, 0.0, 1.0, 1.0) == 2.0 / 3.0 * (1.0 + 0.3) + +# # - Secondary Influences - # +# Ri = 1.0 +# @test dt.longitudinal_curvature_influence(M, Ri) == 8.014 +# @test dt.lateral_strain_influence(ones(6)...) == -8.0 - 1 / 3 +# @test dt.dilation_influence(ones(7)...) == 10 + 1 / 3 + +# # - d2dUedsUeeq0 - # +# @test dt.calculate_d2dUedsUeeq0(ones(4)...) == 0.625 + +# # - CEeq0 - # +# @test dt.calculate_CEeq0(ones(4)...) == 0.0 + +# # - Ctaueq0 - # +# @test isapprox(dt.calculate_Ctaueq0(ones(3)...), 1.936) + +# # - CEeq - # +# @test isapprox(dt.calculate_CEeq(ones(4)...), 0.6907204085147591) + +# # - d2dUedsUeeq - # +# @test dt.calculate_d2dUedsUeeq(ones(4)...) == -0.25 +# end diff --git a/test/data/basic_two_rotor_for_test_NEW.jl b/test/data/basic_two_rotor_for_test_NEW.jl index 650cb905..d77bac35 100644 --- a/test/data/basic_two_rotor_for_test_NEW.jl +++ b/test/data/basic_two_rotor_for_test_NEW.jl @@ -35,7 +35,7 @@ muinf = [1.78e-5] asound = [340.0] Omega = [5000.0, 0.0] * pi / 30 # convert from RPM to rad/s -operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) +operating_point = dt.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) Vref = [10.0] Rref = [Rtip] diff --git a/test/data/dfdc_csor/bin/buildandrun.sh b/test/data/dfdc_csor/bin/buildandrun.sh old mode 100755 new mode 100644 diff --git a/test/data/dfdc_csor/bin/test b/test/data/dfdc_csor/bin/test old mode 100755 new mode 100644 diff --git a/test/data/dfdc_csor/src-driver/hello b/test/data/dfdc_csor/src-driver/hello old mode 100755 new mode 100644 diff --git a/test/data/dfdc_init_iter1/ductape_parameters.jl b/test/data/dfdc_init_iter1/ductape_parameters.jl index 363a7cdd..5bff8f91 100644 --- a/test/data/dfdc_init_iter1/ductape_parameters.jl +++ b/test/data/dfdc_init_iter1/ductape_parameters.jl @@ -85,7 +85,7 @@ reference_parameters = dt.ReferenceParameters(Vref, Rtip) D = 2.0 * Rtip # rotor diameter n = RPM / 60.0 # rotation rate in revolutions per second Vinf = 1.0 * n * D -operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) +operating_point = dt.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) # ##### ----- Solid Body Coordinates (GEOM in the DFDC case file) ----- ##### # # copied from the DFDC case file and put here at end since they're so long diff --git a/test/data/dfdc_panel_integration/bin/test b/test/data/dfdc_panel_integration/bin/test old mode 100755 new mode 100644 diff --git a/test/data/dfdc_panel_integration/src-driver/hello b/test/data/dfdc_panel_integration/src-driver/hello old mode 100755 new mode 100644 diff --git a/test/figures/gengifs.sh b/test/figures/gengifs.sh old mode 100755 new mode 100644 diff --git a/test/iteration_step_tests.jl b/test/iteration_step_tests.jl index 274e4cbb..5c1dd162 100644 --- a/test/iteration_step_tests.jl +++ b/test/iteration_step_tests.jl @@ -20,12 +20,9 @@ println("\nITERATION STEP THROUGH TESTS") include(datapath * "ductape_formatted_dfdc_geometry.jl") include(datapath * "ductape_parameters.jl") #NOTE, for some reason, julia doesn't recognize that this was defined in the above file... - operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) + operating_point = dt.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) ducted_rotor = dt.DuctedRotor( - rp_duct_coordinates, - rp_centerbody_coordinates, - rotor, - paneling_constants, + rp_duct_coordinates, rp_centerbody_coordinates, rotor, paneling_constants ) options = dt.DFDC_options(; @@ -96,7 +93,9 @@ println("\nITERATION STEP THROUGH TESTS") # copy over operating point for f in fieldnames(typeof(operating_point)) - solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + if f != :units + solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + end end # generate inputs diff --git a/test/manual_tests/dfdc_comparisons/convert.sh b/test/manual_tests/dfdc_comparisons/convert.sh old mode 100755 new mode 100644 diff --git a/test/pre_processing_tests.jl b/test/pre_processing_tests.jl index ed5d3b59..f50ff20f 100644 --- a/test/pre_processing_tests.jl +++ b/test/pre_processing_tests.jl @@ -66,21 +66,13 @@ println("\nPRECOMPUTED ROTOR & WAKE INPUTS") rpb4 = copy(rp_duct_coordinates) - dt.place_duct!( - rp_duct_coordinates, - rotor.Rtip[1], - rotor.rotorzloc[1], - rotor.tip_gap[1], - ) + dt.place_duct!(rp_duct_coordinates, rotor.Rtip[1], rotor.rotorzloc[1], rotor.tip_gap[1]) @test rp_duct_coordinates[1, :] == rpb4[1, :] @test rp_duct_coordinates[2, :] == rpb4[2, :] .- 0.75 Rtips, Rhubs = dt.get_blade_ends_from_body_geometry( - rp_duct_coordinates, - rp_centerbody_coordinates, - rotor.tip_gap, - rotor.rotorzloc, + rp_duct_coordinates, rp_centerbody_coordinates, rotor.tip_gap, rotor.rotorzloc ) @test all(Rtips .== 1.0) @@ -170,11 +162,7 @@ println("\nPRECOMPUTED ROTOR & WAKE INPUTS") # rotor blade element objects blade_elements, airfoils = dt.interpolate_blade_elements( - rotor, - Rtips, - Rhubs, - rotor_source_panels.controlpoint[2, :], - problem_dimensions.nbe, + rotor, Rtips, Rhubs, rotor_source_panels.controlpoint[2, :], problem_dimensions.nbe ) @test blade_elements.inner_fraction == [0.75 0.75; 0.25 0.25] @@ -350,58 +338,55 @@ end end @testset "Rotor/Wake Aero Initialization" begin -r1 = [0.25; 0.5; 0.75; 1.0] -Rtip = [1.0, 1.0] -rnondim1 = r1 ./ Rtip[1] -rnondim = [rnondim1 rnondim1] -afparams1 = dt.c4b.DFDCairfoil() -rotorzloc = [0.25, 0.75] -r = rnondim -chords = 0.1 * ones(size(rnondim)) -twists = 20.0 * pi / 180.0 * ones(size(rnondim)) -airfoils = fill(afparams1, 4, 2) -Rhub = [0.25, 0.25] -Rtip = Rtip -tip_gap = [0.0, 0.0] -B = [2, 4] -fliplift = [0.0, 0.0] - -rotor = dt.Rotor( - B, rotorzloc, r, Rhub, Rtip, chords, twists, tip_gap, airfoils, fliplift -) - -ncenterbody_inlet = 1 -nduct_inlet = 1 -nwake_sheets = 3 -wake_length = 1.0 -npanels = [2, 1, 4] -dte_minus_cbte = 0 - -paneling_constants = dt.PanelingConstants( - nduct_inlet, ncenterbody_inlet, npanels, dte_minus_cbte, nwake_sheets, wake_length -) - -Vinf = [10.0] -rhoinf = [1.226] -muinf = [1.78e-5] -asound = [340.0] -Omega = [5000.0, 0.0] * pi / 30 # convert from RPM to rad/s - -operating_point = dt.OperatingPoint(Vinf, rhoinf, muinf, asound, Omega) - -Vref = [10.0] -Rref = [Rtip] -reference_parameters = dt.ReferenceParameters(Vref, Rref) - -duct_coordinates = [1.0 2.0; 0.5 1.5; 0.0 2.0; 0.5 2.5; 1.0 2.0] -centerbody_coordinates = [0.0 0.0; 0.5 0.5; 1.0 0.0] - -ducted_rotor = dt.DuctedRotor( - duct_coordinates, - centerbody_coordinates, - rotor, - paneling_constants, -) + r1 = [0.25; 0.5; 0.75; 1.0] + Rtip = [1.0, 1.0] + rnondim1 = r1 ./ Rtip[1] + rnondim = [rnondim1 rnondim1] + afparams1 = dt.c4b.DFDCairfoil() + rotorzloc = [0.25, 0.75] + r = rnondim + chords = 0.1 * ones(size(rnondim)) + twists = 20.0 * pi / 180.0 * ones(size(rnondim)) + airfoils = fill(afparams1, 4, 2) + Rhub = [0.25, 0.25] + Rtip = Rtip + tip_gap = [0.0, 0.0] + B = [2, 4] + fliplift = [0.0, 0.0] + + rotor = dt.Rotor( + B, rotorzloc, r, Rhub, Rtip, chords, twists, tip_gap, airfoils, fliplift + ) + + ncenterbody_inlet = 1 + nduct_inlet = 1 + nwake_sheets = 3 + wake_length = 1.0 + npanels = [2, 1, 4] + dte_minus_cbte = 0 + + paneling_constants = dt.PanelingConstants( + nduct_inlet, ncenterbody_inlet, npanels, dte_minus_cbte, nwake_sheets, wake_length + ) + + Vinf = [10.0] + rhoinf = [1.226] + muinf = [1.78e-5] + asound = [340.0] + Omega = [5000.0, 0.0] * pi / 30 # convert from RPM to rad/s + + operating_point = dt.OperatingPoint(Vinf, Omega, rhoinf, muinf, asound) + + Vref = [10.0] + Rref = [Rtip] + reference_parameters = dt.ReferenceParameters(Vref, Rref) + + duct_coordinates = [1.0 2.0; 0.5 1.5; 0.0 2.0; 0.5 2.5; 1.0 2.0] + centerbody_coordinates = [0.0 0.0; 0.5 0.5; 1.0 0.0] + + ducted_rotor = dt.DuctedRotor( + duct_coordinates, centerbody_coordinates, rotor, paneling_constants + ) options = dt.set_options() @@ -454,7 +439,9 @@ ducted_rotor = dt.DuctedRotor( # copy over operating point for f in fieldnames(typeof(operating_point)) - solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + if f != :units + solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + end end # - Do preprocessutations - # @@ -566,7 +553,9 @@ ducted_rotor = dt.DuctedRotor( # copy over operating point for f in fieldnames(typeof(operating_point)) + if f != :units solve_parameter_tuple.operating_point[f] .= getfield(operating_point, f) + end end # - Do preprocessutations - # diff --git a/test/runtests.jl b/test/runtests.jl index 10af52b6..75114301 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -39,7 +39,8 @@ include("wake_aero_tests.jl") # - post process related tests - # include("post_processing_tests.jl") - +include("thermodynamics_tests.jl") +include("boundary_layer_tests.jl") ########################################################## #EVERYTHING BELOW THIS POINT NEEDS TO BE UPDATED @@ -48,7 +49,6 @@ include("post_processing_tests.jl") # - Active Development - # # include("solve_checks.jl") - # - Need to update, add, fix, etc. - # # include("body_aero_tests.jl") # # TODO: need to add source and body induced velocities to induced velocity on rotor test diff --git a/test/thermodynamics_tests.jl b/test/thermodynamics_tests.jl new file mode 100644 index 00000000..de4e36ff --- /dev/null +++ b/test/thermodynamics_tests.jl @@ -0,0 +1,35 @@ +@testset "Thermodynamic Relations" begin + @testset "Standard Atmosphere" begin + + # - SI Units - # + T, P, rho, mu = dt.standard_atmosphere(0.0) + @test isapprox(T, 273.15 + 14.99, atol=1e-12) + @test isapprox(P, 101400.93090454886, atol=1e-12) + @test isapprox(rho, 1.22661378741097, atol=1e-12) + @test isapprox(mu, 1.7889517587366847e-5, atol=1e-12) + + # - Imperial Units - # + T, P, rho, mu = dt.standard_atmosphere(dt.Imperial(), 0.0) + @test isapprox(T, 58.982, atol=1e-12) + @test isapprox(P, 2117.8024776319244, atol=1e-12) + @test isapprox(rho, 0.002380023263989253, atol=1e-12) + @test isapprox(mu, 3.7363034244069305e-7, atol=1e-12) + end + + @testset "OperatingPoint" begin + + # - SI Units - # + si = dt.OperatingPoint(1.0, 1.0) + @test isapprox(si.Ttot[], 288.14049793357566, atol=1e-12) + @test isapprox(si.Ptot[], 101401.54421276736, atol=1e-12) + @test isapprox(si.asound[], 340.1974608958744, atol=1e-12) + @test isapprox(si.Minf[], si.Vinf[] / si.asound[], atol=1e-12) + + # - Imperial Units - # + imp = dt.OperatingPoint(dt.Imperial(), 1.0, 1.0) + @test isapprox(imp.Ttot[], 58.98200946928549, atol=1e-12) + @test isapprox(imp.Ptot[], 2117.8036676437946, atol=1e-12) + @test isapprox(imp.asound[], 1116.133498437807, atol=1e-12) + @test isapprox(imp.Minf[], imp.Vinf[] / imp.asound[], atol=1e-12) + end +end