Add the possibility to change the position of the scaling for the cho…

…rd design variable (#410)

* chord scaling position

* flake8

* flake8 balc

* git pb

* pb git

* add test chord scaling

* change documentation

* corrected test

* add warnings and check chord scaling

* stack level

* change version init

* Separate chord_scaling_pos test into two separate tests

---------

Co-authored-by: Eytan Adler <eytana@umich.edu>

openaerostruct/__init__.py

@@ -1 +1 @@
-__version__ = "2.6.1"
+__version__ = "2.6.2"

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

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

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
 

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,19 +171,19 @@ 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"]
         chord_dist = inputs["chord"]
 
         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):

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))

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)