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

[andrewning]

The physics is quite tenuous for this method. Seems like this should be fine to change.
I assume this means that you are changing the airfoil shape during the optimization? If not, I'd definitely use precomputed polars rather than on-the-fly ones.

[dingraha]

I assume this means that you are changing the airfoil shape during the optimization? If not, I'd definitely use precomputed polars rather than on-the-fly ones.

Actually I am using precomputed polars. I create a single CCBlade.jl airfoil struct and use the same Reynolds correction, hub/tip correction, mach number correction, and rotation correction structs for the entire optimization. It appears the optimizer is able manipulate the inputs to the rotation correction in a way that reduces the drag to < 0, despite the DuSeligEggers parameters (a, b, d, m, alpha0) never changing.

[andrewning]

Ah, I mean I always precompute the rotational corrections also.

[andrewning]

This is common practice in the wind energy field. Rotational and other corrections (including extrapolation), especially when combined together can easily go off the rails. So we always create fully precomputed polars. This does mean that you have to apply the rotational correction at a notional section (typically around 70% radius as that is around the midpoint of the blade based on area). The corrections are semi-empirical and can't really be justified all the way towards the root and towards the balde.

[dingraha]

Gotcha, thanks for the explanation. I see there's an example of precomputing the rotational correction in the docs: https://flow.byu.edu/CCBlade.jl/stable/howto/#Airfoil-Data. I'll give that a shot in the future.