peterdsharpe / AeroSandbox

Aircraft design optimization made fast through modern automatic differentiation. Composable analysis tools for aerodynamics, propulsion, structures, trajectory design, and much more.
https://peterdsharpe.github.io/AeroSandbox/
MIT License
690 stars 111 forks source link

Set CL instead of Alpha in OperatingPoint #21

Closed KikeM closed 4 years ago

KikeM commented 4 years ago

Describe the solution you'd like Set a lifting coefficient instead of an angle of attack to perform the calculations.

Describe alternatives you've considered For the moment I am iteratively changing the angle of attack until I converge to the desired CL.

peterdsharpe commented 4 years ago

Ah! This is a great question and definitely possible in the latest version, at least with CasLL1 (might be possible with CasVLM1 too, but not tested)!

You'll basically want to implement alpha as a free variable and set a constraint such that ap.CL == 0.7 or whatever your desired value is. Then, it'll solve the entire system of nonlinear equations together.

Here's a sample runfile! I've highlighted some key lines below it:

import copy

from aerosandbox import *
from aerosandbox.library.airfoils import generic_airfoil, generic_cambered_airfoil

opti = cas.Opti()  # Initialize an analysis/optimization environment

# Define the 3D geometry you want to analyze/optimize.
# Here, all distances are in meters and all angles are in degrees.
airplane = Airplane(
    name="Peter's Glider",
    x_ref=0,  # CG location
    y_ref=0,  # CG location
    z_ref=0,  # CG location
    wings=[
        Wing(
            name="Main Wing",
            x_le=0,  # Coordinates of the wing's leading edge
            y_le=0,  # Coordinates of the wing's leading edge
            z_le=0,  # Coordinates of the wing's leading edge
            symmetric=True,
            xsecs=[  # The wing's cross ("X") sections
                WingXSec(  # Root
                    x_le=0,  # Coordinates of the XSec's leading edge, relative to the wing's leading edge.
                    y_le=0,  # Coordinates of the XSec's leading edge, relative to the wing's leading edge.
                    z_le=0,  # Coordinates of the XSec's leading edge, relative to the wing's leading edge.
                    chord=0.18,
                    twist=2,  # degrees
                    airfoil=generic_cambered_airfoil,  # Airfoils are blended between a given XSec and the next one.
                    control_surface_type='symmetric',
                    # Flap # Control surfaces are applied between a given XSec and the next one.
                    control_surface_deflection=0,  # degrees
                ),
                WingXSec(  # Mid
                    x_le=0.01,
                    y_le=0.5,
                    z_le=0,
                    chord=0.16,
                    twist=0,
                    airfoil=generic_cambered_airfoil,
                    control_surface_type='asymmetric',  # Aileron
                    control_surface_deflection=0,
                ),
                WingXSec(  # Tip
                    x_le=0.08,
                    y_le=1,
                    z_le=0.1,
                    chord=0.08,
                    twist=-2,
                    airfoil=generic_cambered_airfoil,
                ),
            ]
        ),
        Wing(
            name="Horizontal Stabilizer",
            x_le=0.6,
            y_le=0,
            z_le=0.1,
            symmetric=True,
            xsecs=[
                WingXSec(  # root
                    x_le=0,
                    y_le=0,
                    z_le=0,
                    chord=0.1,
                    twist=-10,
                    airfoil=generic_airfoil,
                    control_surface_type='symmetric',  # Elevator
                    control_surface_deflection=0,
                ),
                WingXSec(  # tip
                    x_le=0.02,
                    y_le=0.17,
                    z_le=0,
                    chord=0.08,
                    twist=-10,
                    airfoil=generic_airfoil
                )
            ]
        ),
        Wing(
            name="Vertical Stabilizer",
            x_le=0.6,
            y_le=0,
            z_le=0.15,
            symmetric=False,
            xsecs=[
                WingXSec(
                    x_le=0,
                    y_le=0,
                    z_le=0,
                    chord=0.1,
                    twist=0,
                    airfoil=generic_airfoil,
                    control_surface_type='symmetric',  # Rudder
                    control_surface_deflection=0,
                ),
                WingXSec(
                    x_le=0.04,
                    y_le=0,
                    z_le=0.15,
                    chord=0.06,
                    twist=0,
                    airfoil=generic_airfoil
                )
            ]
        )
    ]
)
airplane.set_spanwise_paneling_everywhere(30)  # Set the resolution of your analysis

alpha = opti.variable()
opti.set_initial(alpha, 0)

ap = Casll1(  # Set up the AeroProblem
    airplane=airplane,
    op_point=OperatingPoint(
        density=1.225,  # kg/m^3
        viscosity=1.81e-5,  # kg/m-s
        velocity=10,  # m/s
        mach=0,  # Freestream mach number
        alpha=alpha,  # In degrees
        beta=0,  # In degrees
        p=0,  # About the body x-axis, in rad/sec
        q=0,  # About the body y-axis, in rad/sec
        r=0,  # About the body z-axis, in rad/sec
    ),
    opti=opti  # Pass it an optimization environment to work in
)

opti.subject_to([
    ap.CL == 0.7
])

# Solver options
p_opts = {}
s_opts = {}
# s_opts["mu_strategy"] = "adaptive"
opti.solver('ipopt', p_opts, s_opts)

# Solve
try:
    sol = opti.solve()
except RuntimeError:
    sol = opti.debug
    raise Exception("An error occurred!")

# Postprocess

# Create solved object
ap_sol = copy.deepcopy(ap)
ap_sol.substitute_solution(sol)

# ap_sol.draw(show=True)  # Generates a pretty picture!

print("CL:", ap_sol.CL)
print("CD:", ap_sol.CD)
print("CY:", ap_sol.CY)
print("Cl:", ap_sol.Cl)
print("Cm:", ap_sol.Cm)
print("Cn:", ap_sol.Cn)

Snippet of output from the file above:

       total  |   1.18 s (  1.18 s)   1.18 s (  1.18 s)         1
CL: 0.699999999999972
CD: 0.027870764828169817
CY: 0.0
Cl: 0.0
Cm: 0.018719090893452202
Cn: 0.0

Key lines from the file above:

alpha = opti.variable()
opti.set_initial(alpha, 0)

ap = Casll1(  # Set up the AeroProblem
    airplane=airplane,
    op_point=OperatingPoint(
        density=1.225,  # kg/m^3
        viscosity=1.81e-5,  # kg/m-s
        velocity=10,  # m/s
        mach=0,  # Freestream mach number
        alpha=alpha,  # In degrees
        beta=0,  # In degrees
        p=0,  # About the body x-axis, in rad/sec
        q=0,  # About the body y-axis, in rad/sec
        r=0,  # About the body z-axis, in rad/sec
    ),
    opti=opti  # Pass it an optimization environment to work in
)

opti.subject_to([
    ap.CL == 0.7
])
peterdsharpe commented 4 years ago

Closing if there are no other questions, but feel free to reopen!