Getting Started

Getting Started

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

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.

Paneling Constants

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.

DuctAPE.PanelingConstants — Type

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 casing inlet.
ncenterbody_inlet::Int : The number of panels to use for the centerbody inlet.
npanels::AbstractVector{Int} : A vector containing the number of panels between discrete locations inside the wake. Specifically, the number of panels between the rotors, between the last rotor and the first body trailing edge, between the body trailing edges (if different), and between the last body trailing edge and the end of the wake. The length of this vector should be N+1 (where N is the number of rotors) if the duct and centerbody trailing edges are aligned, and N+2 if not.
dte_minus_cbte::Float : An indicator concerning the hub and duct trailing edge relative locations. Should be set to -1 if the duct trailing edge axial position minus the centerbody trailing edge axial position is negative, +1 if positive (though any positive or negative number will suffice), and zero if the trailing edges are aligned.
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
)

Assembling the DuctedRotor

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,
)

Operating Point

Next we will assemble the operating point which contains information about the freestream as well as the rotor rotation rate(s).

DuctAPE.OperatingPoint — Type

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.

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)

source

# Freestream
Vinf = 30.0
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, Omega, rhoinf, muinf, asound)

Reference Parameters

The reference parameters are used in the post-processing non-dimensionalizations.

DuctAPE.ReferenceParameters — Type

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])

Set Options

The default options should be sufficient for just starting out and are set through the set_options function.

DuctAPE.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.

source

options = DuctAPE.set_options()

Options{Bool, DuctAPE.HeadsBoundaryLayerOptions{Bool, Float64, Float64, UnionAll, Int64, Float64, Float64, Float64, DuctAPE.DiffEq, Nothing, Nothing}, Vector{Bool}, Float64, Vector{Int64}, Vector{String}, Vector{String}, DuctAPE.var"#57#62", IntegrationOptions{GaussLegendre{Vector{Float64}, Vector{Float64}}, GaussLegendre{Vector{Float64}, Vector{Float64}}}, ChainSolverOptions{Bool, Int64, DuctAPE.ExternalSolverOptions}, GridSolverOptions{Bool, Float64, Int64, Symbol}}(false, true, [1], DuctAPE.var"#57#62"(), 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])), DuctAPE.HeadsBoundaryLayerOptions{Bool, Float64, Float64, UnionAll, Int64, Float64, Float64, Float64, DuctAPE.DiffEq, Nothing, Nothing}(false, true, true, false, 500, 1.0e-6, nothing, nothing, 0.001, DuctAPE.DiffEq(), RadauIIA5, 3.0, 10, 10, 0.2, 0.2, false, 0.0, 0.0001, 0.0), Bool[0], ["outputs.jl"], false, ["outs"], GridSolverOptions{Bool, Float64, Int64, Symbol}(30, 3.0e-10, :newton, :forward, false, 3, Bool[0], [0], [0.0]), ChainSolverOptions{Bool, Int64, DuctAPE.ExternalSolverOptions}(DuctAPE.ExternalSolverOptions[NLsolveOptions{Bool, Float64, Int64, UnionAll, @NamedTuple{}, Symbol}(:anderson, 1.0e-10, 200, LineSearches.MoreThuente, NamedTuple(), Bool[0], [0]), MinpackOptions{Bool, Float64, Int64, Symbol}(:hybr, 1.0e-10, 100, Bool[0], [0]), NLsolveOptions{Bool, Float64, Int64, UnionAll, @NamedTuple{}, Symbol}(:trust_region, 1.0e-10, 20, LineSearches.MoreThuente, NamedTuple(), Bool[0], [0])], Bool[0, 0, 0], [0, 0, 0]), true)

For more advanced option selection, see the examples and API reference.

Run a Single Analysis

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.

DuctAPE.analyze — Method

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

source

outs, success_flag = DuctAPE.analyze(
    ducted_rotor, operating_point, reference_parameters, options
)

Single Run Outputs

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.49951204968284824

# Total Torque Coefficient
outs.totals.CQ

0.09522170321610814

Run a Multi-Point Analysis

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.

DuctAPE.analyze — Method

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 object
operating_point::AbstractVector{OperatingPoint} : Vector of Operating Points at which to analyze the ducted_rotor
reference_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

source

Running a multi-point analysis on the example geometry given there, it might look something like this:

# - Advance Ratio Range - #
Js = 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 = 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.

Multi-point Outputs

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
c_p = (p->p.totals.CP[1]).(outs_vec)

# thrust
c_t = (p->p.totals.CT[1]).(outs_vec)

And then we can plot the data to compare DFDC and DuctAPE.

# set up efficiency plot
pe = plot(; xlabel="Advance Ratio", ylabel="Efficiency")

# plot DFDC data
plot!(
    pe,
    dfdc_J,
    dfdc_eta;
    seriestype=:scatter,
    markersize=5,
    label="DFDC",
)

# Plot DuctAPE outputs
plot!(
    pe,
    Js,
    eta;
    linewidth=2,
    label="DuctAPE",
)

# setup c_p/c_t plot
ppt = plot(; xlabel="Advance Ratio")

# plot DFDC data
plot!(
    ppt,
    dfdc_J,
    dfdc_cp;
    seriestype=:scatter,
    markersize=5,
    markerstrokewidth=2,
    label="DFDC Cp",
)

plot!(
    ppt,
    dfdc_J,
    dfdc_ct;
    seriestype=:scatter,
    markersize=5,
    markerstrokewidth=2,
    label="DFDC Ct",
)

# plot DuctAPE outputs
plot!(
    ppt,
    Js,
    c_p;
    linewidth=1.5,
    label="DuctAPE c_p",
)

plot!(
    ppt,
    Js,
    c_t;
    linewidth=1.5,
    label="DuctAPE Ct",
)

plot(
    pe,
    ppt;
    size=(700, 350),
    layout=(1, 2),
)

Getting Started

Assemble Inputs

Body Geometry

Rotor Geometry

Paneling Constants

Assembling the DuctedRotor

Operating Point

Reference Parameters

Set Options

Run a Single Analysis

Single Run Outputs

Run a Multi-Point Analysis

Multi-point Outputs