Creating a C-shape wing using B-spline control points

[Cham96]

Type of issue

What types of issue is it?

• New feature (non-breaking change which adds functionality)
• Bugfix (non-breaking change which fixes an issue)
• Other (please describe)

Description

Hello Everyone, I am attempting to alter the wing mesh using Y-shear and Z-shear control points to create a C shaped wing or winglets that come back towards the root. Overall I am trying to use these B-spline control points to control the geometry of the mesh and in turn, alter the shape of the wing, however, the Y-shear doesn't seem to have any effect on the shape of the mesh. The closest I have gotten is to create winglets using the z-shear.

Steps to reproduce the issue

The following is the number of control points, the control points, and the mesh dictionary. The control points shown below are adjusted to create a C-shaped wing using y and z control points:

mesh_dict = {'num_y': 9, 'num_x': 9, 'wing_type': 'CRM', 'symmetry': True, 'chord_cos_spacing': 0, 'span_cos_spacing': 0, 'num_twist_cp': 9, 'num_yshear_cp': 9}

Generate the aerodynamic mesh based on the previous dictionary

mesh, twist_cp = generate_mesh(mesh_dict)

Create a dictionary with info and options about the aerodynamic lifting surface

surface = {

Wing definition

        'name': 'wing',          # name of the surface
        'symmetry': True,        # if true, model one half of wing
                                 # reflected across the plane y = 0
        'S_ref_type': 'wetted',  # how we compute the wing area,
                                 # can be 'wetted' or 'projected'
        'fem_model_type': 'tube',
        'mesh': mesh,
        'twist_cp': np.zeros(9),
        'yshear_cp': np.array([8.2, 8.5, 9.5, 9., 7.5, 6., 2., 1., 0.]),  # set y shear control point distribution
        'zshear_cp': np.array([3.0, 2.7, 0.5, 0.2, 0.05, 0., 0., 0., 0.0]),  # set z shear control point distribution
        'CL0': 0.0,            # CL of the surface at alpha=0
        'CD0': 0.015,            # CD of the surface at alpha=0
        # Airfoil properties for viscous drag calculation
        'k_lam': 0.05,                          # percentage of chord with laminar
                                                # flow, used for viscous drag
        't_over_c_cp': np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]),
        'c_max_t': .303,        # chordwise location of maximum (NACA0015)
                                # thickness
        'with_viscous': True,   # if true, compute viscous drag
        'with_wave': False,     # if true, compute wave drag
        }

Current behavior

The Y-shear doesn't seem to have any effect at all when the control points are changed. There's no effect even when the z-shear is turned off.

Expected behavior

The only change to the mesh is effect by the z-shear control points to give winglets at the tip. Any kind of input will be appreciated, thank you.

Code version (if relevant)

**Python version: 3.7.4

[johnjasa]

Hey @Cham96, first and foremost, thank you for your detailed issue, that really helps us understand what's going on. I'll try to spell out my suggestions here, please let me know if anything is unclear.

The solution

I have replicated your script and was able to achieve a C-wing by changing 'wing_type' : 'CRM' in the surface dictionary to 'wing_type' : 'rect'. To be clear, here's the surface dictionary code I used:

mesh_dict = {'num_y': 41,
    'num_x': 9,
    'wing_type': 'rect',
    'symmetry': True,
    'chord_cos_spacing': 0,
    'span_cos_spacing': 0,
    'num_twist_cp': 9,
    'num_yshear_cp': 9}

This produces a wing that looks like this:

Explanation

In short, the 'CRM' option prescribed a geometry that was overriding the y-shear values you set. By using the 'rect' option, which is a simpler mesh, the shear values actually affect the mesh. That being said, it loses the definition from the CRM mesh, especially the taper and sweep. If you want to add those in, you'd have to add that definition into the surface dictionary just like you did for the shear values.

Feel free to comment if you're facing other issues!

[Cham96]

Hello @johnjasa, sorry for the delayed reply, but I would like to thank you for your input. It helped me in regards to going in the right direction and saving a lot of time (very valuable right now). I am performing an Aerostructural optimization for a C wing by having Z shear as a design variable.

I would like to get your input on another issue I'm facing (could be a bug). 1) The optimization results for Lift in the C-wing from the aero structural does not have a zero lift distribution at the wingtips. 2) Utilizing wave drag seems to give me a value error (ValueError: array must not contain infs or NaNs) which is important since I'm running the wing at 0.84 Mach

