The following is a basic tutorial on how to set up and run an analysis of a ducted fan in DuctAPE.
We begin by loading the package:
using DuctAPE
The next step is to create the input object of type DuctedRotor
.
DuctedRotor(duct_coordinates, centerbody_coordinates, rotor, paneling_constants)
Arguments
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.
sourceWe begin by defining a matrix of coordinates for the duct and another for the centerbody geometries. For example:
duct_coordinates = [
+ 0.304466 0.158439
+ 0.294972 0.158441
+ 0.28113 0.158423
+ 0.266505 0.158365
+ 0.251898 0.158254
+ 0.237332 0.158088
+ 0.222751 0.157864
+ 0.208123 0.157586
+ 0.193399 0.157258
+ 0.178507 0.156897
+ 0.16349 0.156523
+ 0.148679 0.156177
+ 0.134222 0.155902
+ 0.12 0.155721
+ 0.106044 0.155585
+ 0.092531 0.155498
+ 0.079836 0.155546
+ 0.067995 0.155792
+ 0.057025 0.156294
+ 0.046983 0.157103
+ 0.037937 0.158256
+ 0.029956 0.159771
+ 0.02311 0.161648
+ 0.017419 0.163862
+ 0.012842 0.166404
+ 0.009324 0.169289
+ 0.006854 0.172546
+ 0.005484 0.176154
+ 0.005242 0.180005
+ 0.006112 0.184067
+ 0.00809 0.188086
+ 0.011135 0.192004
+ 0.015227 0.19579
+ 0.020339 0.199393
+ 0.026403 0.202735
+ 0.033312 0.205736
+ 0.040949 0.208332
+ 0.049193 0.210487
+ 0.057935 0.212174
+ 0.067113 0.21339
+ 0.076647 0.214136
+ 0.086499 0.214421
+ 0.09661 0.214255
+ 0.10695 0.213649
+ 0.117508 0.212618
+ 0.12838 0.211153
+ 0.139859 0.209267
+ 0.151644 0.207051
+ 0.163586 0.204547
+ 0.175647 0.201771
+ 0.187807 0.198746
+ 0.20002 0.19549
+ 0.212269 0.192017
+ 0.224549 0.188335
+ 0.236794 0.18447
+ 0.249026 0.180416
+ 0.261206 0.176188
+ 0.273301 0.171796
+ 0.28524 0.16727
+ 0.29644 0.162842
+ 0.304542 0.159526
+]
centerbody_coordinates = [
+ 0.0 0.0
+ 0.000586 0.005293
+ 0.002179 0.010047
+ 0.004736 0.014551
+ 0.008231 0.018825
+ 0.012632 0.022848
+ 0.01788 0.026585
+ 0.023901 0.030001
+ 0.030604 0.033068
+ 0.0379 0.035771
+ 0.045705 0.038107
+ 0.053933 0.040075
+ 0.06254 0.04169
+ 0.071451 0.042966
+ 0.08063 0.043916
+ 0.090039 0.044561
+ 0.09968 0.044922
+ 0.109361 0.044999
+ 0.12 0.044952
+ 0.135773 0.04495
+ 0.151899 0.04493
+ 0.16806 0.044913
+ 0.184232 0.044898
+ 0.200407 0.044882
+ 0.21658 0.044866
+ 0.232723 0.044847
+ 0.248578 0.044839
+ 0.262095 0.044564
+ 0.274184 0.043576
+ 0.285768 0.041795
+ 0.296701 0.039168
+ 0.306379 0.035928
+]
The body geometry coordinates must be input as columns of z (axial) and r (radial) coordinates, in that order.
The next step is to assemble an object of type Rotor
which contains the geometric information required to define the rotor(s) and their respective blade elements.
Rotor(
+ B, rotorzloc, r, Rhub, Rtip, chords, twists, tip_gap, airfoils, fliplift
+)
Composite type containing the rotor(s) geometric properties.
Note that the actual struct requires the inputs to be arrays, but there is a constructor function that will take in scalars and automatically build the array-based struct.
Arguments
B::AbstractVector{Float}
: The number of blades for each rotor. May not be an integer, but usually is.rotorzloc::AbstractVector{Float}
: Dimensional, axial position of each rotor.r::AbstractArray{Float}
: Non-dimensional radial locations of each blade element.Rhub::AbstractVector{Float}
: Dimensional hub radius of rotor. (may be changed if it does not match the radial position of the centerbody geometry at the selected rotorzloc
.Rtip::AbstractVector{Float}
: Dimensional tip radius of rotor. Is used to determine the radial position of the duct if the autoshiftduct
option is selected.chords::AbstractArray{Float}
: Dimensional chord lengths of the blade elements.twists::AbstractArray{Float}
: Blade element angles, in radians.tip_gap::AbstractVector{Float}
: Currently unused, do not set to anything other than zeros.airfoils::AbstractArray{AFType}
: Airfoil types describing the airfoil polars for each blade element. Currently only fully tested with C4Blade.DFDCairfoil
types.fliplift::AbstractVector{Bool}
: Flag to indicate if the airfoil lift values should be flipped or not.
sourceIn this example, we have a single rotor defined as follows.
# number of rotors
+B = 5
+
+# rotor axial location
+rotorzloc = 0.12
+
+# rotor tip radius
+Rtip = 0.15572081487373543
+
+# rotor hub radius
+Rhub = 0.04495252299071941
+
+# non-dimensional blade element radial stations
+r = [
+ 0.050491
+ 0.061567
+ 0.072644
+ 0.083721
+ 0.094798
+ 0.10587
+ 0.11695
+ 0.12803
+ 0.13911
+ 0.15018
+] ./ Rtip
+
+# dimensional chord lengths
+chords = [
+ 0.089142
+ 0.079785
+ 0.0713
+ 0.063979
+ 0.057777
+ 0.052541
+ 0.048103
+ 0.044316
+ 0.041061
+ 0.038243
+]
+
+# twist angles (from plane of rotation) in radians
+twists = [
+ 69.012
+ 59.142
+ 51.825
+ 46.272
+ 41.952
+ 38.509
+ 35.699
+ 33.354
+ 31.349
+ 29.596
+] .* pi / 180.0
+
+# DFDC-type airfoil object
+afparams = DuctAPE.c4b.DFDCairfoil(;
+ alpha0=0.0,
+ clmax=1.5,
+ clmin=-1.0,
+ dclda=6.28,
+ dclda_stall=0.5,
+ dcl_stall=0.2,
+ cdmin=0.012,
+ clcdmin=0.1,
+ dcddcl2=0.005,
+ cmcon=0.0,
+ Re_ref=2e5,
+ Re_exp=0.35,
+ mcrit=0.7,
+)
+
+# all airfoils are the same
+airfoils = fill(afparams, length(r)) # specify the airfoil array
+
+# assemble rotor parameters
+rotor = DuctAPE.Rotor(
+ [B],
+ [rotorzloc],
+ r,
+ [Rhub],
+ [Rtip],
+ chords,
+ twists,
+ [0.0], # currently only zero tip gaps work.
+ airfoils,
+ [0.0], # can flip the cl lookups on the fly if desired, say, for stator sections
+)
Airfoil types for DuctAPE are currently contained in the C4Blade (Cascade Compatible CCBlade) sub-module of DuctAPE which is exported as c4b
and also contains the various airfoil evaluation functions used for the blade element lookups. The available airfoil types include all the airfoil types from CCBlade, as well as DFDCairfoil
which is an XROTOR-like parametric cascade polar used in DFDC. In addition there are untested cascade types with similar structure to CCBlades airfoil types called DTCascade
. Furthermore, there is an experimental actuator disk model implemented via the ADM
airfoil type in C4Blade.
The PanelingConstants
object contains the constants required for DuctAPE to re-panel the provided geometry into a format compatible with the solve structure. Specifically, the DuctAPE solver makes some assumptions on the relative positioning of the body surfaces relative to the wakes and each other; and this is most easily guarenteed by a re-paneling of the provided body surface geometry. The PanelingConstants
object is also used to build all of the preallocated caches inside DuctAPE, which can be done up-front if desired. Note that there is some functionality in place for cases when the user wants to keep their own specified geometry, but this functionality should be used with caution and only by users who are certain their provided geometry is in the compatible format. See the Examples for an example.
PanelingConstants(
+ nduct_inlet,
+ ncenterbody_inlet,
+ npanels,
+ dte_minus_cbte,
+ nwake_sheets,
+ wake_length=1.0,
+)
Constants used in re-paneling geometry.
Note that unlike other input structures, this one, in general, does not define fields as vectors. This is because these values should not change throughout an optimization, even if the geometry may change. Otherwise, discontinuities could be experienced.
Arguments
nduct_inlet::Int
: The number of panels to use for the duct inlet (this number is used for both the casing and nacelle re-paneling)ncenterbody_inlet::Int
: The number of panels to use for the centerbody inlet.npanels::AbstractVector{Int}
: A vector containing the number of panels between discrete locations inside the wake. Specifically, the number of panels between the rotors, between the last rotor and the first body trailing edge, between the body trailing edges (if different), and between the last body trailing edge and the end of the wake. The length of this vector should be N+1 (where N is the number of rotors) if the duct and centerbody trailing edges are aligned, and N+2 if not.dte_minus_cbte::Float
: An indicator concerning the hub and duct trailing edge relative locations. Should be set to -1 if the duct trailing edge axial position minus the centerbody trailing edge axial position is negative, +1 if positive (though any positive or negative number will suffice), and zero if the trailing edges are aligned.nwake_sheets::Int
: The number of wake sheets to use. Note this will also be setting the number of blade elements to use.wake_length::Float=1.0
: Non-dimensional (based on the length from the foremost body leading edge and the aftmost body trailing edge) length of the wake extending behind the aftmost body trailing edge.
source# number of panels for the duct inlet
+nduct_inlet = 30
+
+# number of panels for the center body inlet
+ncenterbody_inlet = 30
+
+# number of panels from:
+# - rotor to duct trailing edge
+# - duct trailing edge to center body trailing edge
+# - center body trailing edge to end of wake
+npanels = [30, 1, 30]
+
+# the duct trailing edge is ahead of the centerbody trailing edge.
+dte_minus_cbte = -1.0
+
+# number of wake sheets (one more than blade elements to use)
+nwake_sheets = 11
+
+# non-dimensional wake length aft of rear-most trailing edge
+wake_length = 0.8
+
+# assemble paneling constants
+paneling_constants = DuctAPE.PanelingConstants(
+ nduct_inlet, ncenterbody_inlet, npanels, dte_minus_cbte, nwake_sheets, wake_length
+)
We are now posed to construct the DuctedRotor
input type.
# assemble ducted_rotor object
+ducted_rotor = DuctAPE.DuctedRotor(
+ duct_coordinates,
+ centerbody_coordinates,
+ rotor,
+ paneling_constants,
+)
Next we will assemble the operating point which contains information about the freestream as well as the rotor rotation rate(s).
OperatingPoint(Vinf, rhoinf, muinf, asound, Omega)
DuctedRotor operating point information.
Arguments
Vinf::AbstractVector{Float}
: Freestream velocity magnitude (which is only in the axial direction).rhoinf::AbstractVector{Float}
: Freestream densitymuinf::AbstractVector{Float}
: Freestream viscosityasound::AbstractVector{Float}
: Freestream speed of soundOmega::AbstractVector{Float}
: Rotor rototation rate(s)
source# Freestream
+Vinf = 0.0 # hover condition
+rhoinf = 1.226
+asound = 340.0
+muinf = 1.78e-5
+
+# Rotation Rate
+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)
The reference parameters are used in the post-processing non-dimensionalizations.
ReferenceParameters(Vref, Rref)
Reference parameters for post-process non-dimensionalization.
Note that the actual struct requires the inputs to be arrays, but there is a constructor function that will take in scalars and automatically build the array-based struct.
Arguments
Vref::AbstractVector{Float}
: Reference velocity.Rref::AbstractVector{Float}
: Reference rotor tip radius.
source# reference velocity (close to average axial velocity at rotor in this case)
+Vref = 50.0
+
+# reference radius (usually tip radius of rotor)
+Rref = Rtip
+
+# assemble reference parameters
+reference_parameters = DuctAPE.ReferenceParameters([Vref], [Rref])
The default options should be sufficient for just starting out and are set through the set_options
function.
set_options(; kwargs...)
+set_options(multipoint; kwargs...)
Set the options for DuctAPE to use.
Note that the vast majority of the available options are defined through keyword arguments. See the documentation for the various option types for more information.
Arguments
multipoint::AbstractArray{OperatingPoint}
: a vector of operating points to use if running a multi-point analysis.
sourceoptions = DuctAPE.set_options()
Options{Bool, Vector{Bool}, Float64, Vector{Int64}, Vector{String}, Vector{String}, DuctAPE.var"#42#47", IntegrationOptions{GaussLegendre{Vector{Float64}, Vector{Float64}}, GaussLegendre{Vector{Float64}, Vector{Float64}}}, ChainSolverOptions{Bool, Int64, DuctAPE.ExternalSolverOptions}, SLORGridSolverOptions{Bool, Float64, Int64}}(false, true, [1], DuctAPE.var"#42#47"(), true, false, 0.05, 1.0e-15, 2.220446049250313e-15, IntegrationOptions{GaussLegendre{Vector{Float64}, Vector{Float64}}, GaussLegendre{Vector{Float64}, Vector{Float64}}}(GaussLegendre{Vector{Float64}, Vector{Float64}}([0.019855071751231856, 0.10166676129318664, 0.2372337950418355, 0.4082826787521751, 0.591717321247825, 0.7627662049581645, 0.8983332387068134, 0.9801449282487682], [0.05061426814518838, 0.11119051722668723, 0.15685332293894372, 0.18134189168918097, 0.18134189168918097, 0.15685332293894372, 0.11119051722668723, 0.05061426814518838]), GaussLegendre{Vector{Float64}, Vector{Float64}}([0.019855071751231856, 0.10166676129318664, 0.2372337950418355, 0.4082826787521751, 0.591717321247825, 0.7627662049581645, 0.8983332387068134, 0.9801449282487682], [0.05061426814518838, 0.11119051722668723, 0.15685332293894372, 0.18134189168918097, 0.18134189168918097, 0.15685332293894372, 0.11119051722668723, 0.05061426814518838])), Bool[0], ["outputs.jl"], false, ["outs"], SLORGridSolverOptions{Bool, Float64, Int64}(200, 2.220446049250313e-16, Bool[0], [0]), ChainSolverOptions{Bool, Int64, DuctAPE.ExternalSolverOptions}(DuctAPE.ExternalSolverOptions[NLsolveOptions{Bool, Float64, Int64, UnionAll, @NamedTuple{}, Symbol}(:anderson, 1.0e-12, 200, LineSearches.MoreThuente, NamedTuple(), Bool[0], [0]), MinpackOptions{Bool, Float64, Int64, Symbol}(:hybr, 1.0e-12, 100, Bool[0], [0]), NLsolveOptions{Bool, Float64, Int64, UnionAll, @NamedTuple{}, Symbol}(:newton, 1.0e-12, 20, LineSearches.MoreThuente, NamedTuple(), Bool[0], [0])], Bool[0, 0, 0], [0, 0, 0]))
For more advanced option selection, see the examples and API reference.
With the ducted_rotor input build, and the options selected, we are now ready to run an analysis. This is done simply with the analyze
function which dispatches the appropriate analysis, solve, and post-processing functions based on the selected options.
analyze(
+ ducted_rotor::DuctedRotor,
+ operating_point::OperatingPoint,
+ reference_parameters::ReferenceParameters,
+ options::Options=set_options();
+ prepost_container_caching=nothing,
+ solve_parameter_caching=nothing,
+ solve_container_caching=nothing,
+ return_inputs=false,
+)
Analyze ducted_rotor, including preprocessing.
Arguments
ducted_rotor::DuctedRotor
: DuctedRotor input object (see docstring for DuctedRotor
type)operating_point::OperatingPoint
: OperatingPoint input object (see docstring for OperatingPoint
type)reference_parameters::ReferenceParameters
: ReferenceParameters input object (see docstring for ReferenceParameters
type)options::Options=set_options()
: Options object (see set_options
and related functions)
Keyword Arguments
prepost_container_caching=nothing
: Output of allocate_prepost_container_cache
solve_parameter_caching=nothing
: Output of allocate_solve_parameter_container_cache
solve_container_caching=nothing
: Output of allocate_solve_container_cache
return_inputs=false
: flag as to whether or not to return the pre-processed inputs
Returns
outs::NamedTuple
: Named Tuple of various analysis outputs (see docstring for postprocess for details), note, if linear system decomposition fails, no solve is performed and an empty vector is returned.ins::NamedTuple
: Named Tuple of various pre-processed inputs (e.g. panels and body linear system), will only be returned if return_inputs=true
convergence_flag
: Flag for successful solve convergence
sourceouts, success_flag = DuctAPE.analyze(ducted_rotor, operating_point, reference_parameters, options)
There are many outputs contained in the named tuple output from the analyze
function (see the post_process docstring), but some that may be of immediate interest include:
# Total Thrust Coefficient
+outs.totals.CT
0.9693509063670719
# Total Torque Coefficient
+outs.totals.CQ
0.10306885180217912
In the case that one wants to run the same geometry at several different operating points, for example: for a range of advance ratios, there is another dispatch of the analyze
function that accepts an input, multipoint
, that is a vector of operating points.
analyze(
+ ducted_rotor::DuctedRotor,
+ operating_point::AbstractVector{OperatingPoint},
+ reference_parameters::ReferenceParameters,
+ options::Options=set_options();
+ prepost_container_caching=nothing,
+ solve_parameter_caching=nothing,
+ solve_container_caching=nothing,
+ return_inputs=false,
+)
Analyze ducted_rotor, including preprocessing, for a set of operating points.
Arguments
ducted_rotor::DuctedRotor
: DuctedRotor input objectoperating_point::AbstractVector{OperatingPoint}
: Vector of Operating Points at which to analyze the ducted_rotorreference_parameters::ReferenceParameters
: ReferenceParameters input object (see docstring for ReferenceParameters
type)options::Options=set_options()
: Options object
Keyword Arguments
prepost_container_caching=nothing
: Output of allocate_prepost_container_cache
solve_parameter_caching=nothing
: Output of allocate_solve_parameter_container_cache
solve_container_caching=nothing
: Output of allocate_solve_container_cache
return_inputs=false
: flag as to whether or not to return the pre-processed inputs
Returns
outs::Vector{NamedTuple}
: Vector of named tuples of various analysis outputs (see docstring for postprocess for details), note, if linear system decomposition fails, no solve is performed and an empty vector is returned.ins::NamedTuple
: Named Tuple of various pre-processed inputs (e.g. panels and body linear system), will only be returned if return_inputs=true
convergence_flag
: Flag for successful solve convergence
sourceRunning a multi-point analysis on the example geometry given there, it might look something like this:
# - Advance Ratio Range - #
+Js = range(0.0, 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 = Js * 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_vec, success_flags = DuctAPE.analyze(ducted_rotor, operating_points, reference_parameters, DuctAPE.set_options(operating_points))
There are a few things to note here.
- We want to make sure that the operating point objects we put into the input vector are unique instances.
- We need to use the dispatch of
set_options
that accepts the operating point vector to set up the right number of things in the background (like convergence flags for each operating point). - The outputs of the analysis are vectors of the same outputs for a single analysis.
For multi-point analysis outputs, which are given as a vector of output objects, we might access and plot things as follows. We also take the opportunity to present some verification against DFDC, showing that DuctAPE matches remarkably well (within 0.5%) of DFDC. We therefore first provide data from DFDC analyses of the above example geometry at various advance ratios.
# Verification Data From DFDC
+
+dfdc_jept = [
+ 0.0 0.0 0.64763 0.96692
+ 0.1 0.1366 0.64716 0.88394
+ 0.2 0.2506 0.6448 0.80785
+ 0.3 0.3457 0.64044 0.73801
+ 0.4 0.4251 0.63401 0.67382
+ 0.5 0.4915 0.62534 0.61468
+ 0.6 0.547 0.61428 0.56001
+ 0.7 0.5935 0.6006 0.50925
+ 0.8 0.6326 0.58411 0.46187
+ 0.9 0.6654 0.56452 0.41738
+ 1.0 0.693 0.54158 0.37531
+ 1.1 0.716 0.51499 0.33522
+ 1.2 0.7349 0.48446 0.2967
+ 1.3 0.7499 0.44966 0.25937
+ 1.4 0.7606 0.41031 0.2229
+ 1.5 0.7661 0.36604 0.18694
+ 1.6 0.7643 0.31654 0.15121
+ 1.7 0.7506 0.26153 0.11547
+ 1.8 0.7126 0.20061 0.07941
+ 1.9 0.61 0.13355 0.04287
+ 2.0 0.1861 0.05993 0.00558
+]
+
+dfdc_J = dfdc_jept[:,1]
+dfdc_eta = dfdc_jept[:,2]
+dfdc_cp = dfdc_jept[:,3]
+dfdc_ct = dfdc_jept[:,4]
We can then access the various multi-point analysis outputs however is convenient, we choose a broadcasting approach here:
# - extract efficiency, power, and thrust coefficients - #
+# efficiency
+eta = (p->p.totals.total_efficiency[1]).(outs_vec)
+# power
+cp = (p->p.totals.CP[1]).(outs_vec)
+# thrust
+ct = (p->p.totals.CT[1]).(outs_vec)
And then we can plot the data to compare DFDC and DuctAPE.
using Plots
+
+# set up efficiency plot
+pe = plot(; xlabel="Advance Ratio", ylabel="Efficiency")
+
+# plot DFDC data
+plot!(
+ pe,
+ dfdc_J,
+ dfdc_eta;
+ seriestype=:scatter,
+ markersize=5,
+ markercolor=plotsgray,
+ markerstrokecolor=plotsgray,
+ label="DFDC"
+)
+
+# Plot DuctAPE outputs
+plot!(pe, Js, eta; linewidth=2, color=primary, label = "DuctAPE")
+
+# setup cp/ct plot
+ppt = plot(; xlabel="Advance Ratio")
+
+# plot DFDC data
+plot!(
+ ppt,
+ dfdc_J,
+ dfdc_cp;
+ seriestype=:scatter,
+ markersize=5,
+ markercolor=plotsgray,
+ markerstrokecolor=primary,
+ markerstrokewidth=2,
+ label="DFDC Cp"
+)
+plot!(
+ ppt,
+ dfdc_J,
+ dfdc_ct;
+ seriestype=:scatter,
+ markersize=5,
+ markercolor=plotsgray,
+ markerstrokecolor=secondary,
+ markerstrokewidth=2,
+ label="DFDC Ct"
+)
+
+# plot DuctAPE outputs
+plot!(
+ ppt,
+ Js,
+ cp;
+ linewidth=1.5,
+ color=primary,
+ label="DuctAPE Cp"
+)
+plot!(
+ ppt,
+ Js,
+ ct;
+ linewidth=1.5,
+ color=secondary,
+ label="DuctAPE Ct"
+)
+
+plot(pe, ppt; size=(700,350), layout=(1,2), margin=2mm)