Infeasible Problem Detected - CasLL1 alpha sweep with custom Cl, Cd and Cm functions

[alexandroslessis]

Hey Peter, I'm trying to use casll1 as a solver instead of optimizing and I'm performing an alpha sweep analysis on a conventional configuration using my own functions for the sectional Cl, Cd and Cm. These functions are just polynomial fits to data obtained from CFD analyses on the airfoil for the Reynolds and Mach number of interest. I'm defining the aircraft and running the alpha sweep the same way presented in the examples with the addition of the custom airfoil functions. However, depending on the operating point, the solver raises an error, which can be read using the cas.Opti().debug method:

sol
Opti(Opti {
  instance #0
  #variables: 2 (nx = 144)
  #parameters: 1 (np = 1)
  #constraints: 2 (ng = 144)
  CasADi solver allocated.
  CasADi solver was called: Infeasible_Problem_Detected
})

I played around with the solver settings and various paneling combinations, however, the issue persisted. Is this an issue with the casll1 solver or does it have to do with CasADI? Is there a way to use the casll1 solver without setting up an optimization environment?

I am also attaching the code I'm using along with a json file containing all the data required for the aircraft and operating point definition. aircraft_casll1.zip

Sorry for the long question. I really think what you are building here is great!

[peterdsharpe]

Thanks for writing in, and especially thanks for providing the steps to reproduce this - taking a look at this now!

One quick observation - you can comment out line 441 of aero_solver.py:

    # prob.setup(run_symmetric_if_possible=True)

This setup happens as part of the Casll1 constructor (you have this in line 425), so cutting this redundancy will save you a hair of computational time :)

Looking over the rest now!

[peterdsharpe]

Cool - I think I've found the issue!

Why it's failing: Local alphas are too high

If I'm interpreting your setup correctly, this is what the 2D lift curve for your wing airfoil looks like:

While one of CasLL1's biggest strengths is that it can solve with nonlinear 2D data, this lift curve slope is exceptionally nonlinear (significant decambering by 6 degrees; CL_max at 8.5 degrees). It's possible that the solver is getting stuck due to this strong nonlinearity. At the alpha = 10 case, this is what the local angles of attack look like at each section:

Noting that some sections are reporting local angles of attack as high as 13.8 degrees, it's possible that these are nonlinear enough to cause solution failure above this point.

Why local alphas are too high

So, this is puzzling, because you don't typically see local angles of attack that ever significantly exceed the local angle of attack - if anything, they're usually significantly lower than the freestream angle of attack due to downwash.

So what's happening? Why are the local alphas peaking so high (13.8 degrees, -14.8 degrees)? There's something in the solution that's causing a massive unhandled circulation jump - let's see if we can find it. Here's what the local angles of attack look like:

At first, this looks normal - but if we zoom in on the wing, we find this:

It looks like flaps are being implemented here, but their circulation jump is just too sharp to be resolved correctly by our discretization.

You might try spreading out the flap transition length (in the spanwise direction) or reducing the nonlinearity of the lift curves to see if that makes a difference. You could also try upping the paneling resolution (a quick and dirty way is by adding something like airplane.set_spanwise_paneling_everywhere(12) after constructing the airplane).

One hot fix

One thing that you can instantly due to (partially) alleviate some of these issue is to adjust how the analysis initializes itself. By supplying a better initial guess, you'll deal with nonlinearities better. It turns out that a great initial guess is just the last solution, so you can just pop these two lines into your code as an instant upgrade: (near line 475 of aero_solver.py, new lines denoted with # NEW LINE)

        # Solve
        sol = opti.solve()
        opti.set_initial(sol.value_variables()) # NEW LINE: Initialize primals from last solution
        opti.set_initial(opti.lam_g, sol.value(opti.lam_g)) # NEW LINE: Initialize duals from last solution
        # Create solved object
        prob_sol = copy.deepcopy(prob)
        prob_sol.substitute_solution(sol)

By adding these two lines, I was able to keep this converging to alpha = 12 or so! That gives this happy lil airplane (drawn with prob_sol.draw()):

Beyond 12 degrees, my hunch is that this is just too nonlinear with the existing lift curves - but the steps listed above may help!

[alexandroslessis]

Thanks a lot for your prompt and very detailed response! I suppose I have to accept the limitations of the lifting line theory when it comes to very non-linear cases. Below I am attaching the lift curve of my double-slotted flapped airfoil. I guess these excessively non-linear cases cannot be handled by lifting line theory. It would not even converge at 0 degrees.

On the other hand, removing the intermediate spanwise wing sections and applying a finer paneling (25 spanwise panels) allowed me to extend the solution up to 15 degrees, even using the lift curve you showed above. My issue might have something to do with the transition between the flapped and unflapped regions. Increasing the panel density near these regions didn't really make a difference, even if I spread out the flap transition length. Since I was able to get converged solutions with the unflapped wing, I will just implement the flap effects using empirical methods provided by DATCOM. I will also try and simplify the polar fits to something like an exponential and see if that makes a difference.

One more question regarding the definition of the airfoil functions. Is the lambda method the only way to define those functions? And if so, I was trying to use if statements to define different fits depending on the angle of attack, however "alpha" is not compatible with boolean operations, as it's type is "MX", created using the CasADI package. Is there a way to bypass that?

Traceback (most recent call last):
  File ".\aero_solver.py", line 535, in <module>
    data = aero_solver('aero_generate.json', 'cruise')
  File ".\aero_solver.py", line 459, in aero_solver
    opti=opti
  File "C:\Python37\lib\site-packages\aerosandbox\aerodynamics\casll1.py", line 23, in __init__      
    self.setup()
  File "C:\Python37\lib\site-packages\aerosandbox\aerodynamics\casll1.py", line 66, in setup
    self.calculate_vortex_strengths()
  File "C:\Python37\lib\site-packages\aerosandbox\aerodynamics\casll1.py", line 414, in calculate_vortex_strengths
    ) for i in range(self.n_panels)
  File "C:\Python37\lib\site-packages\aerosandbox\aerodynamics\casll1.py", line 414, in <listcomp>   
    ) for i in range(self.n_panels)
  File "C:\Python37\lib\site-packages\aerosandbox\aerodynamics\casll1.py", line 151, in <lambda>     
    deflection=inner_xsec.control_surface_deflection
  File ".\aero_solver.py", line 197, in <lambda>
    airfoil_pC(pClw, alpha)
  File ".\aero_solver.py", line 9, in airfoil_pC
    if aoa <= 10:
  File "C:\Python37\lib\site-packages\casadi\casadi.py", line 9651, in __bool__
    return _casadi.MX___bool__(self, *args)
RuntimeError: .../casadi/core/mx_node.cpp:181: Can only determine truth value of a numeric MX.

[peterdsharpe]

Hi @alexandroslessis !

Will respond in more detail later, but to do a conditional statement, use:

asb.cas.if_else(conditional, value_if_true, value_if_false)

like:

asb.cas.if_else(alpha > 5, 0.1*x, 0.5)

Where asb is your AeroSandbox import name.