Allow user to set npts, backwardsearch in solve

[dingraha]

While using CCBlade.jl for some propeller optimization work, I ran into this interesting problem where, for a certain combination of inputs, multiple solutions to the BEMT equation exist. Usually CCBlade.jl will find the desired solution (the one that corresponds to the most-reasonable local angle of attack value) but every once in a while it gets unlucky and finds a "bad" solution, which significantly changes the usual integrated outputs (thrust, torque etc) for that particular iteration. This bothers the optimizer and slows convergence of the optimization.

Here's a script that makes two small changes to one of the CCBlade.jl test cases and shows the presence of the multiple solutions:

module MultiSolutionTestCase

using CCBlade: CCBlade
using ColorSchemes: colorschemes
using GLMakie

# Connector character for filenames.
const conn = "-"

function doit()
    colormap = colorschemes[:tab10]
    # --- rotor definition ---
    D = 1.6
    Rhub = 0.0
    Rtip = D/2
    Rhub_eff = 1e-6  # something small to eliminate hub effects
    Rtip_eff = 100.0  # something large to eliminate tip effects
    B = 2  # number of blades

    function affunc2(alpha, Re, M)

        cl = 6.2*alpha
        cd = 0.008 - 0.003*cl + 0.01*cl*cl

        clmax = 1.2
        return clamp(cl, -clmax, clmax), cd
    end 

    rotation = CCBlade.DuSeligEggers()
    rotor_no_F = CCBlade.Rotor(Rhub_eff, Rtip_eff, B; rotation)

    # --- section definitions ---

    R = D/2.0
    r = range(R/10, stop=R, length=11)
    pitch = 1.0  # pitch distance in meters.
    theta = atan.(pitch./(2*pi*r))
    chord = 0.10

    sections = CCBlade.Section.(r, chord, theta, Ref(affunc2))

    # --- inflow definitions ---
    RPM = 2100
    rho = 1.225

    phis = (-5.0:0.01:95.0) .* (pi/180)

    num_radial = length(sections)
    Vinfs = [1.0, 15.0, 30.0]

    for idxr in 1:num_radial
        fig = Figure(size=(2*600, 2*450))

        section = sections[idxr]
        for (i, Vinf) in enumerate(Vinfs)
            ax_residual = fig[i+1,1] = Axis(fig; xlabel="φ, deg", ylabel="residual, Vinf = $(Vinf) m/s")
            ax_residual_alpha = fig[i+1,2] = Axis(fig; xlabel="α = θ - φ, deg", ylabel="residual, Vinf = $(Vinf) m/s")

            Omega = RPM * pi/30 
            op = CCBlade.simple_op(Vinf, Omega, r[idxr], rho)

            # package up variables and parameters for residual
            x = [section.r, section.chord, section.theta, rotor_no_F.Rhub, rotor_no_F.Rtip, op.Vx, op.Vy, op.rho, op.pitch, op.mu, op.asound]
            p = (section.af, rotor_no_F.B, rotor_no_F.turbine, rotor_no_F.re, rotor_no_F.mach, rotor_no_F.rotation, rotor_no_F.tip)
            # pull out first argument
            residual(phi) = CCBlade.residual_and_outputs(phi, x, p)[1]

            # Now evaluate the residual for a bunch of phis.
            residuals = residual.(phis)

            l_residuals = lines!(ax_residual, phis .* (180/pi), residuals; color=colormap[i])
            l_residuals_alpha = lines!(ax_residual_alpha, (theta[idxr] .- phis) .* (180/pi), residuals; color=colormap[i])
            ylims!(ax_residual, -10, 10)
            ylims!(ax_residual_alpha, -10, 10)
        end

        Label(fig[1, :], "radial index = $(idxr) of $(num_radial)", tellwidth=false, tellheight=true)
        Vinf_label = mapreduce(v->"Vinf=$(v)$(conn)", *, Vinfs)
        fname = "residual$(conn)$(Vinf_label)idxr=$(idxr).png"
        save(fname, fig)
    end

end

end # module

The two changes to the test case are:

• The airfoil CL is clamped between -1.2 and 1.2
• The CCBlade.DuSeligEggers() rotation correction is used

The script plots the BEMT residual for a range of phi values for each radial station and for three different axial velocities. It happens that the inner-most radial station shows the multiple solution issue for the 30.0 m/s axial velocity case. Here's that plot:

The left-hand column shows the residual as a function of phi, and the right-hand column shows the residual as a function of alpha, the local angle of attack. You can see that the bottom plot shows three solutions exist.

The workaround that I found was to allow the user to set a few parameters that control the sign-change search that CCBlade.jl does in the solve function. Specifically:

• npts allows the user to set the number of phi points evaluated during the sign change search
• forcebackwardsearch, when true, forces solve to search from high phi values to low
• epsilon_everywhere, when true, tells solve to avoid evaluating the residual at pi/2, -pi/2 etc.

I set the default values for these parameters to what they are currently, FYI, so hopefully this PR wouldn't change the current behavior.

[andrewning]

looks good, thanks. yeah, I run into scenarios with multiple solutions from time to time. running into this right now actually with vortex ring state (adding a model for this)