I am using the following mesh dictionary:

`mesh_dict = {'num_y' : 21, 'num_x' : 9, 'wing_type': 'rect', 'span': 60.9, 'root_chord': 8.75, 'symmetry': True, 'chord_cos_spacing': 0, 'span_cos_spacing': 0, 'num_twist_cp': 5, 'num_yshear_cp': 9, 'num_zshear_cp': 9}

mesh = generate_mesh(mesh_dict)`

With the following surface dictionary for the aerodynamics and structural lifting surface

` surface = {

Wing definition

        'name' : 'wing',        # name of the surface
        'symmetry' : True,     # if true, model one half of wing
                                # reflected across the plane y = 0
        'S_ref_type' : 'wetted', # how we compute the wing area,
                                 # can be 'wetted' or 'projected'
        'fem_model_type' : 'tube',

        'thickness_cp' : np.array([.1, .2, .3]),

        'taper': 0.149,
        'sweep': 31.6,
        'twist_cp': np.zeros(5),
        'yshear_cp': np.array([9.0, 0., 0., 0., 0., 0., 0., 0., 0.]),  # set y shear control point distribution
        'zshear_cp': np.array([3.0, 2.7, 0.5, 0.2, 0.05, 0., 0., 0., 0.]),  # set z shear control point distribution
        'mesh' : mesh,

        'CL0' : 0.0,            # CL of the surface at alpha=0
        'CD0' : 0.015,            # CD of the surface at alpha=0

        # Airfoil properties for viscous drag calculation
        'k_lam' : 0.05,         # percentage of chord with laminar
                                # flow, used for viscous drag
        # 't_over_c_cp' : np.array([0.15]),      # thickness over chord ratio (NACA0015)
        't_over_c_cp': np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]),
        'c_max_t' : .303,       # chordwise location of maximum (NACA0015)
                                # thickness
        'with_viscous' : True,
        'with_wave' : False,     # if true, compute wave drag

        # Structural values are based on aluminum 7075
        'E' : 70.e9,            # [Pa] Young's modulus of the spar
        'G' : 30.e9,            # [Pa] shear modulus of the spar
        'yield' : 500.e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case
        'mrho' : 3.e3,          # [kg/m^3] material density
        'fem_origin' : 0.35,    # normalized chordwise location of the spar
        'wing_weight_ratio' : 2.,
        'struct_weight_relief' : True,    # True to add the weight of the structure to the loads on the structure
        'distributed_fuel_weight' : False,
        # Constraints
        'exact_failure_constraint' : False, # if false, use KS function
        }`

With the following constraints, objective, and design variables;

`prob.model.add_design_var('wing.twist_cp', lower=-10., upper=15.) prob.model.add_design_var('wing.geometry.t_over_c_cp', lower=0.06, upper=0.2, scaler=10.) prob.model.add_design_var('wing.geometry.zshear_cp', lower=0, upper=6.)

prob.model.add_constraint('AS_point_0.wing_perf.failure', upper=0.) prob.model.add_constraint('AS_point_0.wing_perf.thickness_intersects', upper=0.) prob.model.add_constraint('AS_point_0.L_equals_W', equals=0.)

prob.model.add_objective('AS_point_0.fuelburn', scaler=1e-5)`

I have attached the intermediate code for your reference;

