`DuSeligEggers` and negative drag coefficients

[dingraha]

When running optimizations with CCBlade.jl, sometimes the optimizer will figure out a way to manipulate the design such that the DuSeligEggers rotation correction will reduce the drag coefficient to < 0. As a workaround, I've been using a tweaked version of the correction:

using CCBlade: CCBlade

"""
    DuSeligEggersNonNegDeltaCD(a, b, d, m, alpha0)
    DuSeligEggersNonNegDeltaCD(a=1.0, b=1.0, d=1.0, m=2*pi, alpha0=0.0)  # uses defaults

DuSelig correction for lift an Eggers correction for drag.

**Arguments**
- `a, b, d::Float64`: parameters in Du-Selig paper.  Normally just 1.0 for each.
- `m::Float64`: lift curve slope.  Defaults to 2 pi for zero argument version.
- `alpha0::Float64`: zero-lift angle of attack.  Defaults to 0 for zero argument version.
"""
struct DuSeligEggersNonNegDeltaCD{TF} <: CCBlade.RotationCorrection
    a::TF
    b::TF
    d::TF
    m::TF
    alpha0::TF
end

DuSeligEggersNonNegDeltaCD() = DuSeligEggersNonNegDeltaCD(1.0, 1.0, 1.0, 2*pi, 0.0)


function CCBlade.rotation_correction(du::DuSeligEggersNonNegDeltaCD, cl, cd, cr, rR, tsr, alpha, phi=alpha, alpha_max_corr=30*pi/180)
    
    # Du-Selig correction for lift
    Lambda = tsr / sqrt(1 + tsr^2)
    expon = du.d / (Lambda * rR)
    fcl = 1.0/du.m*(1.6*cr/0.1267*(du.a-cr^expon)/(du.b+cr^expon)-1)

    # linear lift line
    cl_linear = du.m*(alpha - du.alpha0)

    # adjustment for max correction
    amax = atan(1/0.12) - 5*pi/180  # account for singularity in Eggers (not pi/2)
    if abs(alpha) >= amax 
        adj = 0.0
    elseif abs(alpha) > alpha_max_corr
        adj = ((amax-abs(alpha))/(amax-alpha_max_corr))^2
    else
        adj = 1.0
    end
    
    # increment in cl
    deltacl = adj*fcl*(cl_linear - cl)
    cl += deltacl

    # Eggers correction for drag
    deltacd = deltacl * (sin(phi) - 0.12*cos(phi))/(cos(phi) + 0.12*sin(phi))  # note that we can actually use phi instead of alpha as is done in airfoilprep.py b/c this is done at each iteration
    if deltacd > 0
        cd += deltacd
    end

    return cl, cd
end 

In the code above, the correction is only applied if it actually increases the drag coefficient. But I'm not familiar with the physics and wonder if that's the "right" thing to do (though it definitely prevents the negative drag coefficients I and others see occasionally).