NACA airfoil max. camber as opti.variable

[jannav984]

Dear @peterdsharpe and colleagues,

I'm trying to optimize wing camber distribution along it span. I'm using asb.VortexLatticeMethod to evaluate wing CL and CDi. But there is a problem when I set max_camber of NACA as opti.variable. I get the error message shown below. So I'm wondering if this is implemented in AeroSandBox? (It worked fine when I used wing twist as opti.variable).

raise Exception("Implicit conversion of symbolic CasADi type to numeric matrix not supported.\n"
Exception: Implicit conversion of symbolic CasADi type to numeric matrix not supported.
This may occur when you pass a CasADi object to a numpy function.
Use an equivalent CasADi function instead of that numpy function.

My code is as follows:

import aerosandbox as asb
import aerosandbox.numpy as np


def get_aero(Nsa, camber, alpha):
    # depsg = np.concatenate((0, depsg), axis=0)
    xtip = 0  # m
    ytip = 1.9  # m
    ztip = 0  # m
    croot = 0.6
    ctip = 0.6
    epsg = np.zeros(Nsa)
    xmac = 0
    cmac = 0.6    
    xsa = np.linspace(0, xtip, Nsa)
    ysa = np.linspace(0, ytip, Nsa)
    zsa = np.linspace(0, ztip, Nsa)
    csa = np.linspace(croot, ctip, Nsa)
    asb.geometry.airfoil.airfoil_families._default_n_points_per_side = 200

    coord = asb.geometry.airfoil.airfoil_families.get_NACA_coordinates(max_camber=camber, camber_loc=0.5, thickness=0.1)

    aero_model = asb.Airplane(
        name="Bwing",
        xyz_ref=[xmac + 0.25 * cmac, 0, 0],
        wings=[
            asb.Wing(
                name="Wing",
                xyz_le=[0, 0, 0],
                symmetric=True,
                xsecs=[
                    asb.WingXSec(
                        xyz_le=[xsa[i], ysa[i], zsa[i]],
                        chord=csa[i],
                        twist=epsg[i],                        
                        airfoil=asb.Airfoil(name="NACA", coordinates=coord)
                    ) for i in range(Nsa)
                ]
            )
        ]
    )
    vlm = asb.VortexLatticeMethod(
        airplane=aero_model,
        op_point=asb.OperatingPoint(
            velocity=21.3,  # m/s
            alpha=alpha  # degree
        ),
        spanwise_resolution=3,
        spanwise_spacing_function=np.linspace,
        verbose=True
    )
    return vlm.run()


def main():
    print(asb.__version__)
    Nsa = 9
    CLw = 0.609
    opti = asb.Opti()
    alpha = opti.variable(init_guess=5.3, lower_bound=2., upper_bound=9.)
    camber_init = 0.02  # * np.ones((Nsa,))
    camber_lo = 0.02 #np.ones((Nsa,)) * 0.02
    camber_up = 0.08 #np.ones((Nsa,)) * 0.08
    camber_sec = opti.variable(init_guess=camber_init, lower_bound=camber_lo, upper_bound=camber_up)

    aero = get_aero(Nsa, camber_sec, alpha)
    opti.subject_to(aero["CL"] == CLw)
    opti.minimize((-1 * aero["CL"] ** 3 / 2) / aero["CD"])
    # opti.set_value()
    sol = opti.solve(max_iter=15, behavior_on_failure="return_last")
    print(sol.value(camber_sec))
    print(sol.value(alpha))

    return


if __name__ == "__main__":
    print("Running test")
    main()

Thanks.

Kind regards,
Jan

[peterdsharpe]

Hey Jan,

Thanks for writing in! A good question.

Cause

The high-level cause of this error is that:

1. asb.VortexLatticeMethod, asb.LiftingLine, and similar aero analyses need to discretize the wing geometry to work (i.e., make a surface mesh).
2. If the wing's airfoils are represented as coordinates (i.e., as an asb.Airfoil object), this meshing logic currently needs to be able to identify which index of the airfoil coordinates array is the leading edge.
3. If the asb.Airfoil.coordinates variable is an optimization variable (as opposed to a non-changing NumPy array), the solver can't "guarantee" that the leading edge index won't change throughout airfoil shape optimization as coordinates move around.
4. If the leading edge index can't be guaranteed to be the same, then the type of automatic differentiation that AeroSandbox uses (static-graph AD) mostly doesn't work.

Solution for Immediate Problem

First, for the immediate question - the fix here is to use asb.KulfanAirfoil instead of asb.Airfoil, which represents the Airfoil in a more shape-optimization-friendly way. More details on asb.KulfanAirfoil here.

Implemented, the fixed code looks like:

import aerosandbox as asb
import aerosandbox.numpy as np

af_thickness_mode = asb.Airfoil(
    name="NACA Thickness",
    coordinates=asb.geometry.airfoil.airfoil_families.get_NACA_coordinates(max_camber=0, camber_loc=0.5, thickness=0.1)
).to_kulfan_airfoil()
af_camber_mode = asb.Airfoil(
    name="NACA Camber",
    coordinates=asb.geometry.airfoil.airfoil_families.get_NACA_coordinates(max_camber=0.01, camber_loc=0.5, thickness=0)
).to_kulfan_airfoil()


def get_aero(Nsa, camber, alpha):
    # depsg = np.concatenate((0, depsg), axis=0)
    xtip = 0  # m
    ytip = 1.9  # m
    ztip = 0  # m
    croot = 0.6
    ctip = 0.6
    epsg = np.zeros(Nsa)
    xmac = 0
    cmac = 0.6
    xsa = np.linspace(0, xtip, Nsa)
    ysa = np.linspace(0, ytip, Nsa)
    zsa = np.linspace(0, ztip, Nsa)
    csa = np.linspace(croot, ctip, Nsa)

    af = asb.KulfanAirfoil(
        **{
            k: camber / 0.01 * af_camber_mode.kulfan_parameters[k] + af_thickness_mode.kulfan_parameters[k]
            for k in af_camber_mode.kulfan_parameters.keys()
        }
    )

    aero_model = asb.Airplane(
        name="Bwing",
        xyz_ref=[xmac + 0.25 * cmac, 0, 0],
        wings=[
            asb.Wing(
                name="Wing",
                xyz_le=[0, 0, 0],
                symmetric=True,
                xsecs=[
                    asb.WingXSec(
                        xyz_le=[xsa[i], ysa[i], zsa[i]],
                        chord=csa[i],
                        twist=epsg[i],
                        airfoil=af
                    ) for i in range(Nsa)
                ]
            )
        ]
    )
    vlm = asb.VortexLatticeMethod(  # NOTE: Recommend substituting this for asb.LiftingLine for this particular problem
        airplane=aero_model,
        op_point=asb.OperatingPoint(
            velocity=21.3,  # m/s
            alpha=alpha  # degree
        ),
        spanwise_resolution=3,
        spanwise_spacing_function=np.linspace,
        verbose=True
    )
    return vlm.run()


def main():
    print(asb.__version__)
    Nsa = 9
    CLw = 0.609
    opti = asb.Opti()
    alpha = opti.variable(init_guess=5.3, lower_bound=2., upper_bound=9.)
    camber_init = 0.02  # * np.ones((Nsa,))
    camber_lo = 0.02  # np.ones((Nsa,)) * 0.02
    camber_up = 0.08  # np.ones((Nsa,)) * 0.08
    camber_sec = opti.variable(init_guess=camber_init, lower_bound=camber_lo, upper_bound=camber_up)

    aero = get_aero(Nsa, camber_sec, alpha)
    opti.subject_to(aero["CL"] == CLw)
    opti.minimize((-1 * aero["CL"] ** 3 / 2) / aero["CD"])
    # opti.set_value()
    sol = opti.solve(max_iter=15, behavior_on_failure="return_last")
    print(sol.value(camber_sec))
    print(sol.value(alpha))

    return


if __name__ == "__main__":
    print("Running test")
    main()

Basically what we're doing here is pre-computing Kulfan-parameterized airfoil mode shapes representing "thickness" and "camber", then linearly combining them to achieve the desired NACA airfoil shape.

Finally, on a physics level, you may want to substitute asb.VortexLatticeMethod for asb.LiftingLine here, as the latter includes profile drag in CD while the former does not. LiftingLine is a drop-in replacement for VortexLatticeMethod here, so it's no extra work to implement. This will create a more accurate shape optimization problem - VLM is inviscid only, while LiftingLine includes both inviscid and viscous effects. The camber vs. alpha tradeoff at a section level is mostly determined by viscous effects, so LiftingLine is better-suited here.

Notes for ASB Devs

From this issue, there are a few actionable outcomes the team can work on:

• We should make more helpful error messages for Airfoil equality checking - the direct error is coming from this equality check, which errors on MX-to-MX array comparison. Breaking commit introduced here: f8fa57c

• We should investigate whether the VLM meshing algorithm can be reworked so that it doesn't need to compute the LE index directly. Perhaps this can be done by internally converting to KulfanAirfoil and then using those parameters to regenerate the mean camber line?

[jannav984]

Dear Peter,

thank you for your comprehesive reply and the working code!

You are doing greatjob with ABS and Neuralfoil!

Jan