`prob = Problem()

indep_var_comp = IndepVarComp() indep_var_comp.add_output('v', val=248.136, units='m/s') indep_var_comp.add_output('alpha', val=5., units='deg') indep_var_comp.add_output('Mach_number', val=0.84) indep_var_comp.add_output('re', val=1.e6, units='1/m') indep_var_comp.add_output('rho', val=0.38, units='kg/m*3') indep_var_comp.add_output('CT', val=grav_constant 17.e-6, units='1/s') indep_var_comp.add_output('R', val=11.165e6, units='m') indep_var_comp.add_output('W0', val=0.481 * 3e5, units='kg') indep_var_comp.add_output('speed_of_sound', val=295.4, units='m/s') indep_var_comp.add_output('load_factor', val=1.) indep_var_comp.add_output('empty_cg', val=np.zeros(3), units='m')

prob.model.add_subsystem('prob_vars', indep_var_comp, promotes=['*'])

aerostruct_group = AerostructGeometry(surface=surface) name = 'wing'

prob.model.add_subsystem(name, aerostruct_group)

point_name = 'AS_point_0'

AS_point = AerostructPoint(surfaces=[surface])

prob.model.add_subsystem(point_name, AS_point, promotes_inputs=['v', 'alpha', 'Mach_number', 're', 'rho', 'CT', 'R', 'W0', 'speed_of_sound', 'empty_cg', 'load_factor'])

com_name = point_name + '.' + name + '_perf' prob.model.connect(name + '.local_stiff_transformed', point_name + '.coupled.' + name + '.local_stiff_transformed') prob.model.connect(name + '.nodes', point_name + '.coupled.' + name + '.nodes')

prob.model.connect(name + '.mesh', point_name + '.coupled.' + name + '.mesh')

prob.model.connect(name + '.radius', com_name + '.radius') prob.model.connect(name + '.nodes', com_name + '.nodes') prob.model.connect(name + '.cg_location', point_name + '.' + 'total_perf.' + name + '_cg_location') prob.model.connect(name + '.structural_mass', point_name + '.' + 'total_perf.' + name + '_structural_mass') prob.model.connect(name + '.t_over_c', com_name + '.t_over_c')

from openmdao.api import ScipyOptimizeDriver prob.driver = ScipyOptimizeDriver() prob.driver.options['tol'] = 1e-9 prob.driver.options['optimizer'] = 'SLSQP' prob.driver.options['debug_print'] = ['nl_cons', 'objs', 'desvars']

recorder = SqliteRecorder("aerostruct.db") prob.driver.add_recorder(recorder) prob.driver.recording_options['record_derivatives'] = True prob.driver.recording_options['includes'] = ['*']`

The results I'm obtaining for the aero structural optimization for a C-wing with a Zshear design variable;

I ran the same for the aerodynamic optimization and I obtained the following;

The Lift distribution doesn't go to zero at the wingtips (the elliptical lift does for AeroStruct but not for Aero) and I can't incorporate wave drag without getting the value error.

Any kind of input will be greatly appreciated, and thank you again for all your help @johnjasa .

[johnjasa]

Hey @Cham96, thanks again for a great write-up. Could you please share the full scripts you're using for each case so I can take a look? I appreciate the code snippets here, but I want to make sure I have everything right.

Regarding 1.), the issue with C-wings and other non-planar wings is that the visualization tool was not designed with them in mind, and some of the readouts are not correct. This is just concerning the visualization, not the actual analysis or optimization, though. For example, this line in plot_wing.py assumes that the last entry in the mesh is the furthest outboard point, but that isn't the case for a C-wing. To process non-planar wings correctly and generally would be a non-trivial amount of work that we're currently not planning on pursuing.

I would recommend postprocessing your results using a custom script or plotting utility to make sure that it's handling the non-planarity as you expect.

For 2.), I'm guessing that there's some value during the analysis or optimization upstream of that calculation that becomes NaN at some point, and that is the only place where the error is picked up. If you're able to share the scripts, I can run it on my end and see if I get the same error. Is it during analysis also or only optimization?

[Cham96]

Hello @johnjasa, Thank you for your reply. I did look more into plot_wing.py and realized this might have been the case, but your explanation reassures this and makes a lot more sense. I did manage to scale the horizontal spanwise direction and get the following results, which are satisfactory for what I'm trying to analyze.

For 2.) I get this error during optimization for incorporating the wave drag into the aero-structural analysis. The obtained the results above by having only the viscous drag turned on. The full code is attached below. Thank you again for your support. I am also using OpenMDAO 2.9.1.

