diff --git a/openaerostruct/__init__.py b/openaerostruct/__init__.py index fab833f38..b482efeb4 100644 --- a/openaerostruct/__init__.py +++ b/openaerostruct/__init__.py @@ -1 +1 @@ -__version__ = "2.6.1" +__version__ = "2.6.2" diff --git a/openaerostruct/docs/user_reference/mesh_surface_dict.rst b/openaerostruct/docs/user_reference/mesh_surface_dict.rst index e866cefc4..d1af3aa25 100644 --- a/openaerostruct/docs/user_reference/mesh_surface_dict.rst +++ b/openaerostruct/docs/user_reference/mesh_surface_dict.rst @@ -91,6 +91,14 @@ The surface dict will be provided to Groups, including ``Geometry``, ``AeroPoint - np.array([0, 5]) - deg - B-spline control points for twist distribution. Array convention is ``[wing tip, ..., root]`` in symmetry cases, and ``[tip, ..., root, ... tip]`` when ``symmetry = False``. + * - chord_cp + - np.array([0.1, 5]) + - m + - B-spline control points for chord distribution. Array convention is the same than ``twist_cp``. + * - chord_scaling_pos + - 0.25 + - + - Chord position at which the chord scaling factor is applied. 1 is the trailing edge, 0 is the leading edge. .. list-table:: Aerodynamics definitions :widths: 20 20 5 55 diff --git a/openaerostruct/examples/rectangular_wing/opt_chord.py b/openaerostruct/examples/rectangular_wing/opt_chord.py index 030ed202e..59038a152 100644 --- a/openaerostruct/examples/rectangular_wing/opt_chord.py +++ b/openaerostruct/examples/rectangular_wing/opt_chord.py @@ -57,6 +57,7 @@ "S_ref_type": "projected", # how we compute the wing area, # can be 'wetted' or 'projected' "chord_cp": np.ones(3), # Define chord using 3 B-spline cp's + "chord_scaling_pos": 0.25, # Define the chord scaling position. 0 is the leading edge, 1 is the trailing edge. # distributed along span "mesh": mesh, # Aerodynamic performance of the lifting surface at diff --git a/openaerostruct/geometry/geometry_mesh.py b/openaerostruct/geometry/geometry_mesh.py index 9e9a01cc9..03dcfd821 100644 --- a/openaerostruct/geometry/geometry_mesh.py +++ b/openaerostruct/geometry/geometry_mesh.py @@ -15,6 +15,7 @@ ShearZ, Rotate, ) +import warnings class GeometryMesh(om.Group): @@ -77,12 +78,23 @@ def setup(self): # 2. Scale X val = np.ones(ny) + chord_scaling_pos = 0.25 # if no scaling position is specified : chord scaling w.r.t quarter of chord if "chord_cp" in surface: promotes = ["chord"] + if "chord_scaling_pos" in surface: + chord_scaling_pos = surface["chord_scaling_pos"] else: + if "chord_scaling_pos" in surface: + warnings.warn( + "Chord_scaling_pos has been specified but no chord design variable available", stacklevel=2 + ) promotes = [] - self.add_subsystem("scale_x", ScaleX(val=val, mesh_shape=mesh_shape), promotes_inputs=promotes) + self.add_subsystem( + "scale_x", + ScaleX(val=val, mesh_shape=mesh_shape, chord_scaling_pos=chord_scaling_pos), + promotes_inputs=promotes, + ) # 3. Sweep diff --git a/openaerostruct/geometry/geometry_mesh_transformations.py b/openaerostruct/geometry/geometry_mesh_transformations.py index 89ac16cec..ea5f8ba5e 100644 --- a/openaerostruct/geometry/geometry_mesh_transformations.py +++ b/openaerostruct/geometry/geometry_mesh_transformations.py @@ -131,11 +131,16 @@ def initialize(self): """ self.options.declare("val", desc="Initial value for chord lengths") self.options.declare("mesh_shape", desc="Tuple containing mesh shape (nx, ny).") + self.options.declare( + "chord_scaling_pos", + default=0.25, + desc="float which indicates the chord fraction where to do the chord scaling. 1. is trailing edge, 0. is leading edge", + ) def setup(self): mesh_shape = self.options["mesh_shape"] val = self.options["val"] - + self.chord_scaling_pos = self.options["chord_scaling_pos"] self.add_input("chord", units="m", val=val) self.add_input("in_mesh", shape=mesh_shape, units="m") @@ -166,9 +171,9 @@ def compute(self, inputs, outputs): te = mesh[-1] le = mesh[0] - quarter_chord = 0.25 * te + 0.75 * le + chord_pos = self.chord_scaling_pos * te + (1 - self.chord_scaling_pos) * le - outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - quarter_chord, chord_dist) + quarter_chord + outputs["mesh"] = np.einsum("ijk,j->ijk", mesh - chord_pos, chord_dist) + chord_pos def compute_partials(self, inputs, partials): mesh = inputs["in_mesh"] @@ -176,9 +181,9 @@ def compute_partials(self, inputs, partials): te = mesh[-1] le = mesh[0] - quarter_chord = 0.25 * te + 0.75 * le + chord_pos = self.chord_scaling_pos * te + (1 - self.chord_scaling_pos) * le - partials["mesh", "chord"] = (mesh - quarter_chord).flatten() + partials["mesh", "chord"] = (mesh - chord_pos).flatten() nx, ny, _ = mesh.shape nn = nx * ny * 3 @@ -187,12 +192,12 @@ def compute_partials(self, inputs, partials): d_qc = (np.einsum("ij,i->ij", np.ones((ny, 3)), 1.0 - chord_dist)).flatten() nnq = (nx - 1) * ny * 3 - partials["mesh", "in_mesh"][nn : nn + nnq] = np.tile(0.25 * d_qc, nx - 1) - partials["mesh", "in_mesh"][nn + nnq :] = np.tile(0.75 * d_qc, nx - 1) + partials["mesh", "in_mesh"][nn : nn + nnq] = np.tile(self.chord_scaling_pos * d_qc, nx - 1) + partials["mesh", "in_mesh"][nn + nnq :] = np.tile((1 - self.chord_scaling_pos) * d_qc, nx - 1) nnq = ny * 3 - partials["mesh", "in_mesh"][nn - nnq : nn] += 0.25 * d_qc - partials["mesh", "in_mesh"][:nnq] += 0.75 * d_qc + partials["mesh", "in_mesh"][nn - nnq : nn] += self.chord_scaling_pos * d_qc + partials["mesh", "in_mesh"][:nnq] += (1 - self.chord_scaling_pos) * d_qc class Sweep(om.ExplicitComponent): diff --git a/openaerostruct/geometry/tests/test_geometry_mesh.py b/openaerostruct/geometry/tests/test_geometry_mesh.py index 633cf45e3..002a68d91 100644 --- a/openaerostruct/geometry/tests/test_geometry_mesh.py +++ b/openaerostruct/geometry/tests/test_geometry_mesh.py @@ -21,6 +21,7 @@ def test(self): # The way this is currently set up, we don't actually use the values here surface["twist_cp"] = np.zeros((5)) surface["chord_cp"] = np.zeros((5)) + surface["chord_scaling_pos"] = 0.25 surface["xshear_cp"] = np.zeros((5)) surface["yshear_cp"] = np.zeros((5)) surface["zshear_cp"] = np.zeros((5)) diff --git a/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py b/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py index c2b0856f7..e75b1633b 100644 --- a/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py +++ b/openaerostruct/geometry/tests/test_geometry_mesh_transformations.py @@ -4,7 +4,7 @@ import unittest import openmdao.api as om -from openmdao.utils.assert_utils import assert_check_partials +from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal from openaerostruct.geometry.geometry_mesh_transformations import ( Taper, @@ -31,7 +31,13 @@ def get_mesh(symmetry): ny = (2 * NY - 1) if symmetry else NY # Create a dictionary to store options about the mesh - mesh_dict = {"num_y": ny, "num_x": NX, "wing_type": "CRM", "symmetry": symmetry, "num_twist_cp": NY} + mesh_dict = { + "num_y": ny, + "num_x": NX, + "wing_type": "CRM", + "symmetry": symmetry, + "num_twist_cp": NY, + } # Generate the aerodynamic mesh based on the previous dictionary mesh, twist_cp = generate_mesh(mesh_dict) @@ -132,6 +138,51 @@ def test_scalex_symmetry(self): check = prob.check_partials(compact_print=True, abs_err_tol=1e-5, rel_err_tol=1e-5) assert_check_partials(check, atol=1e-6, rtol=1e-6) + def test_scalex_chord_scaling_pos_random(self): + symmetry = False + mesh = get_mesh(symmetry) + + # Test for random values of chord_scaling_pos, check derivatives + prob = om.Problem() + group = prob.model + + val = self.rng.random(NY) + chord_scaling_pos = self.rng.random(1) + + comp = ScaleX(val=val, mesh_shape=mesh.shape, chord_scaling_pos=chord_scaling_pos) + group.add_subsystem("comp", comp) + + prob.setup() + + prob["comp.in_mesh"] = mesh + + prob.run_model() + + check = prob.check_partials(compact_print=True, abs_err_tol=1e-5, rel_err_tol=1e-5) + assert_check_partials(check, atol=1e-6, rtol=1e-6) + + def test_scalex_chord_scaling_pos_trailing_edge(self): + symmetry = True + mesh = get_mesh(symmetry) + + # Test for chord_scaling_pos at trailing edge + prob = om.Problem() + group = prob.model + + val = self.rng.random(NY) + + comp = ScaleX(val=val, mesh_shape=mesh.shape, chord_scaling_pos=1) + group.add_subsystem("comp", comp) + + prob.setup() + + prob["comp.in_mesh"] = mesh + + prob.run_model() + + # If chord_scaling_pos = 1, TE should not move + assert_near_equal(mesh[-1, :, :], prob["comp.mesh"][-1, :, :], tolerance=1e-10) + def test_sweep(self): symmetry = False mesh = get_mesh(symmetry)