Note: There seems to have been a code change five days ago to geometry_group.py at openaerostruct/geometry/geometry_group.py. Before that, I did manually change the b-spline components from a 'uniform' distribution to a 'sine' distribution, which I'm hoping doesn't have an impact on the wave drag computations. (I can always change it back if that's causing an issue.)

` import numpy as np

from openaerostruct.geometry.utils import generate_mesh from openaerostruct.integration.aerostruct_groups import AerostructGeometry, AerostructPoint from openaerostruct.utils.constants import grav_constant

from openmdao.api import IndepVarComp, Problem, Group, SqliteRecorder

Create a dictionary to store options about the surface

mesh_dict = {'num_y' : 21, 'num_x' : 7, 'wing_type': 'rect', 'span': 60.9, 'root_chord': 8.75, 'symmetry': True, 'chord_cos_spacing': 0, 'span_cos_spacing': 0, 'num_twist_cp': 5, 'num_yshear_cp': 9, 'num_zshear_cp': 9}

mesh = generate_mesh(mesh_dict)

surface = {

Wing definition

        'name' : 'wing',        # name of the surface
        'symmetry' : True,     # if true, model one half of wing
                                # reflected across the plane y = 0
        'S_ref_type' : 'wetted', # how we compute the wing area,
                                 # can be 'wetted' or 'projected'
        'fem_model_type' : 'tube',
        'thickness_cp' : np.array([.1, .2, .3]),
        'taper': 0.149,
        'sweep': 31.6,
        'twist_cp': np.zeros(5),
        'yshear_cp': np.array([9.0, 0., 0., 0., 0., 0., 0., 0., 0.]),  # set y shear control point distribution
        'zshear_cp': np.array([3.0, 2.7, 0.5, 0.2, 0.05, 0., 0., 0., 0.]),  # set z shear control point distribution
        'mesh' : mesh,

        'CL0' : 0.0,            # CL of the surface at alpha=0
        'CD0' : 0.015,            # CD of the surface at alpha=0

        # Airfoil properties for viscous drag calculation
        'k_lam' : 0.05,         # percentage of chord with laminar
                                # flow, used for viscous drag
        't_over_c_cp': np.array([0.08, 0.08, 0.08, 0.10, 0.10, 0.08]),
        'c_max_t' : .303,       # chordwise location of maximum (NACA0015)
                                # thickness
        'with_viscous' : True,
        'with_wave' : True,     # if true, compute wave drag

        # Structural values are based on aluminum 7075
        'E' : 70.e9,            # [Pa] Young's modulus of the spar
        'G' : 30.e9,            # [Pa] shear modulus of the spar
        'yield' : 500.e6 / 2.5, # [Pa] yield stress divided by 2.5 for limiting case
        'mrho' : 3.e3,          # [kg/m^3] material density
        'fem_origin' : 0.35,    # normalized chordwise location of the spar
        'wing_weight_ratio' : 2.,
        'struct_weight_relief' : True,    # True to add the weight of the structure to the loads on the structure
        'distributed_fuel_weight' : False,
        # Constraints
        'exact_failure_constraint' : False, # if false, use KS function
        }

prob = Problem()

Add problem information as an independent variables component

indep_var_comp = IndepVarComp() indep_var_comp.add_output('v', val=248.136, units='m/s') indep_var_comp.add_output('alpha', val=5., units='deg') indep_var_comp.add_output('Mach_number', val=0.84) indep_var_comp.add_output('re', val=1.e6, units='1/m') indep_var_comp.add_output('rho', val=0.38, units='kg/m*3') indep_var_comp.add_output('CT', val=grav_constant 17.e-6, units='1/s') indep_var_comp.add_output('R', val=11.165e6, units='m') indep_var_comp.add_output('W0', val=0.481 * 3e5, units='kg') indep_var_comp.add_output('speed_of_sound', val=295.4, units='m/s') indep_var_comp.add_output('load_factor', val=1.) indep_var_comp.add_output('empty_cg', val=np.zeros(3), units='m')

prob.model.add_subsystem('prob_vars', indep_var_comp, promotes=['*'])

aerostruct_group = AerostructGeometry(surface=surface)

name = 'wing'

Add tmp_group to the problem with the name of the surface.

prob.model.add_subsystem(name, aerostruct_group)

point_name = 'AS_point_0'

AS_point = AerostructPoint(surfaces=[surface])

prob.model.add_subsystem(point_name, AS_point, promotes_inputs=['v', 'alpha', 'Mach_number', 're', 'rho', 'CT', 'R', 'W0', 'speed_of_sound', 'empty_cg', 'load_factor'])

com_name = point_name + '.' + name + '_perf' prob.model.connect(name + '.local_stiff_transformed', point_name + '.coupled.' + name + '.local_stiff_transformed') prob.model.connect(name + '.nodes', point_name + '.coupled.' + name + '.nodes')

Connect aerodynamic mesh to coupled group mesh

prob.model.connect(name + '.mesh', point_name + '.coupled.' + name + '.mesh')

Connect performance calculation variables

prob.model.connect(name + '.radius', com_name + '.radius') prob.model.connect(name + '.nodes', com_name + '.nodes') prob.model.connect(name + '.cg_location', point_name + '.' + 'total_perf.' + name + '_cg_location') prob.model.connect(name + '.structural_mass', point_name + '.' + 'total_perf.' + name + '_structural_mass') prob.model.connect(name + '.t_over_c', com_name + '.t_over_c')

from openmdao.api import ScipyOptimizeDriver prob.driver = ScipyOptimizeDriver() prob.driver.options['tol'] = 1e-9 prob.driver.options['optimizer'] = 'SLSQP' prob.driver.options['debug_print'] = ['nl_cons', 'objs', 'desvars']

recorder = SqliteRecorder("aerostruct.db") prob.driver.add_recorder(recorder) prob.driver.recording_options['record_derivatives'] = True prob.driver.recording_options['includes'] = ['*']

Setup problem and add design variables, constraint, and objective

prob.model.add_design_var('wing.twist_cp', lower=-10., upper=15.) prob.model.add_design_var('wing.geometry.t_over_c_cp', lower=0.06, upper=0.2, scaler=10.) prob.model.add_design_var('wing.geometry.zshear_cp', lower=0, upper=6.) prob.model.add_constraint('AS_point_0.wing_perf.failure', upper=0.) prob.model.add_constraint('AS_point_0.wing_perf.thickness_intersects', upper=0.)

prob.model.add_constraint('AS_point_0.L_equals_W', equals=0.) prob.model.add_objective('AS_point_0.fuelburn', scaler=1e-5)

Set up the problem

prob.setup(check=True) `

[johnjasa]

Great job scaling the visualization! Nice to see. The changes in the past two weeks are only to update OpenAeroStruct's API to be compatible with OpenMDAO v3.0+. If you keep using OpenMDAO v2.9.1, there's no need to upgrade, as there's no increased capability.

I've got your aerostructural optimization running here with wave_drag, and I too get an error. What's happening is that the aerostructural system cannot converge to a valid solution. This means, physically, that the aerodynamic loads that the structural support is experiencing are too great for that structure. You can tell the coupled system is diverging instead of converging with this output:

Design Vars
{'wing.geometry.indep_vars.t_over_c_cp': array([0.06, 0.06, 0.2 , 0.2 , 0.2 , 0.2 ]),
 'wing.geometry.indep_vars.twist_cp': array([-0.0440333 ,  2.99978183,  6.86420281, 10.60009025,  7.69829164]),
 'wing.geometry.indep_vars.zshear_cp': array([2.83757544, 2.36408935, 0.31119383, 0.11182017, 0.05150805,
       0.14093848, 0.19042922, 0.16124922, 0.28119624])}

==================
AS_point_0.coupled
==================
NL: NLBGS 0 ; 643874.584 1
NL: NLBGS 1 ; 8200420.26 12.7360521
NL: NLBGS 2 ; 113333602 176.018132
NL: NLBGS 3 ; 388030133 602.648626
NL: NLBGS 4 ; 1.19852144e+09 1861.42065
NL: NLBGS 5 ; 2.32627199e+09 3612.92719
NL: NLBGS 6 ; 1.6444751e+10 25540.3014
NL: NLBGS 7 ; 7.65527363e+10 118893.863
NL: NLBGS 8 ; 1.23987019e+12 1925639.27
NL: NLBGS 9 ; 5.36861327e+13 83379797.9
NL: NLBGS 10 ; 1.05703325e+19 1.64167568e+13
NL: NLBGS 11 ; 2.28295415e+57 3.54565035e+51

Those numbers that are being printed are the residuals of the coupled system. When they go to 0, we have a valid aerostructural solution. However, in this case, they "explode," showing that the solver could not find a valid solution. The inclusion of wave drag itself isn't causing it, but that adds something to the analysis that prevents this particular design from converging.

So, there are a lot of things you could do to attempt to fix this. Given the complex wing and that I don't know your goals, I'm not sure which option is best. You could:

• make sure the flight conditions and constraints you're flying at are exactly what you want
• change the thickness of your structural supports
• modify the initial shape of this wing
• use a different number of num_x and num_y for the analysis discretization
• restrict the design space of the optimization so the optimizer doesn't choose an invalid design

I hope this helps! There are definitely a lot of nuances to optimizing a wing this complex.

[Cham96]

Hello, @johnjasa , thank you for all your support. The list of possible fixtures was definitely insightful, and I have already looked into most of them. I believe I have resolved to find a configuration that works the best. I completely agree with you about the complexity and all the different components that take to make the optimization work. I'm glad to have used OpenAeroStruct; it is a great tool to analyze such a complex system! Thank you again.

[johnjasa]

I'm glad to hear you found a good solution for your work! Thanks for the nice feedback, and don't hesitate to reach out with any more questions or ideas.