From 2662803688565ef68579da44d5e0ea6a679b369e Mon Sep 17 00:00:00 2001 From: kanekosh Date: Tue, 20 May 2025 15:25:16 -0400 Subject: [PATCH 01/30] truss-braced structure model --- openaerostruct/structures/disp.py | 5 +- openaerostruct/structures/fem.py | 128 ++++--- openaerostruct/structures/fem_truss_braced.py | 339 ++++++++++++++++++ .../structures/spatial_beam_states.py | 64 +++- 4 files changed, 474 insertions(+), 62 deletions(-) create mode 100644 openaerostruct/structures/fem_truss_braced.py diff --git a/openaerostruct/structures/disp.py b/openaerostruct/structures/disp.py index 4a5a59b89..8d5f1e85f 100644 --- a/openaerostruct/structures/disp.py +++ b/openaerostruct/structures/disp.py @@ -9,14 +9,15 @@ class Disp(om.ExplicitComponent): a 2D array so we can more easily use the results. The solution to the linear system has meaingless entires due to the - constraints on the FEM model. The displacements from this portion of + boundary conditions on the FEM model. The "displacements" from this portion of the linear system are not needed, so we select only the relevant portion of the displacements for further calculations. Parameters ---------- disp_aug[6*(ny+1)] : numpy array - Augmented displacement array. Obtained by solving the system + Augmented displacement array with additional 6 Lagrange multipliers + for clamp boundary conditions at the end. Obtained by solving the system K * disp_aug = forces, where forces is a flattened version of loads. Returns diff --git a/openaerostruct/structures/fem.py b/openaerostruct/structures/fem.py index ef962e33c..4c6492875 100644 --- a/openaerostruct/structures/fem.py +++ b/openaerostruct/structures/fem.py @@ -7,6 +7,74 @@ import openmdao.api as om +def get_drdu_sparsity_pattern(ny, vec_size, symmetry): + """ + return sparsity pattern of partials of residuals w.r.t. displacements = stiffness matrix + """ + base_row = np.repeat(0, 6) + base_col = np.arange(6) + + # Upper diagonal blocks + rows1 = np.tile(base_row, 6 * (ny - 1)) + np.repeat(np.arange(6 * (ny - 1)), 6) + col = np.tile(base_col + 6, 6) + cols1 = np.tile(col, ny - 1) + np.repeat(6 * np.arange(ny - 1), 36) + + # Lower diagonal blocks + rows2 = np.tile(base_row + 6, 6 * (ny - 1)) + np.repeat(np.arange(6 * (ny - 1)), 6) + col = np.tile(base_col, 6) + cols2 = np.tile(col, ny - 1) + np.repeat(6 * np.arange(ny - 1), 36) + + # Main diagonal blocks, root + rows3 = np.tile(base_row, 6) + np.repeat(np.arange(6), 6) + cols3 = np.tile(base_col, 6) + + # Main diagonal blocks, tip + rows4 = np.tile(base_row + (ny - 1) * 6, 6) + np.repeat(np.arange(6), 6) + cols4 = np.tile(base_col + (ny - 1) * 6, 6) + + # Main diagonal blocks, interior + rows5 = np.tile(base_row + 6, 6 * (ny - 2)) + np.repeat(np.arange(6 * (ny - 2)), 6) + col = np.tile(base_col + 6, 6) + cols5 = np.tile(col, ny - 2) + np.repeat(6 * np.arange(ny - 2), 36) + + # Find constrained nodes based on closeness to specified cg point + if symmetry: + idx = ny - 1 + else: + idx = (ny - 1) // 2 + + index = 6 * idx + num_dofs = 6 * ny + arange = np.arange(6) + + # Fixed boundary condition. + rows6 = index + arange + cols6 = num_dofs + arange + + rows = np.concatenate([rows1, rows2, rows3, rows4, rows5, rows6, cols6]) + cols = np.concatenate([cols1, cols2, cols3, cols4, cols5, cols6, rows6]) + + # vectorize for multiple RHS + sp_size = len(rows) + vec_rows = np.tile(rows, vec_size) + np.repeat(sp_size * np.arange(vec_size), sp_size) + vec_cols = np.tile(cols, vec_size) + np.repeat(sp_size * np.arange(vec_size), sp_size) + + return rows, cols, vec_rows, vec_cols + + +def get_drdK_sparsity_pattern(ny): + """ + return sparsity pattern of partials of residuals w.r.t. stiffness matrix + """ + base_row = np.tile(0, 12) + base_col = np.arange(12) + row = np.tile(base_row, 12) + np.repeat(np.arange(12), 12) + col = np.tile(base_col, 12) + np.repeat(12 * np.arange(12), 12) + rows = np.tile(row, ny - 1) + np.repeat(6 * np.arange(ny - 1), 144) + cols = np.tile(col, ny - 1) + np.repeat(144 * np.arange(ny - 1), 144) + return rows, cols + + class FEM(om.ImplicitComponent): """ Component that solves a linear system, Ax=b. @@ -78,64 +146,12 @@ def setup(self): # The derivative of residual wrt displacements is the stiffness matrix K. We can use the # sparsity pattern here and when constucting the sparse matrix, so save rows and cols. - - base_row = np.repeat(0, 6) - base_col = np.arange(6) - - # Upper diagonal blocks - rows1 = np.tile(base_row, 6 * (ny - 1)) + np.repeat(np.arange(6 * (ny - 1)), 6) - col = np.tile(base_col + 6, 6) - cols1 = np.tile(col, ny - 1) + np.repeat(6 * np.arange(ny - 1), 36) - - # Lower diagonal blocks - rows2 = np.tile(base_row + 6, 6 * (ny - 1)) + np.repeat(np.arange(6 * (ny - 1)), 6) - col = np.tile(base_col, 6) - cols2 = np.tile(col, ny - 1) + np.repeat(6 * np.arange(ny - 1), 36) - - # Main diagonal blocks, root - rows3 = np.tile(base_row, 6) + np.repeat(np.arange(6), 6) - cols3 = np.tile(base_col, 6) - - # Main diagonal blocks, tip - rows4 = np.tile(base_row + (ny - 1) * 6, 6) + np.repeat(np.arange(6), 6) - cols4 = np.tile(base_col + (ny - 1) * 6, 6) - - # Main diagonal blocks, interior - rows5 = np.tile(base_row + 6, 6 * (ny - 2)) + np.repeat(np.arange(6 * (ny - 2)), 6) - col = np.tile(base_col + 6, 6) - cols5 = np.tile(col, ny - 2) + np.repeat(6 * np.arange(ny - 2), 36) - - # Find constrained nodes based on closeness to specified cg point - symmetry = self.options["surface"]["symmetry"] - if symmetry: - idx = self.ny - 1 - else: - idx = (self.ny - 1) // 2 - - index = 6 * idx - num_dofs = 6 * ny - arange = np.arange(6) - - # Fixed boundary condition. - rows6 = index + arange - cols6 = num_dofs + arange - - self.k_rows = rows = np.concatenate([rows1, rows2, rows3, rows4, rows5, rows6, cols6]) - self.k_cols = cols = np.concatenate([cols1, cols2, cols3, cols4, cols5, cols6, rows6]) - - sp_size = len(rows) - vec_rows = np.tile(rows, vec_size) + np.repeat(sp_size * np.arange(vec_size), sp_size) - vec_cols = np.tile(cols, vec_size) + np.repeat(sp_size * np.arange(vec_size), sp_size) - + rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, self.options["surface"]["symmetry"]) + self.k_rows = rows + self.k_cols = cols self.declare_partials(of="disp_aug", wrt="disp_aug", rows=vec_rows, cols=vec_cols) - base_row = np.tile(0, 12) - base_col = np.arange(12) - row = np.tile(base_row, 12) + np.repeat(np.arange(12), 12) - col = np.tile(base_col, 12) + np.repeat(12 * np.arange(12), 12) - rows = np.tile(row, ny - 1) + np.repeat(6 * np.arange(ny - 1), 144) - cols = np.tile(col, ny - 1) + np.repeat(144 * np.arange(ny - 1), 144) - + rows, cols = get_drdK_sparsity_pattern(ny) self.declare_partials("disp_aug", "local_stiff_transformed", rows=rows, cols=cols) def apply_nonlinear(self, inputs, outputs, residuals): diff --git a/openaerostruct/structures/fem_truss_braced.py b/openaerostruct/structures/fem_truss_braced.py new file mode 100644 index 000000000..716e81fb4 --- /dev/null +++ b/openaerostruct/structures/fem_truss_braced.py @@ -0,0 +1,339 @@ +"""Define the LinearSystemComp class.""" + +import numpy as np +from scipy.sparse import coo_matrix +from scipy.sparse.linalg import splu + +import openmdao.api as om + +from .fem import get_drdu_sparsity_pattern, get_drdK_sparsity_pattern + + +class FEMTrussBraced(om.ImplicitComponent): + """ + Component that solves an FEM linear system, K u = f. + This component is a customization of original fem.py to solve truss-braced wing structure. + + This solve the following linear system for the wing-truss joint problem: + + [ K1 | 0 | C1^T ] [u1] [f1] + [ 0 | K2 | C2^T ] [u2] = [f2] + [ C1 | C2 | 0 ] [lambda] [0 ] + + where: + K1 = the stiffness matrix of the wing + additional entries for wing root boundary conditions + u1 = augmented displacement of the wing (includes Lagrange multipliers for boundary conditions) + f1 = augmented force vector of the wing (includes zeros at the end for boundary conditions) + K2 = the stiffness matrix of the truss + additional entries for truss root boundary conditions + u2 = augmented displacement of the truss (includes Lagrange multipliers for boundary conditions) + f2 = augmented force vector of the truss (includes zeros at the end for boundary conditions) + C1, C2 = constraint matrix for the wing-truss joint + lambda = Lagrange multipliers for the joint constraints (physically these are joint reaction forces) + + Attributes + ---------- + _lup : None or list(object) + matrix factorizations returned from scipy.linag.lu_factor for each A matrix + k_cols : ndarray + Cached column indices for sparse representation of stiffness matrix. + k_rows : ndarray + Cached row indices for sparse representation of stiffness matrix. + k_data : ndarray + Cached values for sparse representation of stiffness matrix. + """ + + def __init__(self, **kwargs): + """ + Intialize the LinearSystemComp component. + + Parameters + ---------- + kwargs : dict of keyword arguments + Keyword arguments that will be mapped into the Component options. + """ + super(FEMTrussBraced, self).__init__(**kwargs) + self._lup = None + self.k_cols = None + self.k_rows = None + self.k_data = None + + def initialize(self): + """ + Declare options. + """ + self.options.declare("surfaces", types=list) + self.options.declare("vec_size", types=int, default=1, desc="Number of linear systems to solve.") + + def setup(self): + """ + Matrix and RHS are inputs, solution vector is the output. + """ + vec_size = self.options["vec_size"] + self.surfaces = self.options["surfaces"] + if len(self.surfaces) != 2: + raise ValueError("FEMTrussBraced component requires exactly two surfaces (wing and truss).") + self.surface_names = [self.surfaces[0]["name"], self.surfaces[1]["name"]] + self.ny = [] # number of spanwise nodes + self.size = [] # size of assembled K matrix for each surface + self.num_k_data = [] # number of non-zero entries for assembled K matrix for each surface + self._lup = [] + k_cols = [] # column indices for stiffness matrix sparsity pattern + k_rows = [] # row indices for stiffness matrix sparsity pattern + joint_indices = [] # node index for the joint between wing and truss + + for i, surface in enumerate(self.surfaces): + name = surface["name"] + + # --- inputs, outputs, partials of FEM system for each surface --- + ny = surface["mesh"].shape[1] + size = int(6 * ny + 6) + self.ny.append(ny) + self.size.append(size) + + full_size = size * vec_size + shape = (vec_size, size) if vec_size > 1 else (size,) + init_locK = np.tile(np.eye(12).flatten(), ny - 1).reshape(ny - 1, 12, 12) + self.add_input(f"local_stiff_transformed_{name}", val=init_locK) + self.add_input(f"forces_{name}", val=np.ones(shape), units="N") + self.add_output(f"disp_aug_{name}", shape=shape, val=0.1, units="m") + + # Set up the derivatives. + row_col = np.arange(full_size, dtype="int") + self.declare_partials(f"disp_aug_{name}", f"forces_{name}", val=np.full(full_size, -1.0), rows=row_col, cols=row_col) + + # The derivative of residual wrt displacements is the stiffness matrix K. We can use the + # sparsity pattern here and when constucting the sparse matrix, so save rows and cols. + rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"]) + if i == 1: + # the second surface (truss) goes to the lower-right block of the entire stiffness matrix + rows += self.size[0] + cols += self.size[0] + k_rows.append(rows) + k_cols.append(cols) + self.num_k_data.append(len(rows)) + self.declare_partials(of=f"disp_aug_{name}", wrt=f"disp_aug_{name}", rows=vec_rows, cols=vec_cols) + + rows, cols = get_drdK_sparsity_pattern(ny) + self.declare_partials(f"disp_aug_{name}", f"local_stiff_transformed_{name}", rows=rows, cols=cols) + + # --- joint constraints (coupling terms between two surfaces) --- + joint_type = self.surfaces[0]["joint_type"] + if joint_type == 'pin': + self.n_con = n_con = 3 # number of joint constraints: translational only + elif joint_type == 'rigid': + self.n_con = n_con = 6 # number of joint constraints: translational and rotational + else: + raise ValueError("Joint type must be either 'pin' or 'rigid'.") + + # add Lagrange multipliers for joint constraints as state variables + self.add_output("joint_Lag", shape=(n_con,), val=0.0) + + k_rows_joint = [] + k_cols_joint = [] + self.k_data_joint = [] + for i, surface in enumerate(self.surfaces): + name = surface["name"] + + # find the node index for the joint + joint_y = surface["joint_y"] + if not surface["symmetry"]: + raise NotImplementedError("Truss-braced wing joint is only implemented for symmetric surfaces.") + # check sign convention for y (spanwise coordinate) + if joint_y * surface["mesh"][0, 0, 1] < 0: + joint_y *= -1 + + joint_idx = np.argmin(np.abs(surface["mesh"][0, :, 1] - joint_y)) + joint_indices.append(joint_idx) + print("Surface", name, "joint y:", surface["mesh"][0, joint_idx, 1], "joint index:", joint_idx) + + # partials of joint constraint residuals w.r.t. displacement + rows = np.arange(n_con) + cols = np.arange(n_con) + 6 * joint_idx + if i == 0: + vals = np.ones(n_con) * 1e9 + else: + vals = np.ones(n_con) * -1e9 + self.declare_partials("joint_Lag", f"disp_aug_{name}", rows=rows, cols=cols, val=vals) + + # partials of FEM residuals (r = Ku - f) w.r.t. joint Lagrange multipliers: transpose of above + self.declare_partials(f"disp_aug_{name}", "joint_Lag", rows=cols, cols=rows, val=vals) + + # append these entries to the bottom and right of the entire K matrix + rows += self.size[0] + self.size[1] + if i == 1: + cols += self.size[0] + + # partials of joint constraint residuals w.r.t. local stiffness matrix + k_rows_joint.append(rows) + k_cols_joint.append(cols) + self.k_data_joint.append(vals) + # partials of FEM residuals wrt joint Lagrange multipliers + k_rows_joint.append(cols) + k_cols_joint.append(rows) + self.k_data_joint.append(vals) + + # row and col indices for the entire stiffness matrix (that combines wing, truss, and joint constraints) + self.k_cols = np.concatenate(k_cols + k_cols_joint) + self.k_rows = np.concatenate(k_rows + k_rows_joint) + self.total_size = sum(self.size) + n_con # size of total stiffness matrix (wing + truss + constraints) + + def apply_nonlinear(self, inputs, outputs, residuals): + """ + R = Ax - b. + + Parameters + ---------- + inputs : Vector + unscaled, dimensional input variables read via inputs[key] + outputs : Vector + unscaled, dimensional output variables read via outputs[key] + residuals : Vector + unscaled, dimensional residuals written to via residuals[key] + """ + K = self.assemble_CSC_K(inputs) + disp = np.concatenate([outputs[f"disp_aug_{name}"] for name in self.surface_names]) + joint_Lag = outputs["joint_Lag"] + u = np.concatenate([disp, joint_Lag]) + force = np.concatenate([inputs[f"forces_{name}"] for name in self.surface_names]) + rhs = np.concatenate([force, np.zeros(self.n_con)]) + + r = K.dot(u) - rhs + + size0, size1 = self.size[0], self.size[1] # size of each residuals + residuals[f"disp_aug_{self.surface_names[0]}"] = r[:size0] + residuals[f"disp_aug_{self.surface_names[1]}"] = r[size0:size0 + size1] + residuals["joint_Lag"] = r[size0 + size1:] + + def solve_nonlinear(self, inputs, outputs): + """ + Use numpy to solve Ax=b for x. + + Parameters + ---------- + inputs : Vector + unscaled, dimensional input variables read via inputs[key] + outputs : Vector + unscaled, dimensional output variables read via outputs[key] + """ + # lu factorization for use with solve_linear + K = self.assemble_CSC_K(inputs) + self._lup = splu(K) + force = np.concatenate([inputs[f"forces_{name}"] for name in self.surface_names]) + rhs = np.concatenate([force, np.zeros(self.n_con)]) + + u = self._lup.solve(rhs) + + size0, size1 = self.size[0], self.size[1] # size of each displacement vectors + outputs[f"disp_aug_{self.surface_names[0]}"] = u[:size0] + outputs[f"disp_aug_{self.surface_names[1]}"] = u[size0:size0 + size1] + outputs["joint_Lag"] = u[size0 + size1:] + + def linearize(self, inputs, outputs, J): + """ + Compute the non-constant partial derivatives. + + Parameters + ---------- + inputs : Vector + unscaled, dimensional input variables read via inputs[key] + outputs : Vector + unscaled, dimensional output variables read via outputs[key] + J : Jacobian + sub-jac components written to jacobian[output_name, input_name] + """ + vec_size = self.options["vec_size"] + + name0, name1 = self.surface_names[0], self.surface_names[1] + ny0, ny1 = self.ny[0], self.ny[1] + idx0 = np.tile(np.tile(np.arange(12), 12), ny0 - 1) + np.repeat(6 * np.arange(ny0 - 1), 144) + idx1 = np.tile(np.tile(np.arange(12), 12), ny1 - 1) + np.repeat(6 * np.arange(ny1 - 1), 144) + disp0 = outputs[f"disp_aug_{name0}"] + disp1 = outputs[f"disp_aug_{name1}"] + J[f"disp_aug_{name0}", f"local_stiff_transformed_{name0}"] = np.tile(disp0[idx0], vec_size) + J[f"disp_aug_{name1}", f"local_stiff_transformed_{name1}"] = np.tile(disp1[idx1], vec_size) + + nk0, nk1 = self.num_k_data[0], self.num_k_data[1] + k_data0 = self.k_data[:nk0] + k_data1 = self.k_data[nk0:nk0 + nk1] + J[f"disp_aug_{name0}", f"disp_aug_{name0}"] = np.tile(k_data0, vec_size) + J[f"disp_aug_{name1}", f"disp_aug_{name1}"] = np.tile(k_data1, vec_size) + + def solve_linear(self, d_outputs, d_residuals, mode): + r""" + Back-substitution to solve the derivatives of the linear system. + + If mode is: + 'fwd': d_residuals \|-> d_outputs + + 'rev': d_outputs \|-> d_residuals + + Parameters + ---------- + d_outputs : Vector + unscaled, dimensional quantities read via d_outputs[key] + d_residuals : Vector + unscaled, dimensional quantities read via d_residuals[key] + mode : str + either 'fwd' or 'rev' + """ + vec_size = self.options["vec_size"] + size0, size1 = self.size[0], self.size[1] + name0, name1 = self.surface_names[0], self.surface_names[1] + + # TODO: verify (check_partials doesn't cover this) + + if mode == "fwd": + if vec_size > 1: + for j in range(vec_size): + rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"][j], d_residuals[f"disp_aug_{name1}"][j], d_residuals["joint_Lag"][j])) + sol = self._lup.solve(rhs) + d_outputs[f"disp_aug_{name0}"][j] = sol[:size0] + d_outputs[f"disp_aug_{name1}"][j] = sol[size0:size0 + size1] + d_outputs["joint_Lag"][j] = sol[size0 + size1:] + else: + rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"], d_residuals[f"disp_aug_{name1}"], d_residuals["joint_Lag"])) + sol = self._lup.solve(rhs) + d_outputs[f"disp_aug_{name0}"] = sol[:size0] + d_outputs[f"disp_aug_{name1}"] = sol[size0:size0 + size1] + d_outputs["joint_Lag"] = sol[size0 + size1:] + else: + if vec_size > 1: + for j in range(vec_size): + rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"][j], d_outputs[f"disp_aug_{name1}"][j], d_outputs["joint_Lag"][j])) + sol = self._lup.solve(rhs) + d_residuals[f"disp_aug_{name0}"][j] = sol[:size0] + d_residuals[f"disp_aug_{name1}"][j] = sol[size0:size0 + size1] + d_residuals["joint_Lag"][j] = sol[size0 + size1:] + else: + rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"], d_outputs[f"disp_aug_{name1}"], d_outputs["joint_Lag"])) + sol = self._lup.solve(rhs) + d_residuals[f"disp_aug_{name0}"] = sol[:size0] + d_residuals[f"disp_aug_{name1}"] = sol[size0:size0 + size1] + d_residuals["joint_Lag"] = sol[size0 + size1:] + + def assemble_CSC_K(self, inputs): + """ + Assemble the stiffness matrix in sparse CSC format. + + Returns + ------- + ndarray + Stiffness matrix as dense ndarray. + """ + data = {} + for surface in self.surfaces: + name = surface["name"] + k_loc = inputs[f"local_stiff_transformed_{name}"] + + data1 = k_loc[:, :6, 6:].flatten() + data2 = k_loc[:, 6:, :6].flatten() + data3 = k_loc[0, :6, :6].flatten() + data4 = k_loc[-1, 6:, 6:].flatten() + data5 = (k_loc[0:-1, 6:, 6:] + k_loc[1:, :6, :6]).flatten() + data6 = np.full((6,), 1e9) + data[name] = np.concatenate([data1, data2, data3, data4, data5, data6, data6]) + + # combine matrices for wing, truss, and joint constraints + self.k_data = np.concatenate(list(data.values()) + self.k_data_joint) + + return coo_matrix((self.k_data, (self.k_rows, self.k_cols)), shape=(self.total_size, self.total_size)).tocsc() diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index 4206d619b..13f0ef4cf 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -1,6 +1,7 @@ import openmdao.api as om from openaerostruct.structures.create_rhs import CreateRHS from openaerostruct.structures.fem import FEM +from openaerostruct.structures.fem_truss_braced import FEMTrussBraced from openaerostruct.structures.disp import Disp from openaerostruct.structures.wing_weight_loads import StructureWeightLoads from openaerostruct.structures.fuel_loads import FuelLoads @@ -9,8 +10,35 @@ from openaerostruct.structures.compute_thrust_loads import ComputeThrustLoads -class SpatialBeamStates(om.Group): - """Group that contains the spatial beam states.""" +class FEMloads(om.Group): + """ + Compute FEM right-hand side force vector + + Parameters + ---------- + nodes : ndarray + The coordinates of the FEM nodes. + load_factor : float + The load factor applied to the structure. + element_mass : ndarray + The mass of each element in the structure for bending relief. + fuel_mass : ndarray + Total fuel mass for bendling relief + fuel_vols : ndarray + The feul volume of each element + point_mass_locations : ndarray + The locations of the point masses. + point_masses : ndarray + Point masses as point load + engine_thrusts : ndarray + Engine thrusts as point load + + Returns + ------- + forces : ndarray + The right-hand side force vector for the FEM. This is a flattened array of loads + plus zeros at the end for boundary condition + """ def initialize(self): self.options.declare("surface", types=dict) @@ -65,6 +93,34 @@ def setup(self): "create_rhs", CreateRHS(surface=surface), promotes_inputs=["total_loads"], promotes_outputs=["forces"] ) - self.add_subsystem("fem", FEM(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) - self.add_subsystem("disp", Disp(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) +class SpatialBeamStates(om.Group): + """Group that contains the spatial beam states.""" + + def initialize(self): + self.options.declare("surface", types=(dict, list)) + self.options.declare("truss_braced", default=False, types=bool) + + def setup(self): + surface = self.options["surface"] + + if not self.options["truss_braced"]: + self.add_subsystem('forces', FEMloads(surface=surface), promotes_inputs=["*"], promotes_outputs=["forces"]) + self.add_subsystem("fem", FEM(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) + self.add_subsystem("disp", Disp(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) + + else: + # compute RHS force vector for wing and truss, respectively + for surf in surface: + name = surf["name"] + self.add_subsystem(f'forces_{name}', FEMloads(surface=surf)) + self.connect(f"forces_{name}.forces", f"fem.forces_{name}") + + # FEM of wing-truss system + self.add_subsystem("fem", FEMTrussBraced(surfaces=surface), promotes_inputs=["local_stiff_transformed_*"]) + + # reshape displacements and remove Lagrange multipliers from disp_aug + for surf in surface: + name = surf["name"] + self.add_subsystem(f"disp_{name}", Disp(surface=surf), promotes_outputs=[("disp", f"disp_{name}")]) + self.connect(f"fem.disp_aug_{name}", f"disp_{name}.disp_aug") From ed9dcbbe3256073267be93ced4a73f7ee2f1e5d2 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Tue, 20 May 2025 18:13:42 -0400 Subject: [PATCH 02/30] renamed truss to strut --- ...em_truss_braced.py => fem_strut_braced.py} | 34 +++++++++++-------- .../structures/spatial_beam_states.py | 12 +++---- 2 files changed, 25 insertions(+), 21 deletions(-) rename openaerostruct/structures/{fem_truss_braced.py => fem_strut_braced.py} (93%) diff --git a/openaerostruct/structures/fem_truss_braced.py b/openaerostruct/structures/fem_strut_braced.py similarity index 93% rename from openaerostruct/structures/fem_truss_braced.py rename to openaerostruct/structures/fem_strut_braced.py index 716e81fb4..6487a3610 100644 --- a/openaerostruct/structures/fem_truss_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -8,13 +8,17 @@ from .fem import get_drdu_sparsity_pattern, get_drdK_sparsity_pattern +# TODO: +# - implement ball and pin joint (current pin joint should be renamed to ball) +# - implement pin/ball joint for the root symmetric boundary -class FEMTrussBraced(om.ImplicitComponent): + +class FEMStrutBraced(om.ImplicitComponent): """ Component that solves an FEM linear system, K u = f. - This component is a customization of original fem.py to solve truss-braced wing structure. + This component is a customization of original fem.py to solve strut-braced wing structure. - This solve the following linear system for the wing-truss joint problem: + This solve the following linear system for the wing-strut joint problem: [ K1 | 0 | C1^T ] [u1] [f1] [ 0 | K2 | C2^T ] [u2] = [f2] @@ -24,10 +28,10 @@ class FEMTrussBraced(om.ImplicitComponent): K1 = the stiffness matrix of the wing + additional entries for wing root boundary conditions u1 = augmented displacement of the wing (includes Lagrange multipliers for boundary conditions) f1 = augmented force vector of the wing (includes zeros at the end for boundary conditions) - K2 = the stiffness matrix of the truss + additional entries for truss root boundary conditions - u2 = augmented displacement of the truss (includes Lagrange multipliers for boundary conditions) - f2 = augmented force vector of the truss (includes zeros at the end for boundary conditions) - C1, C2 = constraint matrix for the wing-truss joint + K2 = the stiffness matrix of the strut + additional entries for strut root boundary conditions + u2 = augmented displacement of the strut (includes Lagrange multipliers for boundary conditions) + f2 = augmented force vector of the strut (includes zeros at the end for boundary conditions) + C1, C2 = constraint matrix for the wing-strut joint lambda = Lagrange multipliers for the joint constraints (physically these are joint reaction forces) Attributes @@ -51,7 +55,7 @@ def __init__(self, **kwargs): kwargs : dict of keyword arguments Keyword arguments that will be mapped into the Component options. """ - super(FEMTrussBraced, self).__init__(**kwargs) + super(FEMStrutBraced, self).__init__(**kwargs) self._lup = None self.k_cols = None self.k_rows = None @@ -71,7 +75,7 @@ def setup(self): vec_size = self.options["vec_size"] self.surfaces = self.options["surfaces"] if len(self.surfaces) != 2: - raise ValueError("FEMTrussBraced component requires exactly two surfaces (wing and truss).") + raise ValueError("FEMStrutBraced component requires exactly two surfaces (wing and strut).") self.surface_names = [self.surfaces[0]["name"], self.surfaces[1]["name"]] self.ny = [] # number of spanwise nodes self.size = [] # size of assembled K matrix for each surface @@ -79,7 +83,7 @@ def setup(self): self._lup = [] k_cols = [] # column indices for stiffness matrix sparsity pattern k_rows = [] # row indices for stiffness matrix sparsity pattern - joint_indices = [] # node index for the joint between wing and truss + joint_indices = [] # node index for the joint between wing and strut for i, surface in enumerate(self.surfaces): name = surface["name"] @@ -105,7 +109,7 @@ def setup(self): # sparsity pattern here and when constucting the sparse matrix, so save rows and cols. rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"]) if i == 1: - # the second surface (truss) goes to the lower-right block of the entire stiffness matrix + # the second surface (strut) goes to the lower-right block of the entire stiffness matrix rows += self.size[0] cols += self.size[0] k_rows.append(rows) @@ -137,7 +141,7 @@ def setup(self): # find the node index for the joint joint_y = surface["joint_y"] if not surface["symmetry"]: - raise NotImplementedError("Truss-braced wing joint is only implemented for symmetric surfaces.") + raise NotImplementedError("Strut-braced wing joint is only implemented for symmetric surfaces.") # check sign convention for y (spanwise coordinate) if joint_y * surface["mesh"][0, 0, 1] < 0: joint_y *= -1 @@ -172,10 +176,10 @@ def setup(self): k_cols_joint.append(rows) self.k_data_joint.append(vals) - # row and col indices for the entire stiffness matrix (that combines wing, truss, and joint constraints) + # row and col indices for the entire stiffness matrix (that combines wing, strut, and joint constraints) self.k_cols = np.concatenate(k_cols + k_cols_joint) self.k_rows = np.concatenate(k_rows + k_rows_joint) - self.total_size = sum(self.size) + n_con # size of total stiffness matrix (wing + truss + constraints) + self.total_size = sum(self.size) + n_con # size of total stiffness matrix (wing + strut + constraints) def apply_nonlinear(self, inputs, outputs, residuals): """ @@ -333,7 +337,7 @@ def assemble_CSC_K(self, inputs): data6 = np.full((6,), 1e9) data[name] = np.concatenate([data1, data2, data3, data4, data5, data6, data6]) - # combine matrices for wing, truss, and joint constraints + # combine matrices for wing, strut, and joint constraints self.k_data = np.concatenate(list(data.values()) + self.k_data_joint) return coo_matrix((self.k_data, (self.k_rows, self.k_cols)), shape=(self.total_size, self.total_size)).tocsc() diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index 13f0ef4cf..c2c63b2cd 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -1,7 +1,7 @@ import openmdao.api as om from openaerostruct.structures.create_rhs import CreateRHS from openaerostruct.structures.fem import FEM -from openaerostruct.structures.fem_truss_braced import FEMTrussBraced +from openaerostruct.structures.fem_strut_braced import FEMStrutBraced from openaerostruct.structures.disp import Disp from openaerostruct.structures.wing_weight_loads import StructureWeightLoads from openaerostruct.structures.fuel_loads import FuelLoads @@ -99,25 +99,25 @@ class SpatialBeamStates(om.Group): def initialize(self): self.options.declare("surface", types=(dict, list)) - self.options.declare("truss_braced", default=False, types=bool) + self.options.declare("strut_braced", default=False, types=bool) def setup(self): surface = self.options["surface"] - if not self.options["truss_braced"]: + if not self.options["strut_braced"]: self.add_subsystem('forces', FEMloads(surface=surface), promotes_inputs=["*"], promotes_outputs=["forces"]) self.add_subsystem("fem", FEM(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) self.add_subsystem("disp", Disp(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) else: - # compute RHS force vector for wing and truss, respectively + # compute RHS force vector for wing and strut, respectively for surf in surface: name = surf["name"] self.add_subsystem(f'forces_{name}', FEMloads(surface=surf)) self.connect(f"forces_{name}.forces", f"fem.forces_{name}") - # FEM of wing-truss system - self.add_subsystem("fem", FEMTrussBraced(surfaces=surface), promotes_inputs=["local_stiff_transformed_*"]) + # FEM of wing-strut system + self.add_subsystem("fem", FEMStrutBraced(surfaces=surface), promotes_inputs=["local_stiff_transformed_*"]) # reshape displacements and remove Lagrange multipliers from disp_aug for surf in surface: From 4e0c5bacf26f50c80f0198cca40837d4e1ae1d23 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Tue, 20 May 2025 23:53:29 -0400 Subject: [PATCH 03/30] modified AS group for strut-braced system --- .../integration/aerostruct_groups.py | 42 ++++++++++++++----- .../structures/spatial_beam_states.py | 2 +- 2 files changed, 33 insertions(+), 11 deletions(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 7eb1278be..1bd6de487 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -104,8 +104,10 @@ def setup(self): class CoupledAS(om.Group): + """ FEM + displacement transfer + VLM geometry property computations for single surface """ def initialize(self): self.options.declare("surface", types=dict) + self.options.declare("include_FEM", default=True, types=bool, desc="If False, skip the FEM subsystem under this group; FEM should be added external to this group") def setup(self): surface = self.options["surface"] @@ -120,12 +122,13 @@ def setup(self): set(["point_mass_locations", "point_masses", "nodes", "load_factor", "engine_thrusts"]) ) - self.add_subsystem( - "struct_states", - SpatialBeamStates(surface=surface), - promotes_inputs=["local_stiff_transformed", "forces", "loads"] + promotes, - promotes_outputs=["disp"], - ) + if self.options["include_FEM"]: + self.add_subsystem( + "struct_states", + SpatialBeamStates(surface=surface), + promotes_inputs=["local_stiff_transformed", "forces", "loads"] + promotes, + promotes_outputs=["disp"], + ) self.add_subsystem( "def_mesh", @@ -215,6 +218,7 @@ def initialize(self): self.options.declare( "rotational", False, types=bool, desc="Set to True to turn on support for computing angular velocities" ) + self.options.declare("strut_braced", default=False, types=bool) def setup(self): surfaces = self.options["surfaces"] @@ -222,13 +226,29 @@ def setup(self): coupled = om.Group() + if self.options["strut_braced"]: + if len(surfaces) != 2: + raise NameError("Strut-braced wing requires exactly two surfaces.") + # add FEM directly under the coupled group because multiple surfaces are coupled in the FEM level + coupled.add_subsystem( + "FEM_SBW", + SpatialBeamStates(surface=surfaces, strut_braced=True), + promotes_inputs=["*"], + promotes_outputs=["disp_*"] + ) + for surface in surfaces: name = surface["name"] # Connect the output of the loads component with the FEM # displacement parameter. This links the coupling within the coupled # group that necessitates the subgroup solver. - coupled.connect(name + "_loads.loads", name + ".loads") + loads_target = f"forces_{name}.loads" if self.options["strut_braced"] else f"{name}.loads" + coupled.connect(name + "_loads.loads", loads_target) + + # Connect displacement from FEM to transfer model. This is automatically done by promotion for non-SBW cases + if self.options["strut_braced"]: + coupled.connect(f"disp_{name}", f"{name}.disp") # Perform the connections with the modified names within the # 'aero_states' group. @@ -250,7 +270,8 @@ def setup(self): # Connect parameters from the 'coupled' group to the performance # groups for the individual surfaces. - self.connect("coupled." + name + ".disp", name + "_perf.disp") + disp_source = f"coupled.disp_{name}" if self.options["strut_braced"] else f"coupled.{name}.disp" + self.connect(disp_source, name + "_perf.disp") self.connect("coupled." + name + ".S_ref", name + "_perf.S_ref") self.connect("coupled." + name + ".widths", name + "_perf.widths") # self.connect('coupled.' + name + '.chords', name + '_perf.chords') @@ -266,13 +287,14 @@ def setup(self): # Add components to the 'coupled' group for each surface. # The 'coupled' group must contain all components and parameters # needed to converge the aerostructural system. - coupled_AS_group = CoupledAS(surface=surface) + include_FEM = False if self.options["strut_braced"] else True # For strut-braced case, FEM is already added outside of coupled AS + coupled_AS_group = CoupledAS(surface=surface, include_FEM=include_FEM) if ( surface["distributed_fuel_weight"] or "n_point_masses" in surface.keys() or surface["struct_weight_relief"] - ): + ) and not self.options["strut_braced"]: prom_in = ["load_factor"] else: prom_in = [] diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index c2c63b2cd..95157a7dc 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -113,7 +113,7 @@ def setup(self): # compute RHS force vector for wing and strut, respectively for surf in surface: name = surf["name"] - self.add_subsystem(f'forces_{name}', FEMloads(surface=surf)) + self.add_subsystem(f'forces_{name}', FEMloads(surface=surf), promotes_inputs=["load_factor"]) self.connect(f"forces_{name}.forces", f"fem.forces_{name}") # FEM of wing-strut system From f067b8999aa6f403d7438c1cd8427bf457715df8 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Thu, 3 Jul 2025 17:57:13 +0900 Subject: [PATCH 04/30] added buckling components --- .../integration/aerostruct_groups.py | 21 ++- .../structures/failure_buckling_ks.py | 128 ++++++++++++++++++ .../structures/spatial_beam_functionals.py | 14 ++ openaerostruct/structures/vonmises_wingbox.py | 28 ++++ 4 files changed, 178 insertions(+), 13 deletions(-) create mode 100644 openaerostruct/structures/failure_buckling_ks.py diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 1bd6de487..37d7c8cdc 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -183,22 +183,17 @@ def setup(self): ) elif surface["fem_model_type"] == "wingbox": + promotes_inputs = ["Qz", "J", "A_enc", "spar_thickness", "htop", "hbottom", "hfront", "hrear", "nodes", "disp"] + promotes_outputs = ["vonmises", "failure"] + if "buckling" in surface and surface["buckling"]: + promotes_inputs += ["skin_thickness", "t_over_c", "fem_chords"] + promotes_outputs += ["failure_buckling"] + self.add_subsystem( "struct_funcs", SpatialBeamFunctionals(surface=surface), - promotes_inputs=[ - "Qz", - "J", - "A_enc", - "spar_thickness", - "htop", - "hbottom", - "hfront", - "hrear", - "nodes", - "disp", - ], - promotes_outputs=["vonmises", "failure"], + promotes_inputs=promotes_inputs, + promotes_outputs=promotes_outputs, ) else: raise NameError("Please select a valid `fem_model_type` from either `tube` or `wingbox`.") diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py new file mode 100644 index 000000000..d723ea739 --- /dev/null +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -0,0 +1,128 @@ +import numpy as np +import openmdao.api as om + + +class FailureBucklingKS(om.ExplicitComponent): + """ + Buckling stress constraints for skin compression and spar shear buckling. + Returns KS-aggregated failure metric that should be constrained to be <= 0. + + Source: Michael Niu, Airframe Stress Analysis and Sizing, Hong Kong Conmilit Press Ltd., 1997 + + The rho inputeter controls how conservatively the KS function aggregates + the failure constraints. A lower value is more conservative while a greater + value is more aggressive (closer approximation to the max() function). + + parameters + ---------- + upper_skin_comp_stress[ny-1] : numpy array + Compression stress in the upper skin for each FEM element. + lower_skin_comp_stress[ny-1] : numpy array + Compression stress in the lower skin for each FEM element. + front_spar_shear_stress[ny-1] : numpy array + Shear stress in the front spar for each FEM element. + rear_spar_shear_stress[ny-1] : numpy array + Shear stress in the rear spar for each FEM element. + + Returns + ------- + failure_buckling : float + KS aggregation quantity by combining the failure constraints from upper skin buckling, + front spar shear buckling, and rear spar shear buckling for all FEM elements. + upper_skin_buckling_margin[ny-1] : numpy array + Margin for upper skin buckling for each FEM element, should be >= 0 + Margin = (buckling_limit - stress) / buckling_limit + front_spar_buckling_margin[ny-1] : numpy array + Margin for front spar shear buckling for each FEM element, should be >= 0 + rear_spar_buckling_margin[ny-1] : numpy array + Margin for rear spar shear buckling for each FEM element, should be >= 0 + + """ + + def initialize(self): + self.options.declare("surface", types=dict) + self.options.declare("rho", types=float, default=100.0) + + def setup(self): + surface = self.options["surface"] + rho = self.options["rho"] + + if surface["fem_model_type"] == "tube": + raise NotImplementedError("FailureBucklingKS is not implemented for tube FEM model type") + + self.ny = surface["mesh"].shape[1] + + # stress of each element + self.add_input("upper_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_input("lower_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_input("front_spar_shear_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_input("rear_spar_shear_stress", val=np.zeros(self.ny - 1), units="N/m**2") + # wingbox geometry + self.add_input("skin_thickness", val=np.zeros(self.ny - 1), units="m") + self.add_input("spar_thickness", val=np.zeros(self.ny - 1), units="m") + self.add_input("t_over_c", val=np.ones((self.ny - 1))) + self.add_input("fem_chords", val=np.zeros(self.ny - 1), units="m") + + self.add_output("failure_buckling", val=0.0) + # also output margins for each element/buckling for post-processing + # but don't declare partials for these outputs because we only use aggregated failure + self.add_output("upper_skin_buckling_margin", val=np.zeros(self.ny - 1), desc='(buckling_limit - stress) / buckling_limit >= 0?') + self.add_output("lower_skin_buckling_margin", val=np.zeros(self.ny - 1), desc='(buckling_limit - stress) / buckling_limit >= 0?') + self.add_output("front_spar_buckling_margin", val=np.zeros(self.ny - 1), desc='(buckling_limit - stress) / buckling_limit >= 0?') + self.add_output("rear_spar_buckling_margin", val=np.zeros(self.ny - 1), desc='(buckling_limit - stress) / buckling_limit >= 0?') + + self.rho = rho + + self.declare_partials("failure_buckling", "*", method="cs") + + def compute(self, inputs, outputs): + surface = self.options["surface"] + E = self.options["surface"]["E"] + rib_pitch = self.options["surface"]["rib_pitch"] + rho = self.options["rho"] # KS aggregation parameter + + # --- compute buckling stress --- + # buckling model parameters + eff_factor = 0.88 # stiffened panel efficiency for skin compression buckling. From Fig. 14.4.4 in Niu 1997 + Ks = 10.0 # model parameter for spar shear buckling. From Eq. 11.3.4 in Niu 1997 + + # buckling stress for skin compression. Eq. 14.4.2, Page 618 from Niu 1997 + # Positive stress = compression. + load_intensity_upper = inputs["skin_thickness"] * inputs["upper_skin_comp_stress"] # Pa*m + load_positive_upper = np.maximum(load_intensity_upper, 1.) # NOTE: this is C1 discontinuous so may negatively affect convergence + sigma_max_upper = eff_factor * (load_positive_upper * E / rib_pitch)**0.5 + load_intensity_lower = inputs["skin_thickness"] * inputs["lower_skin_comp_stress"] + load_positive_lower = np.maximum(load_intensity_lower, 1.) + sigma_max_lower = eff_factor * (load_positive_lower * E / rib_pitch)**0.5 + + # shear buckling stress of spars. Eq. 14.4.3, Page 618 from Niu 1997 + t_over_c_orig = surface["original_wingbox_airfoil_t_over_c"] + spar_height_front_nondim = (surface["data_y_upper"][0] - surface["data_y_lower"][0]) # nondim height from wingbox cross-section + spar_height_rear_nondim = (surface["data_y_upper"][-1] - surface["data_y_lower"][-1]) + spar_height_front = spar_height_front_nondim * inputs["fem_chords"] * (inputs["t_over_c"] / t_over_c_orig) + spar_height_rear = spar_height_rear_nondim * inputs["fem_chords"] * (inputs["t_over_c"] / t_over_c_orig) + tau_max_front = Ks * E * (inputs["spar_thickness"] / spar_height_front)**2 + tau_max_rear = Ks * E * (inputs["spar_thickness"] / spar_height_rear)**2 + + # --- compute failure and stress margins --- + outputs["upper_skin_buckling_margin"] = (sigma_max_upper - inputs["upper_skin_comp_stress"]) / sigma_max_upper + outputs["lower_skin_buckling_margin"] = (sigma_max_lower - inputs["lower_skin_comp_stress"]) / sigma_max_lower + # apply 40% to shear stress based on Elham 2014, Effect of wing-box structure on the optimum wing outer shape, The Aeronautical Journal + outputs["front_spar_buckling_margin"] = (tau_max_front - inputs["front_spar_shear_stress"] * 0.4) / tau_max_front + outputs["rear_spar_buckling_margin"] = (tau_max_rear - inputs["rear_spar_shear_stress"] * 0.4) / tau_max_rear + + # Failure criteria (must be <= 0) + failure_all = np.concatenate( + [ + inputs["upper_skin_comp_stress"] / sigma_max_upper - 1, + inputs["lower_skin_comp_stress"] / sigma_max_lower - 1, + inputs["front_spar_shear_stress"] * 0.4 / tau_max_front - 1, + inputs["rear_spar_shear_stress"] * 0.4 / tau_max_rear - 1, + ] + ) + + fmax = np.max(failure_all) + + nlog, nsum, nexp = np.log, np.sum, np.exp + ks = 1 / rho * nlog(nsum(nexp(rho * (failure_all - fmax)))) + outputs["failure_buckling"] = fmax + ks diff --git a/openaerostruct/structures/spatial_beam_functionals.py b/openaerostruct/structures/spatial_beam_functionals.py index cff8cdd76..447e97cb9 100644 --- a/openaerostruct/structures/spatial_beam_functionals.py +++ b/openaerostruct/structures/spatial_beam_functionals.py @@ -8,6 +8,7 @@ from openaerostruct.structures.non_intersecting_thickness import NonIntersectingThickness from openaerostruct.structures.failure_exact import FailureExact from openaerostruct.structures.failure_ks import FailureKS +from openaerostruct.structures.failure_buckling_ks import FailureBucklingKS class SpatialBeamFunctionals(om.Group): @@ -74,3 +75,16 @@ def setup(self): self.add_subsystem( "failure", FailureKS(surface=surface), promotes_inputs=["vonmises"], promotes_outputs=["failure"] ) + + # compute buckling failure + if "buckling" in surface and surface["buckling"]: + self.add_subsystem( + "failure_buckling", + FailureBucklingKS(surface=surface), + promotes_inputs=["skin_thickness", "spar_thickness", "t_over_c", "fem_chords"], + promotes_outputs=["failure_buckling"] + ) + self.connect("vonmises.upper_skin_comp_stress", "failure_buckling.upper_skin_comp_stress") + self.connect("vonmises.lower_skin_comp_stress", "failure_buckling.lower_skin_comp_stress") + self.connect("vonmises.front_spar_shear_stress", "failure_buckling.front_spar_shear_stress") + self.connect("vonmises.rear_spar_shear_stress", "failure_buckling.rear_spar_shear_stress") diff --git a/openaerostruct/structures/vonmises_wingbox.py b/openaerostruct/structures/vonmises_wingbox.py index 052bfbfef..115e6e404 100644 --- a/openaerostruct/structures/vonmises_wingbox.py +++ b/openaerostruct/structures/vonmises_wingbox.py @@ -42,6 +42,18 @@ class VonMisesWingbox(om.ExplicitComponent): ------- vonmises[ny-1, 4] : numpy array von Mises stresses for 4 stress combinations for each FEM element. + upper_skin_comp_stress[ny-1] : numpy array + Optional output for buckling failure check. + Compression stress in the upper skin for each FEM element. + lower_skin_comp_stress[ny-1] : numpy array + Optional output for buckling failure check. + Compression stress in the lower skin for each FEM element. + front_spar_shear_stress[ny-1] : numpy array + Optional output for buckling failure check. + Shear stress in the front spar for each FEM element. + rear_spar_shear_stress[ny-1] : numpy array + Optional output for buckling failure check. + Shear stress in the rear spar for each FEM element. """ @@ -66,6 +78,14 @@ def setup(self): self.add_output("vonmises", val=np.zeros((self.ny - 1, 4)), units="N/m**2") + # also output upper skin's compression stress and spar shear stress for buckling + if "buckling" in self.surface and self.surface["buckling"]: + # positive = compression + self.add_output("upper_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_output("lower_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_output("front_spar_shear_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.add_output("rear_spar_shear_stress", val=np.zeros(self.ny - 1), units="N/m**2") + self.E = surface["E"] self.G = surface["G"] @@ -164,3 +184,11 @@ def compute(self, inputs, outputs): np.sqrt((rear_bending_stress + axial_stress) ** 2 + 3 * (torsion_stress + vertical_shear) ** 2) / self.tssf ) + + # compute upper skin's compression stress and spar shear stress for buckling failure check + if "buckling" in self.surface and self.surface["buckling"]: + # flip sign to let positive = compression + outputs["upper_skin_comp_stress"][ielem] = -(axial_stress + top_bending_stress) + outputs["lower_skin_comp_stress"][ielem] = -(axial_stress + bottom_bending_stress) + outputs["front_spar_shear_stress"][ielem] = torsion_stress - vertical_shear + outputs["rear_spar_shear_stress"][ielem] = torsion_stress + vertical_shear From 36fddefbbc35802bbdfb01e11eaa0897a73635f8 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Thu, 3 Jul 2025 21:27:59 +0900 Subject: [PATCH 05/30] added safety factor for buckling stress --- openaerostruct/structures/failure_buckling_ks.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py index d723ea739..968196129 100644 --- a/openaerostruct/structures/failure_buckling_ks.py +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -81,6 +81,8 @@ def compute(self, inputs, outputs): rib_pitch = self.options["surface"]["rib_pitch"] rho = self.options["rho"] # KS aggregation parameter + safety_factor = 1.5 + # --- compute buckling stress --- # buckling model parameters eff_factor = 0.88 # stiffened panel efficiency for skin compression buckling. From Fig. 14.4.4 in Niu 1997 @@ -90,10 +92,10 @@ def compute(self, inputs, outputs): # Positive stress = compression. load_intensity_upper = inputs["skin_thickness"] * inputs["upper_skin_comp_stress"] # Pa*m load_positive_upper = np.maximum(load_intensity_upper, 1.) # NOTE: this is C1 discontinuous so may negatively affect convergence - sigma_max_upper = eff_factor * (load_positive_upper * E / rib_pitch)**0.5 + sigma_max_upper = (eff_factor * (load_positive_upper * E / rib_pitch)**0.5) / safety_factor load_intensity_lower = inputs["skin_thickness"] * inputs["lower_skin_comp_stress"] load_positive_lower = np.maximum(load_intensity_lower, 1.) - sigma_max_lower = eff_factor * (load_positive_lower * E / rib_pitch)**0.5 + sigma_max_lower = (eff_factor * (load_positive_lower * E / rib_pitch)**0.5) / safety_factor # shear buckling stress of spars. Eq. 14.4.3, Page 618 from Niu 1997 t_over_c_orig = surface["original_wingbox_airfoil_t_over_c"] @@ -101,8 +103,8 @@ def compute(self, inputs, outputs): spar_height_rear_nondim = (surface["data_y_upper"][-1] - surface["data_y_lower"][-1]) spar_height_front = spar_height_front_nondim * inputs["fem_chords"] * (inputs["t_over_c"] / t_over_c_orig) spar_height_rear = spar_height_rear_nondim * inputs["fem_chords"] * (inputs["t_over_c"] / t_over_c_orig) - tau_max_front = Ks * E * (inputs["spar_thickness"] / spar_height_front)**2 - tau_max_rear = Ks * E * (inputs["spar_thickness"] / spar_height_rear)**2 + tau_max_front = (Ks * E * (inputs["spar_thickness"] / spar_height_front)**2) / safety_factor + tau_max_rear = (Ks * E * (inputs["spar_thickness"] / spar_height_rear)**2) / safety_factor # --- compute failure and stress margins --- outputs["upper_skin_buckling_margin"] = (sigma_max_upper - inputs["upper_skin_comp_stress"]) / sigma_max_upper From fc02b9063b2932fab7d6bc82d48424f748b15cf7 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 16 Jul 2025 17:17:57 -0400 Subject: [PATCH 06/30] added Euler column buckling for strut --- .../integration/aerostruct_groups.py | 6 +- .../structures/failure_buckling_ks.py | 94 ++++++++++++++++++- .../structures/spatial_beam_functionals.py | 26 +++-- 3 files changed, 112 insertions(+), 14 deletions(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 37d7c8cdc..a9bf12f43 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -187,7 +187,11 @@ def setup(self): promotes_outputs = ["vonmises", "failure"] if "buckling" in surface and surface["buckling"]: promotes_inputs += ["skin_thickness", "t_over_c", "fem_chords"] - promotes_outputs += ["failure_buckling"] + promotes_outputs += ["failure_local_buckling"] + if surface["name"] == "strut": + # additional inputs/outputs for column buckling + promotes_inputs += ["Iz", "joint_load"] + promotes_outputs += ["failure_column_buckling"] self.add_subsystem( "struct_funcs", diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py index 968196129..ef68e0c92 100644 --- a/openaerostruct/structures/failure_buckling_ks.py +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -2,9 +2,9 @@ import openmdao.api as om -class FailureBucklingKS(om.ExplicitComponent): +class PanelLocalBucklingFailureKS(om.ExplicitComponent): """ - Buckling stress constraints for skin compression and spar shear buckling. + Panel local buckling stress constraints for skin compression and spar shear buckling. Returns KS-aggregated failure metric that should be constrained to be <= 0. Source: Michael Niu, Airframe Stress Analysis and Sizing, Hong Kong Conmilit Press Ltd., 1997 @@ -63,7 +63,7 @@ def setup(self): self.add_input("t_over_c", val=np.ones((self.ny - 1))) self.add_input("fem_chords", val=np.zeros(self.ny - 1), units="m") - self.add_output("failure_buckling", val=0.0) + self.add_output("failure_local_buckling", val=0.0, desc="Local buckling failure metric, should be <= 0") # also output margins for each element/buckling for post-processing # but don't declare partials for these outputs because we only use aggregated failure self.add_output("upper_skin_buckling_margin", val=np.zeros(self.ny - 1), desc='(buckling_limit - stress) / buckling_limit >= 0?') @@ -73,7 +73,7 @@ def setup(self): self.rho = rho - self.declare_partials("failure_buckling", "*", method="cs") + self.declare_partials("failure_local_buckling", "*", method="cs") def compute(self, inputs, outputs): surface = self.options["surface"] @@ -127,4 +127,88 @@ def compute(self, inputs, outputs): nlog, nsum, nexp = np.log, np.sum, np.exp ks = 1 / rho * nlog(nsum(nexp(rho * (failure_all - fmax)))) - outputs["failure_buckling"] = fmax + ks + outputs["failure_local_buckling"] = fmax + ks + + +class EulerColumnBucklingFailureKS(om.ExplicitComponent): + """ + Euler column buckling. + + Given the EI (which varies along the strut span) and strut length, it computes the buckling critical load. + The critical load varies along the strut span. + We them impose: + critical_load >= column_compression_load (which comes from the joint reaction force) + This is converted to a failure metric: + failure = (column_compression_load - critical_load) / column_compression_load <= 0. + and will be KS-aggregated. + + Parameters + ---------- + nodes[ny, 3] : numpy array + FEM node coordinates. + joint_load[3] : numpy array + Force vector applied to the strut at the joint. + Iz[ny-1] : numpy array + Moment of inertia of the strut about the z-axis. + + Returns + ------- + failure_Euler_buckling : float + Euler column buckling failure metric, should be <= 0 + Euler_buckling_margin[ny-1] : numpy array + """ + + def initialize(self): + self.options.declare("surface", types=dict) + self.options.declare("rho", types=float, default=100.0, desc="KS aggregation smoothness parameter") + + def setup(self): + surface = self.options["surface"] + ny = surface["mesh"].shape[1] + + self.add_input("nodes", shape=(ny, 3), units="m", desc="FEM node coordinates") + self.add_input("joint_load", shape=3, desc="load vector applied to the strut at the joint") + self.add_input("Iz", shape=(ny - 1), units="m**4") # use Iz because Iz < Iy + + self.add_output("failure_column_buckling", val=0.0, desc="Euler column buckling failure metric, should be <= 0") + self.add_output("column_buckling_margin", shape=(ny - 1), desc="Euler column buckling margin, should be >= 0") + # this buckling margin is only for post-processing, so don't declare partials for it + + self.declare_partials("failure_column_buckling", ["joint_load", "Iz"], method="cs") + # we only use nodes[0, :] and nodes[-1, :], so declare the sparse partials + col_indices = [0, 1, 2, (ny * 3) - 3, (ny * 3) - 2, (ny * 3) - 1] + self.declare_partials("column_buckling_margin", "nodes", method="cs", rows=np.zeros(6), cols=col_indices) + + def compute(self, inputs, outputs): + nodes = inputs["nodes"] + E = self.options["surface"]["E"] + joint_load = inputs["joint_load"] + Iz = inputs["Iz"] + + # convert joint load to N (FEM normalizes this by 1e9, so we apply it back) + joint_load *= 1e9 + + # strut length + # The length we compute from the FEM nodes are longer than the actual length because e.g. we don't model fuselage so the strut root is at the symmetry plane. + # To account for this, apply a factor < 1 to the strut length. + column_length_factor = self.options["surface"]["column_length_factor"] + strut_len = np.linalg.norm(nodes[0, :] - nodes[-1, :]) * column_length_factor + + # compression load + strut_direction = (nodes[-1, :] - nodes[0, :]) + strut_direction /= np.linalg.norm(strut_direction) # unit magniture + comp_load = np.dot(joint_load, strut_direction) # flip sign to make compression = positive + + # critical load with safety factor of 1.5 + sf = 1.5 + P_crit = np.pi**2 * E * Iz / strut_len**2 / sf + + # failure metric + outputs["column_buckling_margin"] = (P_crit - comp_load) / P_crit + + # KS-aggregated failure metric (should be <= 0) + failure_all = comp_load / P_crit - 1 + fmax = np.max(failure_all) + rho = self.options["rho"] + ks = 1 / rho * np.log(np.sum(np.exp(rho * (failure_all - fmax)))) + outputs["failure_column_buckling"] = fmax + ks \ No newline at end of file diff --git a/openaerostruct/structures/spatial_beam_functionals.py b/openaerostruct/structures/spatial_beam_functionals.py index 447e97cb9..212b88423 100644 --- a/openaerostruct/structures/spatial_beam_functionals.py +++ b/openaerostruct/structures/spatial_beam_functionals.py @@ -8,7 +8,7 @@ from openaerostruct.structures.non_intersecting_thickness import NonIntersectingThickness from openaerostruct.structures.failure_exact import FailureExact from openaerostruct.structures.failure_ks import FailureKS -from openaerostruct.structures.failure_buckling_ks import FailureBucklingKS +from openaerostruct.structures.failure_buckling_ks import PanelLocalBucklingFailureKS, EulerColumnBucklingFailureKS class SpatialBeamFunctionals(om.Group): @@ -78,13 +78,23 @@ def setup(self): # compute buckling failure if "buckling" in surface and surface["buckling"]: + # skin panel buckling and spar shear buckling self.add_subsystem( - "failure_buckling", - FailureBucklingKS(surface=surface), + "local_buckling", + PanelLocalBucklingFailureKS(surface=surface), promotes_inputs=["skin_thickness", "spar_thickness", "t_over_c", "fem_chords"], - promotes_outputs=["failure_buckling"] + promotes_outputs=["failure_local_buckling"] ) - self.connect("vonmises.upper_skin_comp_stress", "failure_buckling.upper_skin_comp_stress") - self.connect("vonmises.lower_skin_comp_stress", "failure_buckling.lower_skin_comp_stress") - self.connect("vonmises.front_spar_shear_stress", "failure_buckling.front_spar_shear_stress") - self.connect("vonmises.rear_spar_shear_stress", "failure_buckling.rear_spar_shear_stress") + self.connect("vonmises.upper_skin_comp_stress", "local_buckling.upper_skin_comp_stress") + self.connect("vonmises.lower_skin_comp_stress", "local_buckling.lower_skin_comp_stress") + self.connect("vonmises.front_spar_shear_stress", "local_buckling.front_spar_shear_stress") + self.connect("vonmises.rear_spar_shear_stress", "local_buckling.rear_spar_shear_stress") + + # global Euler column buckling + if surface["name"] == "strut": + self.add_subsystem( + "column_buckling", + EulerColumnBucklingFailureKS(surface=surface), + promotes_inputs=["nodes", "joint_load", "Iz"], + promotes_outputs=["failure_column_buckling"] + ) From 230fb3601d8bb2723ef8d0123afd5bb481026763 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 16 Jul 2025 18:49:45 -0400 Subject: [PATCH 07/30] support pin boundary condition for strut root --- openaerostruct/structures/create_rhs.py | 10 ++++++++- openaerostruct/structures/disp.py | 12 +++++++++-- openaerostruct/structures/fem.py | 9 +++++--- openaerostruct/structures/fem_strut_braced.py | 21 ++++++++++++++++--- 4 files changed, 43 insertions(+), 9 deletions(-) diff --git a/openaerostruct/structures/create_rhs.py b/openaerostruct/structures/create_rhs.py index 1992e27b2..870f3cd13 100644 --- a/openaerostruct/structures/create_rhs.py +++ b/openaerostruct/structures/create_rhs.py @@ -31,7 +31,15 @@ def setup(self): self.ny = surface["mesh"].shape[1] self.add_input("total_loads", val=np.zeros((self.ny, 6)), units="N") - self.add_output("forces", val=np.ones(((self.ny + 1) * 6)), units="N") + + # shape of forces depends on the root boundary condition type + if "root_BC_type" in surface and surface["root_BC_type"] == "pin": + self.root_BC_pin = True + forces_size = self.ny * 6 + 3 + else: + self.root_BC_pin = False + forces_size = (self.ny + 1) * 6 + self.add_output("forces", val=np.ones((forces_size)), units="N") n = self.ny * 6 arange = np.arange((n)) diff --git a/openaerostruct/structures/disp.py b/openaerostruct/structures/disp.py index 8d5f1e85f..ce3f687f4 100644 --- a/openaerostruct/structures/disp.py +++ b/openaerostruct/structures/disp.py @@ -35,7 +35,15 @@ def setup(self): self.ny = surface["mesh"].shape[1] - self.add_input("disp_aug", val=np.zeros(((self.ny + 1) * 6)), units="m") + # shape of disp_aug depends on the root boundary condition type + if "root_BC_type" in surface and surface["root_BC_type"] == "pin": + self.root_BC_pin = True + disp_aug_size = self.ny * 6 + 3 + else: + self.root_BC_pin = False + disp_aug_size = (self.ny + 1) * 6 + self.add_input("disp_aug", val=np.zeros((disp_aug_size)), units="m") + self.add_output("disp", val=np.zeros((self.ny, 6)), units="m") n = self.ny * 6 @@ -45,4 +53,4 @@ def setup(self): def compute(self, inputs, outputs): # Obtain the relevant portions of disp_aug and store the reshaped # displacements in disp - outputs["disp"] = inputs["disp_aug"][:-6].reshape((-1, 6)) + outputs["disp"] = inputs["disp_aug"][:self.ny * 6].reshape((-1, 6)) diff --git a/openaerostruct/structures/fem.py b/openaerostruct/structures/fem.py index 4c6492875..bb7c42fbd 100644 --- a/openaerostruct/structures/fem.py +++ b/openaerostruct/structures/fem.py @@ -7,7 +7,7 @@ import openmdao.api as om -def get_drdu_sparsity_pattern(ny, vec_size, symmetry): +def get_drdu_sparsity_pattern(ny, vec_size, symmetry, root_BC_pin=False): """ return sparsity pattern of partials of residuals w.r.t. displacements = stiffness matrix """ @@ -45,9 +45,12 @@ def get_drdu_sparsity_pattern(ny, vec_size, symmetry): index = 6 * idx num_dofs = 6 * ny - arange = np.arange(6) - # Fixed boundary condition. + # Fixed boundary condition. Number of additional "DoFs" varies depending on the root boundary condition type. + if root_BC_pin: + arange = np.arange(3) # translation only + else: + arange = np.arange(6) # translation and rotation rows6 = index + arange cols6 = num_dofs + arange diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 6487a3610..9a954f454 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -90,7 +90,14 @@ def setup(self): # --- inputs, outputs, partials of FEM system for each surface --- ny = surface["mesh"].shape[1] - size = int(6 * ny + 6) + if "root_BC_type" in surface and surface["root_BC_type"] == "pin": + # pin boundary condition at the root + size = int(6 * ny + 3) + root_BC_pin = True + else: + # rigid boundary condition at the root + size = int(6 * ny + 6) + root_BC_pin = False self.ny.append(ny) self.size.append(size) @@ -107,7 +114,7 @@ def setup(self): # The derivative of residual wrt displacements is the stiffness matrix K. We can use the # sparsity pattern here and when constucting the sparse matrix, so save rows and cols. - rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"]) + rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"], root_BC_pin) if i == 1: # the second surface (strut) goes to the lower-right block of the entire stiffness matrix rows += self.size[0] @@ -334,7 +341,15 @@ def assemble_CSC_K(self, inputs): data3 = k_loc[0, :6, :6].flatten() data4 = k_loc[-1, 6:, 6:].flatten() data5 = (k_loc[0:-1, 6:, 6:] + k_loc[1:, :6, :6]).flatten() - data6 = np.full((6,), 1e9) + + # data6 corresponds to the root boundary condition + if "root_BC_type" in surface and surface["root_BC_type"] == "pin": + # for pin BC (constraint translation only) + data6 = np.full((3,), 1e9) + else: + # for rigid BC (constraint translation and rotation) + data6 = np.full((6,), 1e9) + data[name] = np.concatenate([data1, data2, data3, data4, data5, data6, data6]) # combine matrices for wing, strut, and joint constraints From 96962df57f0cdf0126a2e8f3318c175f77875265 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 16 Jul 2025 19:18:26 -0400 Subject: [PATCH 08/30] pin and ball boundary conditions for strut root --- openaerostruct/structures/create_rhs.py | 12 ++++---- openaerostruct/structures/disp.py | 12 ++++---- openaerostruct/structures/fem.py | 20 +++++++++---- openaerostruct/structures/fem_strut_braced.py | 29 ++++++++++++------- 4 files changed, 44 insertions(+), 29 deletions(-) diff --git a/openaerostruct/structures/create_rhs.py b/openaerostruct/structures/create_rhs.py index 870f3cd13..20bc56492 100644 --- a/openaerostruct/structures/create_rhs.py +++ b/openaerostruct/structures/create_rhs.py @@ -33,13 +33,13 @@ def setup(self): self.add_input("total_loads", val=np.zeros((self.ny, 6)), units="N") # shape of forces depends on the root boundary condition type - if "root_BC_type" in surface and surface["root_BC_type"] == "pin": - self.root_BC_pin = True - forces_size = self.ny * 6 + 3 + if "root_BC_type" in surface and surface["root_BC_type"] == "ball": + dof_of_boundary = 3 # translation only + elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": + dof_of_boundary = 5 # translation and rotation in y and z else: - self.root_BC_pin = False - forces_size = (self.ny + 1) * 6 - self.add_output("forces", val=np.ones((forces_size)), units="N") + dof_of_boundary = 6 # translation and rotation + self.add_output("forces", val=np.ones((self.ny * 6 + dof_of_boundary)), units="N") n = self.ny * 6 arange = np.arange((n)) diff --git a/openaerostruct/structures/disp.py b/openaerostruct/structures/disp.py index ce3f687f4..44bc0f27f 100644 --- a/openaerostruct/structures/disp.py +++ b/openaerostruct/structures/disp.py @@ -36,13 +36,13 @@ def setup(self): self.ny = surface["mesh"].shape[1] # shape of disp_aug depends on the root boundary condition type - if "root_BC_type" in surface and surface["root_BC_type"] == "pin": - self.root_BC_pin = True - disp_aug_size = self.ny * 6 + 3 + if "root_BC_type" in surface and surface["root_BC_type"] == "ball": + dof_of_boundary = 3 # translation only + elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": + dof_of_boundary = 5 # translation and rotation in y and z else: - self.root_BC_pin = False - disp_aug_size = (self.ny + 1) * 6 - self.add_input("disp_aug", val=np.zeros((disp_aug_size)), units="m") + dof_of_boundary = 6 # translation and rotation + self.add_input("disp_aug", val=np.zeros((self.ny * 6 + dof_of_boundary)), units="m") self.add_output("disp", val=np.zeros((self.ny, 6)), units="m") diff --git a/openaerostruct/structures/fem.py b/openaerostruct/structures/fem.py index bb7c42fbd..2f3323194 100644 --- a/openaerostruct/structures/fem.py +++ b/openaerostruct/structures/fem.py @@ -7,7 +7,7 @@ import openmdao.api as om -def get_drdu_sparsity_pattern(ny, vec_size, symmetry, root_BC_pin=False): +def get_drdu_sparsity_pattern(ny, vec_size, symmetry, root_BC="rigid"): """ return sparsity pattern of partials of residuals w.r.t. displacements = stiffness matrix """ @@ -47,12 +47,20 @@ def get_drdu_sparsity_pattern(ny, vec_size, symmetry, root_BC_pin=False): num_dofs = 6 * ny # Fixed boundary condition. Number of additional "DoFs" varies depending on the root boundary condition type. - if root_BC_pin: - arange = np.arange(3) # translation only + if root_BC == "rigid": # translation and rotation + row_arange = np.arange(6) + col_arange = np.arange(6) + elif root_BC == "pin": # translation and rotation in y and z + row_arange = np.array([0, 1, 2, 4, 5]) # exclude rotation in x + col_arange = np.arange(5) + elif root_BC == "ball": # translation only + row_arange = np.arange(3) + col_arange = np.arange(3) else: - arange = np.arange(6) # translation and rotation - rows6 = index + arange - cols6 = num_dofs + arange + raise ValueError(f"Invalid root boundary condition type: {root_BC}") + + rows6 = index + row_arange + cols6 = num_dofs + col_arange rows = np.concatenate([rows1, rows2, rows3, rows4, rows5, rows6, cols6]) cols = np.concatenate([cols1, cols2, cols3, cols4, cols5, cols6, rows6]) diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 9a954f454..817d68ee1 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -10,7 +10,6 @@ # TODO: # - implement ball and pin joint (current pin joint should be renamed to ball) -# - implement pin/ball joint for the root symmetric boundary class FEMStrutBraced(om.ImplicitComponent): @@ -90,15 +89,20 @@ def setup(self): # --- inputs, outputs, partials of FEM system for each surface --- ny = surface["mesh"].shape[1] - if "root_BC_type" in surface and surface["root_BC_type"] == "pin": - # pin boundary condition at the root - size = int(6 * ny + 3) - root_BC_pin = True + if "root_BC_type" in surface and surface["root_BC_type"] == "ball": + # ball boundary condition at the root. Constrain translation only. + dof_of_boundary = 3 + root_BC = "ball" + elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": + # pin boundary condition at the root. Constrain translation and rotation in y and z. + dof_of_boundary = 5 + root_BC = "pin" else: - # rigid boundary condition at the root - size = int(6 * ny + 6) - root_BC_pin = False + # rigid boundary condition at the root. Constrain translation and rotation. + dof_of_boundary = 6 + root_BC = "rigid" self.ny.append(ny) + size = int(6 * ny + dof_of_boundary) self.size.append(size) full_size = size * vec_size @@ -114,7 +118,7 @@ def setup(self): # The derivative of residual wrt displacements is the stiffness matrix K. We can use the # sparsity pattern here and when constucting the sparse matrix, so save rows and cols. - rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"], root_BC_pin) + rows, cols, vec_rows, vec_cols = get_drdu_sparsity_pattern(ny, vec_size, surface["symmetry"], root_BC) if i == 1: # the second surface (strut) goes to the lower-right block of the entire stiffness matrix rows += self.size[0] @@ -343,9 +347,12 @@ def assemble_CSC_K(self, inputs): data5 = (k_loc[0:-1, 6:, 6:] + k_loc[1:, :6, :6]).flatten() # data6 corresponds to the root boundary condition - if "root_BC_type" in surface and surface["root_BC_type"] == "pin": - # for pin BC (constraint translation only) + if "root_BC_type" in surface and surface["root_BC_type"] == "ball": + # for ball BC (constraint translation only) data6 = np.full((3,), 1e9) + elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": + # for pin BC (constraint translation and rotation in y and z) + data6 = np.full((5,), 1e9) else: # for rigid BC (constraint translation and rotation) data6 = np.full((6,), 1e9) From e400ce5f396bcb0e275cb1d298fe194baa82847f Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 16 Jul 2025 19:35:58 -0400 Subject: [PATCH 09/30] added pin joint; renamed previous joint to ball joint --- openaerostruct/structures/fem_strut_braced.py | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 817d68ee1..2b152d23e 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -8,9 +8,6 @@ from .fem import get_drdu_sparsity_pattern, get_drdK_sparsity_pattern -# TODO: -# - implement ball and pin joint (current pin joint should be renamed to ball) - class FEMStrutBraced(om.ImplicitComponent): """ @@ -133,10 +130,12 @@ def setup(self): # --- joint constraints (coupling terms between two surfaces) --- joint_type = self.surfaces[0]["joint_type"] - if joint_type == 'pin': + if joint_type == 'ball': self.n_con = n_con = 3 # number of joint constraints: translational only + elif joint_type == 'pin': + self.n_con = n_con = 5 # translational and rotational in y and z elif joint_type == 'rigid': - self.n_con = n_con = 6 # number of joint constraints: translational and rotational + self.n_con = n_con = 6 # ranslational and rotational in x, y, and z else: raise ValueError("Joint type must be either 'pin' or 'rigid'.") @@ -163,7 +162,10 @@ def setup(self): # partials of joint constraint residuals w.r.t. displacement rows = np.arange(n_con) - cols = np.arange(n_con) + 6 * joint_idx + if joint_type in ['ball', 'rigid']: + cols = np.arange(n_con) + 6 * joint_idx + elif joint_type == 'pin': + cols = np.array([0, 1, 2, 4, 5]) + 6 * joint_idx # exclude x-rotation if i == 0: vals = np.ones(n_con) * 1e9 else: From 42d5a053d53b0b516ca7916a541c7f0f4ec28ec3 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 16 Jul 2025 21:22:50 -0400 Subject: [PATCH 10/30] fixed Euler buckling derivative --- openaerostruct/structures/failure_buckling_ks.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py index ef68e0c92..881dbdb83 100644 --- a/openaerostruct/structures/failure_buckling_ks.py +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -185,9 +185,6 @@ def compute(self, inputs, outputs): joint_load = inputs["joint_load"] Iz = inputs["Iz"] - # convert joint load to N (FEM normalizes this by 1e9, so we apply it back) - joint_load *= 1e9 - # strut length # The length we compute from the FEM nodes are longer than the actual length because e.g. we don't model fuselage so the strut root is at the symmetry plane. # To account for this, apply a factor < 1 to the strut length. @@ -197,7 +194,8 @@ def compute(self, inputs, outputs): # compression load strut_direction = (nodes[-1, :] - nodes[0, :]) strut_direction /= np.linalg.norm(strut_direction) # unit magniture - comp_load = np.dot(joint_load, strut_direction) # flip sign to make compression = positive + # multiply by 1e9 to make the units to N (FEM normalizes this by 1e9) + comp_load = np.dot(joint_load, strut_direction) * 1e9 # critical load with safety factor of 1.5 sf = 1.5 From 114193007e6c829cf8b904c668caa79f215c85b7 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Wed, 30 Jul 2025 17:10:11 -0400 Subject: [PATCH 11/30] refactored wing-strut joint definitions --- .../integration/aerostruct_groups.py | 11 +-- openaerostruct/structures/fem_strut_braced.py | 72 +++++++++++++------ .../structures/spatial_beam_functionals.py | 20 +++--- openaerostruct/structures/vonmises_wingbox.py | 4 +- 4 files changed, 70 insertions(+), 37 deletions(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index a9bf12f43..040f62228 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -185,13 +185,14 @@ def setup(self): elif surface["fem_model_type"] == "wingbox": promotes_inputs = ["Qz", "J", "A_enc", "spar_thickness", "htop", "hbottom", "hfront", "hrear", "nodes", "disp"] promotes_outputs = ["vonmises", "failure"] - if "buckling" in surface and surface["buckling"]: + + # additional inputs/outputs for buckling + if "panel_buckling" in surface and surface["panel_buckling"]: promotes_inputs += ["skin_thickness", "t_over_c", "fem_chords"] promotes_outputs += ["failure_local_buckling"] - if surface["name"] == "strut": - # additional inputs/outputs for column buckling - promotes_inputs += ["Iz", "joint_load"] - promotes_outputs += ["failure_column_buckling"] + if "column_buckling" in surface and surface["column_buckling"]: + promotes_inputs += ["Iz", "joint_load"] + promotes_outputs += ["failure_column_buckling"] self.add_subsystem( "struct_funcs", diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 2b152d23e..97e1605c5 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -68,11 +68,52 @@ def setup(self): """ Matrix and RHS are inputs, solution vector is the output. """ + # --- input option check and processings --- vec_size = self.options["vec_size"] self.surfaces = self.options["surfaces"] - if len(self.surfaces) != 2: - raise ValueError("FEMStrutBraced component requires exactly two surfaces (wing and strut).") - self.surface_names = [self.surfaces[0]["name"], self.surfaces[1]["name"]] + self.surface_names = [surface["name"] for surface in self.surfaces] + + # surfaces must be either [wing, strut] or [wing, strut, jury] in this order + if len(self.surfaces) == 2: + if self.surface_names != ["wing", "strut"]: + raise ValueError("FEMStrutBraced component requires the surfaces to be in the order of [wing, strut].") + include_jury = False + self.wing_surface = self.surfaces[0] + self.strut_surface = self.surfaces[1] + elif len(self.surfaces) == 3: + if self.surface_names != ["wing", "strut", "jury"]: + raise ValueError("FEMStrutBraced component requires the surfaces to be in the order of [wing, strut, jury].") + include_jury = True + self.wing_surface = self.surfaces[0] + self.strut_surface = self.surfaces[1] + self.jury_surface = self.surfaces[2] + else: + raise ValueError("FEMStrutBraced component requires exactly two or three surfaces (wing, strut, and/or jury).") + + # no support for asymmetric surfaces + for surface in self.surfaces: + if not surface["symmetry"]: + raise NotImplementedError("Strut-braced wing joint is only implemented for symmetric surfaces.") + + # We only support the surface name to be "wing" and "strut", and "jury" for now + for surface in self.surfaces: + if surface["name"] not in ["wing", "strut", "jury"]: + raise ValueError("FEMStrutBraced component does not support surface name: " + surface["name"]) + + # get joint information + wing_strut_joint_type = self.strut_surface["wing_strut_joint_type"] + wing_strut_joint_y = self.strut_surface["wing_strut_joint_y"] + # check sign convention for y (spanwise coordinate) + if wing_strut_joint_y * self.wing_surface["mesh"][0, 0, 1] < 0: + wing_strut_joint_y *= -1 + + if include_jury: + wing_jury_joint_type = self.jury_surface["wing_strut_joint_type"] + wing_jury_joint_y = self.jury_surface["wing_strut_joint_y"] + strut_jury_joint_type = self.jury_surface["strut_joint_type"] + strut_jury_joint_y = self.jury_surface["strut_joint_y"] + + # prepare inputs self.ny = [] # number of spanwise nodes self.size = [] # size of assembled K matrix for each surface self.num_k_data = [] # number of non-zero entries for assembled K matrix for each surface @@ -128,13 +169,12 @@ def setup(self): rows, cols = get_drdK_sparsity_pattern(ny) self.declare_partials(f"disp_aug_{name}", f"local_stiff_transformed_{name}", rows=rows, cols=cols) - # --- joint constraints (coupling terms between two surfaces) --- - joint_type = self.surfaces[0]["joint_type"] - if joint_type == 'ball': + # --- wing-strut joint constraints (coupling terms between two surfaces) --- + if wing_strut_joint_type == 'ball': self.n_con = n_con = 3 # number of joint constraints: translational only - elif joint_type == 'pin': + elif wing_strut_joint_type == 'pin': self.n_con = n_con = 5 # translational and rotational in y and z - elif joint_type == 'rigid': + elif wing_strut_joint_type == 'rigid': self.n_con = n_con = 6 # ranslational and rotational in x, y, and z else: raise ValueError("Joint type must be either 'pin' or 'rigid'.") @@ -145,26 +185,18 @@ def setup(self): k_rows_joint = [] k_cols_joint = [] self.k_data_joint = [] - for i, surface in enumerate(self.surfaces): + for i, surface in enumerate([self.wing_surface, self.strut_surface]): name = surface["name"] - # find the node index for the joint - joint_y = surface["joint_y"] - if not surface["symmetry"]: - raise NotImplementedError("Strut-braced wing joint is only implemented for symmetric surfaces.") - # check sign convention for y (spanwise coordinate) - if joint_y * surface["mesh"][0, 0, 1] < 0: - joint_y *= -1 - - joint_idx = np.argmin(np.abs(surface["mesh"][0, :, 1] - joint_y)) + joint_idx = np.argmin(np.abs(surface["mesh"][0, :, 1] - wing_strut_joint_y)) joint_indices.append(joint_idx) print("Surface", name, "joint y:", surface["mesh"][0, joint_idx, 1], "joint index:", joint_idx) # partials of joint constraint residuals w.r.t. displacement rows = np.arange(n_con) - if joint_type in ['ball', 'rigid']: + if wing_strut_joint_type in ['ball', 'rigid']: cols = np.arange(n_con) + 6 * joint_idx - elif joint_type == 'pin': + elif wing_strut_joint_type == 'pin': cols = np.array([0, 1, 2, 4, 5]) + 6 * joint_idx # exclude x-rotation if i == 0: vals = np.ones(n_con) * 1e9 diff --git a/openaerostruct/structures/spatial_beam_functionals.py b/openaerostruct/structures/spatial_beam_functionals.py index 212b88423..705ad5163 100644 --- a/openaerostruct/structures/spatial_beam_functionals.py +++ b/openaerostruct/structures/spatial_beam_functionals.py @@ -76,8 +76,8 @@ def setup(self): "failure", FailureKS(surface=surface), promotes_inputs=["vonmises"], promotes_outputs=["failure"] ) - # compute buckling failure - if "buckling" in surface and surface["buckling"]: + # compute panel local buckling failure + if "panel_buckling" in surface and surface["panel_buckling"]: # skin panel buckling and spar shear buckling self.add_subsystem( "local_buckling", @@ -90,11 +90,11 @@ def setup(self): self.connect("vonmises.front_spar_shear_stress", "local_buckling.front_spar_shear_stress") self.connect("vonmises.rear_spar_shear_stress", "local_buckling.rear_spar_shear_stress") - # global Euler column buckling - if surface["name"] == "strut": - self.add_subsystem( - "column_buckling", - EulerColumnBucklingFailureKS(surface=surface), - promotes_inputs=["nodes", "joint_load", "Iz"], - promotes_outputs=["failure_column_buckling"] - ) + # global Euler column buckling + if "column_buckling" in surface and surface["column_buckling"]: + self.add_subsystem( + "column_buckling", + EulerColumnBucklingFailureKS(surface=surface), + promotes_inputs=["nodes", "joint_load", "Iz"], + promotes_outputs=["failure_column_buckling"] + ) diff --git a/openaerostruct/structures/vonmises_wingbox.py b/openaerostruct/structures/vonmises_wingbox.py index 115e6e404..68b80dfac 100644 --- a/openaerostruct/structures/vonmises_wingbox.py +++ b/openaerostruct/structures/vonmises_wingbox.py @@ -79,7 +79,7 @@ def setup(self): self.add_output("vonmises", val=np.zeros((self.ny - 1, 4)), units="N/m**2") # also output upper skin's compression stress and spar shear stress for buckling - if "buckling" in self.surface and self.surface["buckling"]: + if "panel_buckling" in self.surface and self.surface["panel_buckling"]: # positive = compression self.add_output("upper_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") self.add_output("lower_skin_comp_stress", val=np.zeros(self.ny - 1), units="N/m**2") @@ -186,7 +186,7 @@ def compute(self, inputs, outputs): ) # compute upper skin's compression stress and spar shear stress for buckling failure check - if "buckling" in self.surface and self.surface["buckling"]: + if "panel_buckling" in self.surface and self.surface["panel_buckling"]: # flip sign to let positive = compression outputs["upper_skin_comp_stress"][ielem] = -(axial_stress + top_bending_stress) outputs["lower_skin_comp_stress"][ielem] = -(axial_stress + bottom_bending_stress) From b4f8ef089cec0051b3f70591fa5368c25cedb2ed Mon Sep 17 00:00:00 2001 From: kanekosh Date: Thu, 31 Jul 2025 20:12:11 -0400 Subject: [PATCH 12/30] working on jury strut, code runs but results seem wrong --- .../functionals/total_performance.py | 32 +- .../integration/aerostruct_groups.py | 60 +++- openaerostruct/structures/create_rhs.py | 2 + openaerostruct/structures/disp.py | 2 + openaerostruct/structures/fem.py | 3 + openaerostruct/structures/fem_strut_braced.py | 285 ++++++++++++++---- .../structures/spatial_beam_states.py | 13 + 7 files changed, 318 insertions(+), 79 deletions(-) diff --git a/openaerostruct/functionals/total_performance.py b/openaerostruct/functionals/total_performance.py index 4145a0cee..d14b0e7e8 100644 --- a/openaerostruct/functionals/total_performance.py +++ b/openaerostruct/functionals/total_performance.py @@ -17,13 +17,33 @@ def initialize(self): self.options.declare("surfaces", types=list) self.options.declare("user_specified_Sref", types=bool) self.options.declare("internally_connect_fuelburn", types=bool, default=True) + self.options.declare("strut_braced", default=False, types=bool) def setup(self): surfaces = self.options["surfaces"] + if self.options["strut_braced"]: + if len(surfaces) == 2: + # wing + strut. Make sure that the given surfaces are in order + if surfaces[0]["name"] == "wing" and surfaces[1]["name"] == "strut": + surfaces_all = surfaces + surfaces_AS = surfaces + pass + else: + raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + elif len(surfaces) == 3: + # wing + strut + jury + if surfaces[0]["name"] == "wing" and surfaces[1]["name"] == "strut" and surfaces[2]["name"] == "jury": + surfaces_all = surfaces + surfaces_AS = [surfaces[0], surfaces[1]] # exclude jury from VLM + else: + raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + else: + raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + if not self.options["user_specified_Sref"]: self.add_subsystem( - "sum_areas", SumAreas(surfaces=surfaces), promotes_inputs=["*S_ref"], promotes_outputs=["S_ref_total"] + "sum_areas", SumAreas(surfaces=surfaces_AS), promotes_inputs=["*S_ref"], promotes_outputs=["S_ref_total"] ) if self.options["internally_connect_fuelburn"]: @@ -33,21 +53,21 @@ def setup(self): self.add_subsystem( "CL_CD", - TotalLiftDrag(surfaces=surfaces), + TotalLiftDrag(surfaces=surfaces_AS), promotes_inputs=["*CL", "*CD", "*S_ref", "S_ref_total", "rho", "v"], promotes_outputs=["CL", "CD", "L", "D"], ) self.add_subsystem( "fuelburn", - BreguetRange(surfaces=surfaces), + BreguetRange(surfaces=surfaces_all), promotes_inputs=["*structural_mass", "CL", "CD", "CT", "speed_of_sound", "R", "Mach_number", "W0"], promotes_outputs=["fuelburn"], ) self.add_subsystem( "L_equals_W", - Equilibrium(surfaces=surfaces), + Equilibrium(surfaces=surfaces_all), promotes_inputs=["CL", "*structural_mass", "S_ref_total", "W0", "load_factor", "rho", "v"] + promote_fuelburn, promotes_outputs=["L_equals_W", "total_weight"], @@ -55,7 +75,7 @@ def setup(self): self.add_subsystem( "CG", - CenterOfGravity(surfaces=surfaces), + CenterOfGravity(surfaces=surfaces_all), promotes_inputs=["*structural_mass", "*cg_location", "total_weight", "W0", "empty_cg", "load_factor"] + promote_fuelburn, promotes_outputs=["cg"], @@ -63,7 +83,7 @@ def setup(self): self.add_subsystem( "moment", - MomentCoefficient(surfaces=surfaces), + MomentCoefficient(surfaces=surfaces_AS), promotes_inputs=["v", "rho", "cg", "S_ref_total", "*b_pts", "*widths", "*chords", "*sec_forces", "*S_ref"], promotes_outputs=["CM"], ) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 040f62228..8ecd2f1f6 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -227,17 +227,48 @@ def setup(self): coupled = om.Group() if self.options["strut_braced"]: - if len(surfaces) != 2: - raise NameError("Strut-braced wing requires exactly two surfaces.") - # add FEM directly under the coupled group because multiple surfaces are coupled in the FEM level + if len(surfaces) == 2: + # wing + strut. Make sure that the given surfaces are in order + if surfaces[0]["name"] == "wing" and surfaces[1]["name"] == "strut": + surfaces_all = surfaces + surfaces_AS = surfaces + pass + else: + raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + elif len(surfaces) == 3: + # wing + strut + jury + if surfaces[0]["name"] == "wing" and surfaces[1]["name"] == "strut" and surfaces[2]["name"] == "jury": + surfaces_all = surfaces + surfaces_AS = [surfaces[0], surfaces[1]] # exclude jury from VLM + pass + else: + raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + + # disable jury strut root boundary condition + surfaces[2]["root_BC_type"] = "none" + surfaces_all[2]["root_BC_type"] = "none" + + # add coupled FEM group coupled.add_subsystem( "FEM_SBW", - SpatialBeamStates(surface=surfaces, strut_braced=True), + SpatialBeamStates(surface=surfaces_all, strut_braced=True), promotes_inputs=["*"], promotes_outputs=["disp_*"] ) + else: + surfaces_all = surfaces + surfaces_AS = surfaces - for surface in surfaces: + # connections for all structural components + for surface in surfaces_all: + name = surface["name"] + + # connect displacements to the performance group + disp_source = f"coupled.disp_{name}" if self.options["strut_braced"] else f"coupled.{name}.disp" + self.connect(disp_source, name + "_perf.disp") + + # connections for all aerostructural surfaces + for surface in surfaces_AS: name = surface["name"] # Connect the output of the loads component with the FEM @@ -270,8 +301,6 @@ def setup(self): # Connect parameters from the 'coupled' group to the performance # groups for the individual surfaces. - disp_source = f"coupled.disp_{name}" if self.options["strut_braced"] else f"coupled.{name}.disp" - self.connect(disp_source, name + "_perf.disp") self.connect("coupled." + name + ".S_ref", name + "_perf.S_ref") self.connect("coupled." + name + ".widths", name + "_perf.widths") # self.connect('coupled.' + name + '.chords', name + '_perf.chords') @@ -308,10 +337,10 @@ def setup(self): ground_effect = True if self.options["compressible"] is True: - aero_states = CompressibleVLMStates(surfaces=surfaces, rotational=rotational) + aero_states = CompressibleVLMStates(surfaces=surfaces_AS, rotational=rotational) prom_in = ["v", "alpha", "beta", "rho", "Mach_number"] else: - aero_states = VLMStates(surfaces=surfaces, rotational=rotational) + aero_states = VLMStates(surfaces=surfaces_AS, rotational=rotational) prom_in = ["v", "alpha", "beta", "rho"] if ground_effect: prom_in.append("height_agl") @@ -322,7 +351,7 @@ def setup(self): # Explicitly connect parameters from each surface's group and the common # 'aero_states' group. - for surface in surfaces: + for surface in surfaces_AS: name = surface["name"] # Add a loads component to the coupled group @@ -363,7 +392,8 @@ def setup(self): # Add the coupled group to the model problem self.add_subsystem("coupled", coupled, promotes_inputs=prom_in) - for surface in surfaces: + # post-AS-analysis computations for aerostructural surfaces + for surface in surfaces_AS: name = surface["name"] # Add a performance group which evaluates the data after solving @@ -374,6 +404,13 @@ def setup(self): name + "_perf", perf_group, promotes_inputs=["rho", "v", "alpha", "beta", "re", "Mach_number"] ) + # post-AS-analysis computations for structure-only components + for surface in [s for s in surfaces_all if s not in surfaces_AS]: + self.add_subsystem( + f"{surface["name"]}_perf", + SpatialBeamFunctionals(surface=surface), + ) + # Add functionals to evaluate performance of the system. # Note that only the interesting results are promoted here; not all # of the parameters. @@ -381,6 +418,7 @@ def setup(self): "total_perf", TotalPerformance( surfaces=surfaces, + strut_braced=self.options["strut_braced"], user_specified_Sref=self.options["user_specified_Sref"], internally_connect_fuelburn=self.options["internally_connect_fuelburn"], ), diff --git a/openaerostruct/structures/create_rhs.py b/openaerostruct/structures/create_rhs.py index 20bc56492..52ced83eb 100644 --- a/openaerostruct/structures/create_rhs.py +++ b/openaerostruct/structures/create_rhs.py @@ -37,6 +37,8 @@ def setup(self): dof_of_boundary = 3 # translation only elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": dof_of_boundary = 5 # translation and rotation in y and z + elif "root_BC_type" in surface and surface["root_BC_type"] == "none": + dof_of_boundary = 0 # no boundary conditions (for jury strut) else: dof_of_boundary = 6 # translation and rotation self.add_output("forces", val=np.ones((self.ny * 6 + dof_of_boundary)), units="N") diff --git a/openaerostruct/structures/disp.py b/openaerostruct/structures/disp.py index 44bc0f27f..60f1052ec 100644 --- a/openaerostruct/structures/disp.py +++ b/openaerostruct/structures/disp.py @@ -40,6 +40,8 @@ def setup(self): dof_of_boundary = 3 # translation only elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": dof_of_boundary = 5 # translation and rotation in y and z + elif "root_BC_type" in surface and surface["root_BC_type"] == "none": + dof_of_boundary = 0 # no boundary conditions (for jury strut) else: dof_of_boundary = 6 # translation and rotation self.add_input("disp_aug", val=np.zeros((self.ny * 6 + dof_of_boundary)), units="m") diff --git a/openaerostruct/structures/fem.py b/openaerostruct/structures/fem.py index 2f3323194..74d884ba8 100644 --- a/openaerostruct/structures/fem.py +++ b/openaerostruct/structures/fem.py @@ -56,6 +56,9 @@ def get_drdu_sparsity_pattern(ny, vec_size, symmetry, root_BC="rigid"): elif root_BC == "ball": # translation only row_arange = np.arange(3) col_arange = np.arange(3) + elif root_BC == "none": # no root BC (e.g. for jury strut) + row_arange = np.array([]) + col_arange = np.array([]) else: raise ValueError(f"Invalid root boundary condition type: {root_BC}") diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 97e1605c5..628cdf698 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -9,6 +9,46 @@ from .fem import get_drdu_sparsity_pattern, get_drdK_sparsity_pattern +def _get_joint_size(joint_type): + """ + Return the number of joint constraints for a given joint type. + """ + if joint_type == 'ball': + n_con = 3 # number of joint constraints: translational only + elif joint_type == 'pin': + n_con = 5 # translational and rotational in y and z + elif joint_type == 'rigid': + n_con = 6 # ranslational and rotational in x, y, and z + else: + raise ValueError("Joint type must be either 'pin' or 'rigid'.") + return n_con + +def _get_joint_idx_from_y(surface, joint_y): + """ + Return the joint index from the joint y (spanwise) coordinate. + This is only good for wing or strut, but do not use for jury (or any vertical surface) + """ + if surface["name"] == "jury": + raise RuntimeError("Cannot use _get_joint_idx_from_y for jury surface.") + joint_idx = np.argmin((np.abs(surface["mesh"][0, :, 1]) - abs(joint_y))**2) + return joint_idx + +def _joint_rows_cols(joint_type, surface, joint_idx): + """ + Return rows and cols indices for joint constraints. + rows and cols are defined for the component-level indexing, but not the global wing-strut-jury system. + """ + # partials of joint constraint residuals w.r.t. displacement + n_con = _get_joint_size(joint_type) + rows = np.arange(n_con) + if joint_type in ['ball', 'rigid']: + cols = np.arange(n_con) + 6 * joint_idx + elif joint_type == 'pin': + cols = np.array([0, 1, 2, 4, 5]) + 6 * joint_idx # exclude x-rotation + + return rows, cols + + class FEMStrutBraced(om.ImplicitComponent): """ Component that solves an FEM linear system, K u = f. @@ -77,13 +117,13 @@ def setup(self): if len(self.surfaces) == 2: if self.surface_names != ["wing", "strut"]: raise ValueError("FEMStrutBraced component requires the surfaces to be in the order of [wing, strut].") - include_jury = False + self.include_jury = False self.wing_surface = self.surfaces[0] self.strut_surface = self.surfaces[1] elif len(self.surfaces) == 3: if self.surface_names != ["wing", "strut", "jury"]: raise ValueError("FEMStrutBraced component requires the surfaces to be in the order of [wing, strut, jury].") - include_jury = True + self.include_jury = True self.wing_surface = self.surfaces[0] self.strut_surface = self.surfaces[1] self.jury_surface = self.surfaces[2] @@ -103,24 +143,20 @@ def setup(self): # get joint information wing_strut_joint_type = self.strut_surface["wing_strut_joint_type"] wing_strut_joint_y = self.strut_surface["wing_strut_joint_y"] - # check sign convention for y (spanwise coordinate) - if wing_strut_joint_y * self.wing_surface["mesh"][0, 0, 1] < 0: - wing_strut_joint_y *= -1 - if include_jury: - wing_jury_joint_type = self.jury_surface["wing_strut_joint_type"] - wing_jury_joint_y = self.jury_surface["wing_strut_joint_y"] - strut_jury_joint_type = self.jury_surface["strut_joint_type"] - strut_jury_joint_y = self.jury_surface["strut_joint_y"] + if self.include_jury: + wing_jury_joint_type = self.jury_surface["wing_jury_joint_type"] + wing_jury_joint_y = self.jury_surface["wing_jury_joint_y"] + strut_jury_joint_type = self.jury_surface["strut_jury_joint_type"] + strut_jury_joint_y = self.jury_surface["strut_jury_joint_y"] # prepare inputs self.ny = [] # number of spanwise nodes - self.size = [] # size of assembled K matrix for each surface + self.size = [] # list of size (num rows) of assembled K matrix for each surface self.num_k_data = [] # number of non-zero entries for assembled K matrix for each surface self._lup = [] k_cols = [] # column indices for stiffness matrix sparsity pattern k_rows = [] # row indices for stiffness matrix sparsity pattern - joint_indices = [] # node index for the joint between wing and strut for i, surface in enumerate(self.surfaces): name = surface["name"] @@ -135,6 +171,10 @@ def setup(self): # pin boundary condition at the root. Constrain translation and rotation in y and z. dof_of_boundary = 5 root_BC = "pin" + elif "root_BC_type" in surface and surface["root_BC_type"] == "none": + # no root boundary conditions for jury strut + dof_of_boundary = 0 + root_BC = "none" else: # rigid boundary condition at the root. Constrain translation and rotation. dof_of_boundary = 6 @@ -170,61 +210,160 @@ def setup(self): self.declare_partials(f"disp_aug_{name}", f"local_stiff_transformed_{name}", rows=rows, cols=cols) # --- wing-strut joint constraints (coupling terms between two surfaces) --- - if wing_strut_joint_type == 'ball': - self.n_con = n_con = 3 # number of joint constraints: translational only - elif wing_strut_joint_type == 'pin': - self.n_con = n_con = 5 # translational and rotational in y and z - elif wing_strut_joint_type == 'rigid': - self.n_con = n_con = 6 # ranslational and rotational in x, y, and z - else: - raise ValueError("Joint type must be either 'pin' or 'rigid'.") + # list of all joing Lagrange multipliers (will use later) + self.joint_list = [] + # number of joint constraints + self.n_con = 0 # total number of joint constraints (wing-strut + wing-jury + strut-jury) + n_con_ws = self.n_con_ws = _get_joint_size(wing_strut_joint_type) + self.n_con += n_con_ws # add Lagrange multipliers for joint constraints as state variables - self.add_output("joint_Lag", shape=(n_con,), val=0.0) + self.add_output("joint_Lag_wing_strut", shape=(n_con_ws,), val=0.0) + self.joint_list.append("wing_strut") - k_rows_joint = [] - k_cols_joint = [] - self.k_data_joint = [] + k_rows_joint_ws = [] + k_cols_joint_ws = [] + self.k_data_joint = [] # for all joints (wing-strut, wing-jury, strut-jury) for i, surface in enumerate([self.wing_surface, self.strut_surface]): name = surface["name"] - joint_idx = np.argmin(np.abs(surface["mesh"][0, :, 1] - wing_strut_joint_y)) - joint_indices.append(joint_idx) - print("Surface", name, "joint y:", surface["mesh"][0, joint_idx, 1], "joint index:", joint_idx) - # partials of joint constraint residuals w.r.t. displacement - rows = np.arange(n_con) - if wing_strut_joint_type in ['ball', 'rigid']: - cols = np.arange(n_con) + 6 * joint_idx - elif wing_strut_joint_type == 'pin': - cols = np.array([0, 1, 2, 4, 5]) + 6 * joint_idx # exclude x-rotation + joint_idx = _get_joint_idx_from_y(surface, wing_strut_joint_y) + print(f'Wing-strut joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') + rows, cols = _joint_rows_cols(wing_strut_joint_type, surface, joint_idx) if i == 0: - vals = np.ones(n_con) * 1e9 + vals = np.ones(n_con_ws) * 1e9 else: - vals = np.ones(n_con) * -1e9 - self.declare_partials("joint_Lag", f"disp_aug_{name}", rows=rows, cols=cols, val=vals) + vals = np.ones(n_con_ws) * -1e9 + self.declare_partials("joint_Lag_wing_strut", f"disp_aug_{name}", rows=rows, cols=cols, val=vals) # partials of FEM residuals (r = Ku - f) w.r.t. joint Lagrange multipliers: transpose of above - self.declare_partials(f"disp_aug_{name}", "joint_Lag", rows=cols, cols=rows, val=vals) + self.declare_partials(f"disp_aug_{name}", "joint_Lag_wing_strut", rows=cols, cols=rows, val=vals) - # append these entries to the bottom and right of the entire K matrix - rows += self.size[0] + self.size[1] - if i == 1: + # rows/cols in the global K matrix (bottom of the matrix) + rows += sum(self.size) + if name == "strut": cols += self.size[0] - # partials of joint constraint residuals w.r.t. local stiffness matrix - k_rows_joint.append(rows) - k_cols_joint.append(cols) + # append these entries to the bottom of the global K matrix and its transpose + k_rows_joint_ws.append(rows) + k_cols_joint_ws.append(cols) self.k_data_joint.append(vals) - # partials of FEM residuals wrt joint Lagrange multipliers - k_rows_joint.append(cols) - k_cols_joint.append(rows) + k_rows_joint_ws.append(cols) + k_cols_joint_ws.append(rows) self.k_data_joint.append(vals) - # row and col indices for the entire stiffness matrix (that combines wing, strut, and joint constraints) - self.k_cols = np.concatenate(k_cols + k_cols_joint) - self.k_rows = np.concatenate(k_rows + k_rows_joint) - self.total_size = sum(self.size) + n_con # size of total stiffness matrix (wing + strut + constraints) + print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') + + # row and col indices for the entire stiffness matrix (that combines wing, strut, (jury), and joint constraints) + self.k_cols = np.concatenate(k_cols + k_cols_joint_ws) + self.k_rows = np.concatenate(k_rows + k_rows_joint_ws) + self.total_size = sum(self.size) + n_con_ws # size (num rows) of total stiffness matrix (K matrices + joint constraints) + + print('wing-strut joint added, total_size', self.total_size, '\n') + + if self.include_jury: + # --- wing-jury joint constraints --- + n_con_wj = self.n_con_wj = _get_joint_size(wing_jury_joint_type) + self.n_con += n_con_wj + + # add Lagrange multipliers for joint constraints as state variables + self.add_output("joint_Lag_wing_jury", shape=(n_con_wj,), val=0.0) + self.joint_list.append("wing_jury") + + k_rows_joint_wj = [] + k_cols_joint_wj = [] + for i, surface in enumerate([self.wing_surface, self.jury_surface]): + name = surface["name"] + + # partials of joint constraint residuals w.r.t. displacement + if name == "wing": + joint_idx = _get_joint_idx_from_y(surface, wing_jury_joint_y) + elif name == "jury": + joint_idx = 0 # jury index starts from 0 at the wing-jury joint + print(f'Wing-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') + rows, cols = _joint_rows_cols(wing_jury_joint_type, surface, joint_idx) + if i == 0: + vals = np.ones(n_con_wj) * 1e9 + else: + vals = np.ones(n_con_wj) * -1e9 + self.declare_partials("joint_Lag_wing_jury", f"disp_aug_{name}", rows=rows, cols=cols, val=vals) + + # partials of FEM residuals (r = Ku - f) w.r.t. joint Lagrange multipliers: transpose of above + self.declare_partials(f"disp_aug_{name}", "joint_Lag_wing_jury", rows=cols, cols=rows, val=vals) + + # append these entries to the global K matrix + rows += self.total_size + if name == "jury": + cols += self.size[0] + self.size[1] + + print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') + + k_rows_joint_wj.append(rows) + k_cols_joint_wj.append(cols) + self.k_data_joint.append(vals) + k_rows_joint_wj.append(cols) + k_cols_joint_wj.append(rows) + self.k_data_joint.append(vals) + + # row and col indices for the entire stiffness matrix + self.k_cols = np.concatenate((self.k_cols, np.concatenate(k_cols_joint_wj))) + self.k_rows = np.concatenate((self.k_rows, np.concatenate(k_rows_joint_wj))) + self.total_size += n_con_wj + + print('wing-jury joint added, total_size', self.total_size, '\n') + + # --- strut-jury joint constraints --- + n_con_sj = self.n_con_sj = _get_joint_size(strut_jury_joint_type) + self.n_con += n_con_sj + + # add Lagrange multipliers for joint constraints as state variables + self.add_output("joint_Lag_strut_jury", shape=(n_con_sj,), val=0.0) + self.joint_list.append("strut_jury") + + k_rows_joint_sj = [] + k_cols_joint_sj = [] + for i, surface in enumerate([self.strut_surface, self.jury_surface]): + name = surface["name"] + + # partials of joint constraint residuals w.r.t. displacement + if name == "strut": + joint_idx = _get_joint_idx_from_y(surface, strut_jury_joint_y) + elif name == "jury": + joint_idx = surface["mesh"].shape[1] - 1 # last index of jury surface + print(f'Strut-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') + rows, cols = _joint_rows_cols(strut_jury_joint_type, surface, joint_idx) + if i == 0: + vals = np.ones(n_con_sj) * 1e9 + else: + vals = np.ones(n_con_sj) * -1e9 + self.declare_partials("joint_Lag_strut_jury", f"disp_aug_{name}", rows=rows, cols=cols, val=vals) + + # partials of FEM residuals (r = Ku - f) w.r.t. joint Lagrange multipliers: transpose of above + self.declare_partials(f"disp_aug_{name}", "joint_Lag_strut_jury", rows=cols, cols=rows, val=vals) + + # append these entries to the bottom of the global K matrix (and the transpose to the right) + rows += self.total_size + if name == "strut": + cols += self.size[0] + if name == "jury": + cols += self.size[0] + self.size[1] + + k_rows_joint_sj.append(rows) + k_cols_joint_sj.append(cols) + self.k_data_joint.append(vals) + k_rows_joint_sj.append(cols) + k_cols_joint_sj.append(rows) + self.k_data_joint.append(vals) + + print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') + + # row and col indices for the entire stiffness matrix + self.k_cols = np.concatenate((self.k_cols, np.concatenate(k_cols_joint_sj))) + self.k_rows = np.concatenate((self.k_rows, np.concatenate(k_rows_joint_sj))) + self.total_size += n_con_sj + + print('strut-jury joint added, total_size', self.total_size, '\n') def apply_nonlinear(self, inputs, outputs, residuals): """ @@ -241,17 +380,22 @@ def apply_nonlinear(self, inputs, outputs, residuals): """ K = self.assemble_CSC_K(inputs) disp = np.concatenate([outputs[f"disp_aug_{name}"] for name in self.surface_names]) - joint_Lag = outputs["joint_Lag"] + joint_Lag = np.concatenate([outputs[f"joint_Lag_{name}"] for name in self.joint_list]) u = np.concatenate([disp, joint_Lag]) force = np.concatenate([inputs[f"forces_{name}"] for name in self.surface_names]) rhs = np.concatenate([force, np.zeros(self.n_con)]) r = K.dot(u) - rhs - size0, size1 = self.size[0], self.size[1] # size of each residuals + size0, size1 = self.size[0], self.size[1] # size of each component's K matrix + sum_size = sum(self.size) residuals[f"disp_aug_{self.surface_names[0]}"] = r[:size0] residuals[f"disp_aug_{self.surface_names[1]}"] = r[size0:size0 + size1] - residuals["joint_Lag"] = r[size0 + size1:] + residuals["joint_Lag_wing_strut"] = r[sum_size:sum_size + self.n_con_ws] + if self.include_jury: + residuals[f"disp_aug_{self.surface_names[2]}"] = r[size0 + size1:sum_size] + residuals["joint_Lag_wing_jury"] = r[sum_size + self.n_con_ws:sum_size + self.n_con_ws + self.n_con_wj] + residuals["joint_Lag_strut_jury"] = r[sum_size + self.n_con_ws + self.n_con_wj:] def solve_nonlinear(self, inputs, outputs): """ @@ -266,16 +410,29 @@ def solve_nonlinear(self, inputs, outputs): """ # lu factorization for use with solve_linear K = self.assemble_CSC_K(inputs) + + # Print the matrix K in a nicely formatted way + print('K (dense):') + print(K.toarray()) + self._lup = splu(K) force = np.concatenate([inputs[f"forces_{name}"] for name in self.surface_names]) rhs = np.concatenate([force, np.zeros(self.n_con)]) + print('\nrhs shape', rhs.shape) + print('K shape', K.shape) + u = self._lup.solve(rhs) size0, size1 = self.size[0], self.size[1] # size of each displacement vectors + sum_size = sum(self.size) outputs[f"disp_aug_{self.surface_names[0]}"] = u[:size0] outputs[f"disp_aug_{self.surface_names[1]}"] = u[size0:size0 + size1] - outputs["joint_Lag"] = u[size0 + size1:] + outputs["joint_Lag_wing_strut"] = u[sum_size:sum_size + self.n_con_ws] + if self.include_jury: + outputs[f"disp_aug_{self.surface_names[2]}"] = u[size0 + size1:sum_size] + outputs["joint_Lag_wing_jury"] = u[sum_size + self.n_con_ws:sum_size + self.n_con_ws + self.n_con_wj] + outputs["joint_Lag_strut_jury"] = u[sum_size + self.n_con_ws + self.n_con_wj:] def linearize(self, inputs, outputs, J): """ @@ -290,6 +447,7 @@ def linearize(self, inputs, outputs, J): J : Jacobian sub-jac components written to jacobian[output_name, input_name] """ + raise RuntimeError("Not implemented") vec_size = self.options["vec_size"] name0, name1 = self.surface_names[0], self.surface_names[1] @@ -325,6 +483,7 @@ def solve_linear(self, d_outputs, d_residuals, mode): mode : str either 'fwd' or 'rev' """ + raise RuntimeError("Not implemented") vec_size = self.options["vec_size"] size0, size1 = self.size[0], self.size[1] name0, name1 = self.surface_names[0], self.surface_names[1] @@ -334,31 +493,31 @@ def solve_linear(self, d_outputs, d_residuals, mode): if mode == "fwd": if vec_size > 1: for j in range(vec_size): - rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"][j], d_residuals[f"disp_aug_{name1}"][j], d_residuals["joint_Lag"][j])) + rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"][j], d_residuals[f"disp_aug_{name1}"][j], d_residuals["joint_Lag_wing_strut"][j])) sol = self._lup.solve(rhs) d_outputs[f"disp_aug_{name0}"][j] = sol[:size0] d_outputs[f"disp_aug_{name1}"][j] = sol[size0:size0 + size1] - d_outputs["joint_Lag"][j] = sol[size0 + size1:] + d_outputs["joint_Lag_wing_strut"][j] = sol[size0 + size1:] else: - rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"], d_residuals[f"disp_aug_{name1}"], d_residuals["joint_Lag"])) + rhs = np.concatenate((d_residuals[f"disp_aug_{name0}"], d_residuals[f"disp_aug_{name1}"], d_residuals["joint_Lag_wing_strut"])) sol = self._lup.solve(rhs) d_outputs[f"disp_aug_{name0}"] = sol[:size0] d_outputs[f"disp_aug_{name1}"] = sol[size0:size0 + size1] - d_outputs["joint_Lag"] = sol[size0 + size1:] + d_outputs["joint_Lag_wing_strut"] = sol[size0 + size1:] else: if vec_size > 1: for j in range(vec_size): - rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"][j], d_outputs[f"disp_aug_{name1}"][j], d_outputs["joint_Lag"][j])) + rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"][j], d_outputs[f"disp_aug_{name1}"][j], d_outputs["joint_Lag_wing_strut"][j])) sol = self._lup.solve(rhs) d_residuals[f"disp_aug_{name0}"][j] = sol[:size0] d_residuals[f"disp_aug_{name1}"][j] = sol[size0:size0 + size1] - d_residuals["joint_Lag"][j] = sol[size0 + size1:] + d_residuals["joint_Lag_wing_strut"][j] = sol[size0 + size1:] else: - rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"], d_outputs[f"disp_aug_{name1}"], d_outputs["joint_Lag"])) + rhs = np.concatenate((d_outputs[f"disp_aug_{name0}"], d_outputs[f"disp_aug_{name1}"], d_outputs["joint_Lag_wing_strut"])) sol = self._lup.solve(rhs) d_residuals[f"disp_aug_{name0}"] = sol[:size0] d_residuals[f"disp_aug_{name1}"] = sol[size0:size0 + size1] - d_residuals["joint_Lag"] = sol[size0 + size1:] + d_residuals["joint_Lag_wing_strut"] = sol[size0 + size1:] def assemble_CSC_K(self, inputs): """ @@ -387,6 +546,8 @@ def assemble_CSC_K(self, inputs): elif "root_BC_type" in surface and surface["root_BC_type"] == "pin": # for pin BC (constraint translation and rotation in y and z) data6 = np.full((5,), 1e9) + elif "root_BC_type" in surface and surface["root_BC_type"] == "none": + data6 = np.array([]) else: # for rigid BC (constraint translation and rotation) data6 = np.full((6,), 1e9) diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index 95157a7dc..4edc3ae95 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -1,3 +1,4 @@ +import numpy as np import openmdao.api as om from openaerostruct.structures.create_rhs import CreateRHS from openaerostruct.structures.fem import FEM @@ -113,6 +114,18 @@ def setup(self): # compute RHS force vector for wing and strut, respectively for surf in surface: name = surf["name"] + if name in ["strut", "jury"]: + # no fuel and point mass loads for struts + surf["distributed_fuel_weight"] = False + surf.pop("n_point_masses", None) + + if name == "jury": + # jury is structure-only component so no external loads + ny = surf["mesh"].shape[1] + indep = self.add_subsystem("zero", om.IndepVarComp()) + indep.add_output("jury_loads", val=np.zeros((ny, 6))) + self.connect("zero.jury_loads", f"forces_{name}.loads") + self.add_subsystem(f'forces_{name}', FEMloads(surface=surf), promotes_inputs=["load_factor"]) self.connect(f"forces_{name}.forces", f"fem.forces_{name}") From f7aa2d9f740c64e97b4f6a30cda55b733b278193 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Mon, 4 Aug 2025 18:39:59 -0400 Subject: [PATCH 13/30] fixed strut nodal coordinate computation --- openaerostruct/structures/spatial_beam_setup.py | 5 ++++- openaerostruct/structures/spatial_beam_states.py | 8 +++++++- 2 files changed, 11 insertions(+), 2 deletions(-) diff --git a/openaerostruct/structures/spatial_beam_setup.py b/openaerostruct/structures/spatial_beam_setup.py index 7c4b490d8..cdfc7be35 100644 --- a/openaerostruct/structures/spatial_beam_setup.py +++ b/openaerostruct/structures/spatial_beam_setup.py @@ -12,11 +12,14 @@ class SpatialBeamSetup(om.Group): def initialize(self): self.options.declare("surface", types=dict) + self.options.declare("strut_braced", default=False, types=bool) def setup(self): surface = self.options["surface"] - self.add_subsystem("nodes", ComputeNodes(surface=surface), promotes_inputs=["mesh"], promotes_outputs=["nodes"]) + if not (self.options["strut_braced"] and surface["name"] == "jury"): + # skip node computation for jury strut because no VLM mesh exists for jury strut + self.add_subsystem("nodes", ComputeNodes(surface=surface), promotes_inputs=["mesh"], promotes_outputs=["nodes"]) self.add_subsystem( "assembly", diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index 4edc3ae95..4c6dc2e32 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -119,6 +119,12 @@ def setup(self): surf["distributed_fuel_weight"] = False surf.pop("n_point_masses", None) + # see if we have load_factor input + if surf["struct_weight_relief"] or surf["distributed_fuel_weight"] or "n_point_masses" in surf.keys(): + promotes_inputs = ["load_factor"] + else: + promotes_inputs = [] + if name == "jury": # jury is structure-only component so no external loads ny = surf["mesh"].shape[1] @@ -126,7 +132,7 @@ def setup(self): indep.add_output("jury_loads", val=np.zeros((ny, 6))) self.connect("zero.jury_loads", f"forces_{name}.loads") - self.add_subsystem(f'forces_{name}', FEMloads(surface=surf), promotes_inputs=["load_factor"]) + self.add_subsystem(f'forces_{name}', FEMloads(surface=surf), promotes_inputs=promotes_inputs) self.connect(f"forces_{name}.forces", f"fem.forces_{name}") # FEM of wing-strut system From 086f271f8296a5ed1caf823f292fb8f8de250efe Mon Sep 17 00:00:00 2001 From: kanekosh Date: Mon, 4 Aug 2025 21:21:38 -0400 Subject: [PATCH 14/30] working on jury FEM, something is wrong --- openaerostruct/structures/fem_strut_braced.py | 29 ++----------- openaerostruct/structures/local_stiff.py | 42 +++++++++++-------- .../structures/spatial_beam_states.py | 5 +++ 3 files changed, 33 insertions(+), 43 deletions(-) diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 628cdf698..64d754657 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -145,9 +145,7 @@ def setup(self): wing_strut_joint_y = self.strut_surface["wing_strut_joint_y"] if self.include_jury: - wing_jury_joint_type = self.jury_surface["wing_jury_joint_type"] wing_jury_joint_y = self.jury_surface["wing_jury_joint_y"] - strut_jury_joint_type = self.jury_surface["strut_jury_joint_type"] strut_jury_joint_y = self.jury_surface["strut_jury_joint_y"] # prepare inputs @@ -253,18 +251,14 @@ def setup(self): k_cols_joint_ws.append(rows) self.k_data_joint.append(vals) - print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') - # row and col indices for the entire stiffness matrix (that combines wing, strut, (jury), and joint constraints) self.k_cols = np.concatenate(k_cols + k_cols_joint_ws) self.k_rows = np.concatenate(k_rows + k_rows_joint_ws) self.total_size = sum(self.size) + n_con_ws # size (num rows) of total stiffness matrix (K matrices + joint constraints) - print('wing-strut joint added, total_size', self.total_size, '\n') - if self.include_jury: # --- wing-jury joint constraints --- - n_con_wj = self.n_con_wj = _get_joint_size(wing_jury_joint_type) + n_con_wj = self.n_con_wj = 6 # need to constrain all DOFs to avoid singularity self.n_con += n_con_wj # add Lagrange multipliers for joint constraints as state variables @@ -282,7 +276,7 @@ def setup(self): elif name == "jury": joint_idx = 0 # jury index starts from 0 at the wing-jury joint print(f'Wing-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') - rows, cols = _joint_rows_cols(wing_jury_joint_type, surface, joint_idx) + rows, cols = _joint_rows_cols("rigid", surface, joint_idx) if i == 0: vals = np.ones(n_con_wj) * 1e9 else: @@ -297,8 +291,6 @@ def setup(self): if name == "jury": cols += self.size[0] + self.size[1] - print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') - k_rows_joint_wj.append(rows) k_cols_joint_wj.append(cols) self.k_data_joint.append(vals) @@ -311,10 +303,8 @@ def setup(self): self.k_rows = np.concatenate((self.k_rows, np.concatenate(k_rows_joint_wj))) self.total_size += n_con_wj - print('wing-jury joint added, total_size', self.total_size, '\n') - # --- strut-jury joint constraints --- - n_con_sj = self.n_con_sj = _get_joint_size(strut_jury_joint_type) + n_con_sj = self.n_con_sj = 6 self.n_con += n_con_sj # add Lagrange multipliers for joint constraints as state variables @@ -332,7 +322,7 @@ def setup(self): elif name == "jury": joint_idx = surface["mesh"].shape[1] - 1 # last index of jury surface print(f'Strut-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') - rows, cols = _joint_rows_cols(strut_jury_joint_type, surface, joint_idx) + rows, cols = _joint_rows_cols("rigid", surface, joint_idx) if i == 0: vals = np.ones(n_con_sj) * 1e9 else: @@ -356,15 +346,11 @@ def setup(self): k_cols_joint_sj.append(rows) self.k_data_joint.append(vals) - print(f'name = {name}, rows = {rows}, cols = {cols}, vals = {vals}') - # row and col indices for the entire stiffness matrix self.k_cols = np.concatenate((self.k_cols, np.concatenate(k_cols_joint_sj))) self.k_rows = np.concatenate((self.k_rows, np.concatenate(k_rows_joint_sj))) self.total_size += n_con_sj - print('strut-jury joint added, total_size', self.total_size, '\n') - def apply_nonlinear(self, inputs, outputs, residuals): """ R = Ax - b. @@ -411,17 +397,10 @@ def solve_nonlinear(self, inputs, outputs): # lu factorization for use with solve_linear K = self.assemble_CSC_K(inputs) - # Print the matrix K in a nicely formatted way - print('K (dense):') - print(K.toarray()) - self._lup = splu(K) force = np.concatenate([inputs[f"forces_{name}"] for name in self.surface_names]) rhs = np.concatenate([force, np.zeros(self.n_con)]) - print('\nrhs shape', rhs.shape) - print('K shape', K.shape) - u = self._lup.solve(rhs) size0, size1 = self.size[0], self.size[1] # size of each displacement vectors diff --git a/openaerostruct/structures/local_stiff.py b/openaerostruct/structures/local_stiff.py index f4b5f321e..3a26d4c37 100644 --- a/openaerostruct/structures/local_stiff.py +++ b/openaerostruct/structures/local_stiff.py @@ -32,6 +32,7 @@ class LocalStiff(om.ExplicitComponent): def initialize(self): self.options.declare("surface", types=dict) + self.options.declare("strut_braced", default=False, types=bool) def setup(self): surface = self.options["surface"] @@ -55,6 +56,11 @@ def setup(self): self.declare_partials("local_stiff", "Iz", rows=rows, cols=cols) self.declare_partials("local_stiff", "element_lengths", rows=rows, cols=cols) + self.scale_EI = 1.0 + if self.options["strut_braced"] and surface["name"] == "jury": + # for jury strut, we set very low bending stiffness to effectively make this a rod element + self.scale_EI = 1.0 + def compute(self, inputs, outputs): surface = self.options["surface"] @@ -76,8 +82,8 @@ def compute(self, inputs, outputs): for i in range(4): for j in range(4): - outputs["local_stiff"][:, 4 + i, 4 + j] = E * Iy / L**3 * coeffs_y[i, j] - outputs["local_stiff"][:, 8 + i, 8 + j] = E * Iz / L**3 * coeffs_z[i, j] + outputs["local_stiff"][:, 4 + i, 4 + j] = E * Iy / L**3 * coeffs_y[i, j] * self.scale_EI + outputs["local_stiff"][:, 8 + i, 8 + j] = E * Iz / L**3 * coeffs_z[i, j] * self.scale_EI for i in [1, 3]: for j in range(4): @@ -122,31 +128,31 @@ def compute_partials(self, inputs, partials): for i in range(4): for j in range(4): - derivs_Iy[:, 4 + i, 4 + j] = E / L**3 * coeffs_y[i, j] - derivs_L[:, 4 + i, 4 + j] = -3 * E * Iy / L**4 * coeffs_y[i, j] + derivs_Iy[:, 4 + i, 4 + j] = E / L**3 * coeffs_y[i, j] * self.scale_EI + derivs_L[:, 4 + i, 4 + j] = -3 * E * Iy / L**4 * coeffs_y[i, j] * self.scale_EI - derivs_Iz[:, 8 + i, 8 + j] = E / L**3 * coeffs_z[i, j] - derivs_L[:, 8 + i, 8 + j] = -3 * E * Iz / L**4 * coeffs_z[i, j] + derivs_Iz[:, 8 + i, 8 + j] = E / L**3 * coeffs_z[i, j] * self.scale_EI + derivs_L[:, 8 + i, 8 + j] = -3 * E * Iz / L**4 * coeffs_z[i, j] * self.scale_EI for i in [1, 3]: for j in range(4): - derivs_Iy[:, 4 + i, 4 + j] = E / L**2 * coeffs_y[i, j] - derivs_L[:, 4 + i, 4 + j] = -2 * E * Iy / L**3 * coeffs_y[i, j] + derivs_Iy[:, 4 + i, 4 + j] = E / L**2 * coeffs_y[i, j] * self.scale_EI + derivs_L[:, 4 + i, 4 + j] = -2 * E * Iy / L**3 * coeffs_y[i, j] * self.scale_EI - derivs_Iz[:, 8 + i, 8 + j] = E / L**2 * coeffs_z[i, j] - derivs_L[:, 8 + i, 8 + j] = -2 * E * Iz / L**3 * coeffs_z[i, j] + derivs_Iz[:, 8 + i, 8 + j] = E / L**2 * coeffs_z[i, j] * self.scale_EI + derivs_L[:, 8 + i, 8 + j] = -2 * E * Iz / L**3 * coeffs_z[i, j] * self.scale_EI for i in range(4): for j in [1, 3]: - derivs_Iy[:, 4 + i, 4 + j] = E / L**2 * coeffs_y[i, j] - derivs_L[:, 4 + i, 4 + j] = -2 * E * Iy / L**3 * coeffs_y[i, j] + derivs_Iy[:, 4 + i, 4 + j] = E / L**2 * coeffs_y[i, j] * self.scale_EI + derivs_L[:, 4 + i, 4 + j] = -2 * E * Iy / L**3 * coeffs_y[i, j] * self.scale_EI - derivs_Iz[:, 8 + i, 8 + j] = E / L**2 * coeffs_z[i, j] - derivs_L[:, 8 + i, 8 + j] = -2 * E * Iz / L**3 * coeffs_z[i, j] + derivs_Iz[:, 8 + i, 8 + j] = E / L**2 * coeffs_z[i, j] * self.scale_EI + derivs_L[:, 8 + i, 8 + j] = -2 * E * Iz / L**3 * coeffs_z[i, j] * self.scale_EI for i in [1, 3]: for j in [1, 3]: - derivs_Iy[:, 4 + i, 4 + j] = E / L * coeffs_y[i, j] - derivs_L[:, 4 + i, 4 + j] = -E * Iy / L**2 * coeffs_y[i, j] + derivs_Iy[:, 4 + i, 4 + j] = E / L * coeffs_y[i, j] * self.scale_EI + derivs_L[:, 4 + i, 4 + j] = -E * Iy / L**2 * coeffs_y[i, j] * self.scale_EI - derivs_Iz[:, 8 + i, 8 + j] = E / L * coeffs_z[i, j] - derivs_L[:, 8 + i, 8 + j] = -E * Iz / L**2 * coeffs_z[i, j] + derivs_Iz[:, 8 + i, 8 + j] = E / L * coeffs_z[i, j] * self.scale_EI + derivs_L[:, 8 + i, 8 + j] = -E * Iz / L**2 * coeffs_z[i, j] * self.scale_EI diff --git a/openaerostruct/structures/spatial_beam_states.py b/openaerostruct/structures/spatial_beam_states.py index 4c6dc2e32..265c8f0d4 100644 --- a/openaerostruct/structures/spatial_beam_states.py +++ b/openaerostruct/structures/spatial_beam_states.py @@ -111,6 +111,11 @@ def setup(self): self.add_subsystem("disp", Disp(surface=surface), promotes_inputs=["*"], promotes_outputs=["*"]) else: + # set no root boundary conditions for jury strut + for i, surf in enumerate(surface): + if surf["name"] == "jury": + surface[i]["root_BC_type"] = "none" + # compute RHS force vector for wing and strut, respectively for surf in surface: name = surf["name"] From 9f87f03df686c503db81de653bfb7ccca7eeea5e Mon Sep 17 00:00:00 2001 From: kanekosh Date: Fri, 8 Aug 2025 00:35:59 -0400 Subject: [PATCH 15/30] jury strut works for aerostructural analysis --- .../integration/aerostruct_groups.py | 109 +++++++++-- openaerostruct/structures/disp.py | 3 + .../structures/failure_buckling_ks.py | 27 ++- openaerostruct/structures/fem_strut_braced.py | 19 +- .../structures/spatial_beam_functionals.py | 11 +- .../spatial_beam_states_truss_braced.py | 181 ++++++++++++++++++ 6 files changed, 309 insertions(+), 41 deletions(-) create mode 100644 openaerostruct/structures/spatial_beam_states_truss_braced.py diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 8ecd2f1f6..95e230bb0 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -1,8 +1,11 @@ from openaerostruct.aerodynamics.geometry import VLMGeometry from openaerostruct.geometry.geometry_group import Geometry from openaerostruct.transfer.displacement_transfer_group import DisplacementTransferGroup +from openaerostruct.structures.failure_buckling_ks import EulerColumnBucklingFailureKS +from openaerostruct.structures.fem_strut_braced import get_joint_idx_from_y from openaerostruct.structures.spatial_beam_setup import SpatialBeamSetup from openaerostruct.structures.spatial_beam_states import SpatialBeamStates +from openaerostruct.structures.spatial_beam_states_truss_braced import SpatialBeamStatesTrussBraced from openaerostruct.aerodynamics.functionals import VLMFunctionals from openaerostruct.structures.spatial_beam_functionals import SpatialBeamFunctionals from openaerostruct.functionals.total_performance import TotalPerformance @@ -190,9 +193,6 @@ def setup(self): if "panel_buckling" in surface and surface["panel_buckling"]: promotes_inputs += ["skin_thickness", "t_over_c", "fem_chords"] promotes_outputs += ["failure_local_buckling"] - if "column_buckling" in surface and surface["column_buckling"]: - promotes_inputs += ["Iz", "joint_load"] - promotes_outputs += ["failure_column_buckling"] self.add_subsystem( "struct_funcs", @@ -235,6 +235,17 @@ def setup(self): pass else: raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + + # add coupled FEM group + coupled.add_subsystem( + "FEM_SBW", + SpatialBeamStates(surface=surfaces_all, strut_braced=True), + promotes_inputs=["*"], + promotes_outputs=["disp_*"] + ) + + include_jury = False + elif len(surfaces) == 3: # wing + strut + jury if surfaces[0]["name"] == "wing" and surfaces[1]["name"] == "strut" and surfaces[2]["name"] == "jury": @@ -248,19 +259,22 @@ def setup(self): surfaces[2]["root_BC_type"] = "none" surfaces_all[2]["root_BC_type"] = "none" - # add coupled FEM group - coupled.add_subsystem( - "FEM_SBW", - SpatialBeamStates(surface=surfaces_all, strut_braced=True), - promotes_inputs=["*"], - promotes_outputs=["disp_*"] - ) + # add coupled FEM group + coupled.add_subsystem( + "FEM_SBW", + SpatialBeamStatesTrussBraced(surfaces=surfaces_all, strut_braced=True), + promotes_inputs=["*"], + promotes_outputs=["disp_*"] + ) + + include_jury = True + surface_jury = surfaces[2] else: surfaces_all = surfaces surfaces_AS = surfaces # connections for all structural components - for surface in surfaces_all: + for surface in surfaces_AS: name = surface["name"] # connect displacements to the performance group @@ -274,7 +288,13 @@ def setup(self): # Connect the output of the loads component with the FEM # displacement parameter. This links the coupling within the coupled # group that necessitates the subgroup solver. - loads_target = f"forces_{name}.loads" if self.options["strut_braced"] else f"{name}.loads" + if self.options["strut_braced"] and include_jury: + loads_target = f"external_loads_{name}.loads" + elif self.options["strut_braced"]: + loads_target = f"forces_{name}.loads" + else: + loads_target = f"{name}.loads" + coupled.connect(name + "_loads.loads", loads_target) # Connect displacement from FEM to transfer model. This is automatically done by promotion for non-SBW cases @@ -404,12 +424,71 @@ def setup(self): name + "_perf", perf_group, promotes_inputs=["rho", "v", "alpha", "beta", "re", "Mach_number"] ) + # column buckling failure for strut + if "column_buckling" in surface and surface["column_buckling"]: + if include_jury: + ny = surface["mesh"].shape[1] + passthru = om.ExecComp( + ["nodes_all = nodes", "Iz_all = Iz"], + nodes_all={'units': 'm', 'shape': (ny, 3)}, + nodes={'units': 'm', 'shape': (ny, 3)}, + Iz_all={'units': 'm**4', 'shape': ny - 1}, + Iz={'units': 'm**4', 'shape': ny - 1} + ) + perf_group.add_subsystem("passthru", passthru, promotes_inputs=["nodes", "Iz"]) + + # strut-jury joint index + joint_idx = get_joint_idx_from_y(surface, surface_jury["strut_jury_joint_y"]) + ny_outer = joint_idx + 1 # number of nodes from wing to strut-jury joint, including the joint node + ny_inner = ny - joint_idx # number of nodes from strut-jury joint to root, including the joint node + + # column buckling for outer strut segment. Use strut-wing joint load as a compression load + perf_group.add_subsystem( + "column_buckling_outer", + EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_outer), + promotes_inputs=[("joint_load", "joint_load_strut_wing")], + promotes_outputs=[("failure_column_buckling", "failure_column_buckling_outer")] + ) + perf_group.connect("passthru.nodes_all", "column_buckling_outer.nodes", src_indices=om.slicer[0:ny_outer, :]) + perf_group.connect("passthru.Iz_all", "column_buckling_outer.Iz", src_indices=om.slicer[0:ny_outer - 1]) + + # column buckling for inner strut segment. Use strut-root reaction force as a compression load + perf_group.add_subsystem( + "column_buckling_inner", + EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_inner), + promotes_inputs=[("joint_load", "reaction_load_strut_root")], + promotes_outputs=[("failure_column_buckling", "failure_column_buckling_inner")] + ) + perf_group.connect("passthru.nodes_all", "column_buckling_inner.nodes", src_indices=om.slicer[joint_idx:, :]) + perf_group.connect("passthru.Iz_all", "column_buckling_inner.Iz", src_indices=om.slicer[joint_idx:]) + else: + # column buckling for the entire strut. Use strut-wing joint load as a compression load + ny = surface["mesh"].shape[1] + perf_group.add_subsystem( + "column_buckling", + EulerColumnBucklingFailureKS(surface=surface, ny=ny, joint_load_type="vector"), + promotes_inputs=["nodes", "Iz", ("joint_load", "joint_load_strut_wing")], + promotes_outputs=["failure_column_buckling"] + ) + # post-AS-analysis computations for structure-only components + # for surface in [s for s in surfaces_all if s not in surfaces_AS]: + # self.add_subsystem( + # f"{surface["name"]}_perf", + # SpatialBeamFunctionals(surface=surface), + # ) + + # column buckling failure for jury strut for surface in [s for s in surfaces_all if s not in surfaces_AS]: - self.add_subsystem( - f"{surface["name"]}_perf", - SpatialBeamFunctionals(surface=surface), + ny = surface["mesh"].shape[1] + perf_group = self.add_subsystem(f"{surface["name"]}_perf", om.Group()) + perf_group.add_subsystem( + "column_buckling", + EulerColumnBucklingFailureKS(surface=surface, ny=ny, joint_load_type="scaler"), + promotes_inputs=["nodes", "Iz", "joint_load"], + promotes_outputs=["failure_column_buckling"] ) + self.connect("coupled.FEM_SBW.joint_axias_force", f"{surface["name"]}_perf.joint_load") # Add functionals to evaluate performance of the system. # Note that only the interesting results are promoted here; not all diff --git a/openaerostruct/structures/disp.py b/openaerostruct/structures/disp.py index 60f1052ec..f6804583e 100644 --- a/openaerostruct/structures/disp.py +++ b/openaerostruct/structures/disp.py @@ -47,12 +47,15 @@ def setup(self): self.add_input("disp_aug", val=np.zeros((self.ny * 6 + dof_of_boundary)), units="m") self.add_output("disp", val=np.zeros((self.ny, 6)), units="m") + self.add_output("boundary_Lag", val=np.zeros(dof_of_boundary)) # Lagrange multipliers for boundary conditions n = self.ny * 6 arange = np.arange((n)) self.declare_partials("disp", "disp_aug", val=1.0, rows=arange, cols=arange) + self.declare_partials("boundary_Lag", "disp_aug", val=1.0, rows=np.arange(dof_of_boundary), cols=np.arange(dof_of_boundary) + n) def compute(self, inputs, outputs): # Obtain the relevant portions of disp_aug and store the reshaped # displacements in disp outputs["disp"] = inputs["disp_aug"][:self.ny * 6].reshape((-1, 6)) + outputs["boundary_Lag"] = inputs["disp_aug"][self.ny * 6:] * 1.0 diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py index 881dbdb83..b2f8dcadf 100644 --- a/openaerostruct/structures/failure_buckling_ks.py +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -146,8 +146,8 @@ class EulerColumnBucklingFailureKS(om.ExplicitComponent): ---------- nodes[ny, 3] : numpy array FEM node coordinates. - joint_load[3] : numpy array - Force vector applied to the strut at the joint. + joint_load : numpy array or float + Force vector applied to the strut at the joint. If `joint_load_type` is `vector`, this is a 3-element vector. If `joint_load_type` is `scaler`, this is a scalar tensile force. Iz[ny-1] : numpy array Moment of inertia of the strut about the z-axis. @@ -160,14 +160,20 @@ class EulerColumnBucklingFailureKS(om.ExplicitComponent): def initialize(self): self.options.declare("surface", types=dict) + self.options.declare("ny", desc="number of nodes in the strut") self.options.declare("rho", types=float, default=100.0, desc="KS aggregation smoothness parameter") + self.options.declare("joint_load_type", types=str, default="vector", desc="type of joint load, either `vector` or `scaler`. If scaler, this should be the axial tension force") def setup(self): - surface = self.options["surface"] - ny = surface["mesh"].shape[1] + ny = self.options["ny"] self.add_input("nodes", shape=(ny, 3), units="m", desc="FEM node coordinates") - self.add_input("joint_load", shape=3, desc="load vector applied to the strut at the joint") + if self.options["joint_load_type"] == "vector": + self.add_input("joint_load", shape=3, desc="load vector applied to the strut at the joint") + elif self.options["joint_load_type"] == "scaler": + self.add_input("joint_load", shape=1, desc="tension force applied to the strut at the joint") + else: + raise ValueError(f"Invalid joint load type: {self.options['joint_load_type']}") self.add_input("Iz", shape=(ny - 1), units="m**4") # use Iz because Iz < Iy self.add_output("failure_column_buckling", val=0.0, desc="Euler column buckling failure metric, should be <= 0") @@ -192,10 +198,13 @@ def compute(self, inputs, outputs): strut_len = np.linalg.norm(nodes[0, :] - nodes[-1, :]) * column_length_factor # compression load - strut_direction = (nodes[-1, :] - nodes[0, :]) - strut_direction /= np.linalg.norm(strut_direction) # unit magniture - # multiply by 1e9 to make the units to N (FEM normalizes this by 1e9) - comp_load = np.dot(joint_load, strut_direction) * 1e9 + if self.options["joint_load_type"] == "vector": + strut_direction = (nodes[-1, :] - nodes[0, :]) + strut_direction /= np.linalg.norm(strut_direction) # unit magniture + # multiply by 1e9 to make the units to N (FEM normalizes this by 1e9) + comp_load = np.dot(joint_load, strut_direction) * 1e9 + elif self.options["joint_load_type"] == "scaler": + comp_load = -joint_load # flip sign to make compression positive # critical load with safety factor of 1.5 sf = 1.5 diff --git a/openaerostruct/structures/fem_strut_braced.py b/openaerostruct/structures/fem_strut_braced.py index 64d754657..b3c3526a2 100644 --- a/openaerostruct/structures/fem_strut_braced.py +++ b/openaerostruct/structures/fem_strut_braced.py @@ -23,13 +23,13 @@ def _get_joint_size(joint_type): raise ValueError("Joint type must be either 'pin' or 'rigid'.") return n_con -def _get_joint_idx_from_y(surface, joint_y): +def get_joint_idx_from_y(surface, joint_y): """ Return the joint index from the joint y (spanwise) coordinate. This is only good for wing or strut, but do not use for jury (or any vertical surface) """ if surface["name"] == "jury": - raise RuntimeError("Cannot use _get_joint_idx_from_y for jury surface.") + raise RuntimeError("Cannot use get_joint_idx_from_y for jury surface.") joint_idx = np.argmin((np.abs(surface["mesh"][0, :, 1]) - abs(joint_y))**2) return joint_idx @@ -127,6 +127,7 @@ def setup(self): self.wing_surface = self.surfaces[0] self.strut_surface = self.surfaces[1] self.jury_surface = self.surfaces[2] + raise NotImplementedError("Jury strut FEM model is still under debugging") else: raise ValueError("FEMStrutBraced component requires exactly two or three surfaces (wing, strut, and/or jury).") @@ -226,7 +227,7 @@ def setup(self): name = surface["name"] # partials of joint constraint residuals w.r.t. displacement - joint_idx = _get_joint_idx_from_y(surface, wing_strut_joint_y) + joint_idx = get_joint_idx_from_y(surface, wing_strut_joint_y) print(f'Wing-strut joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') rows, cols = _joint_rows_cols(wing_strut_joint_type, surface, joint_idx) if i == 0: @@ -272,7 +273,7 @@ def setup(self): # partials of joint constraint residuals w.r.t. displacement if name == "wing": - joint_idx = _get_joint_idx_from_y(surface, wing_jury_joint_y) + joint_idx = get_joint_idx_from_y(surface, wing_jury_joint_y) elif name == "jury": joint_idx = 0 # jury index starts from 0 at the wing-jury joint print(f'Wing-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') @@ -318,7 +319,7 @@ def setup(self): # partials of joint constraint residuals w.r.t. displacement if name == "strut": - joint_idx = _get_joint_idx_from_y(surface, strut_jury_joint_y) + joint_idx = get_joint_idx_from_y(surface, strut_jury_joint_y) elif name == "jury": joint_idx = surface["mesh"].shape[1] - 1 # last index of jury surface print(f'Strut-jury joint ({surface["name"]}): joint y = {surface["mesh"][0, joint_idx, 1]} and joint index = {joint_idx}') @@ -426,7 +427,9 @@ def linearize(self, inputs, outputs, J): J : Jacobian sub-jac components written to jacobian[output_name, input_name] """ - raise RuntimeError("Not implemented") + if self.include_jury: + raise RuntimeError("Not implemented for jury strut") + vec_size = self.options["vec_size"] name0, name1 = self.surface_names[0], self.surface_names[1] @@ -462,7 +465,9 @@ def solve_linear(self, d_outputs, d_residuals, mode): mode : str either 'fwd' or 'rev' """ - raise RuntimeError("Not implemented") + if self.include_jury: + raise RuntimeError("Not implemented for jury strut") + vec_size = self.options["vec_size"] size0, size1 = self.size[0], self.size[1] name0, name1 = self.surface_names[0], self.surface_names[1] diff --git a/openaerostruct/structures/spatial_beam_functionals.py b/openaerostruct/structures/spatial_beam_functionals.py index 705ad5163..5dca4f64a 100644 --- a/openaerostruct/structures/spatial_beam_functionals.py +++ b/openaerostruct/structures/spatial_beam_functionals.py @@ -8,7 +8,7 @@ from openaerostruct.structures.non_intersecting_thickness import NonIntersectingThickness from openaerostruct.structures.failure_exact import FailureExact from openaerostruct.structures.failure_ks import FailureKS -from openaerostruct.structures.failure_buckling_ks import PanelLocalBucklingFailureKS, EulerColumnBucklingFailureKS +from openaerostruct.structures.failure_buckling_ks import PanelLocalBucklingFailureKS class SpatialBeamFunctionals(om.Group): @@ -89,12 +89,3 @@ def setup(self): self.connect("vonmises.lower_skin_comp_stress", "local_buckling.lower_skin_comp_stress") self.connect("vonmises.front_spar_shear_stress", "local_buckling.front_spar_shear_stress") self.connect("vonmises.rear_spar_shear_stress", "local_buckling.rear_spar_shear_stress") - - # global Euler column buckling - if "column_buckling" in surface and surface["column_buckling"]: - self.add_subsystem( - "column_buckling", - EulerColumnBucklingFailureKS(surface=surface), - promotes_inputs=["nodes", "joint_load", "Iz"], - promotes_outputs=["failure_column_buckling"] - ) diff --git a/openaerostruct/structures/spatial_beam_states_truss_braced.py b/openaerostruct/structures/spatial_beam_states_truss_braced.py new file mode 100644 index 000000000..794aa51a5 --- /dev/null +++ b/openaerostruct/structures/spatial_beam_states_truss_braced.py @@ -0,0 +1,181 @@ +import numpy as np +import openmdao.api as om +from openaerostruct.structures.create_rhs import CreateRHS +from openaerostruct.structures.fem import FEM +from openaerostruct.structures.fem_strut_braced import FEMStrutBraced, get_joint_idx_from_y +from openaerostruct.structures.disp import Disp +from openaerostruct.structures.wing_weight_loads import StructureWeightLoads +from openaerostruct.structures.fuel_loads import FuelLoads +from openaerostruct.structures.total_loads import TotalLoads +from openaerostruct.structures.compute_point_mass_loads import ComputePointMassLoads +from openaerostruct.structures.compute_thrust_loads import ComputeThrustLoads +from openaerostruct.structures.spatial_beam_states import FEMloads + + +class SpatialBeamStatesTrussBraced(om.Group): + """ + Group that contains the spatial beam states for truss-braced wing model + that combines FEM for wing and main strut + analytic rod for jury strut + """ + + def initialize(self): + self.options.declare("surfaces", types=list) + self.options.declare("strut_braced", default=True, types=bool) + + def setup(self): + surfaces = self.options["surfaces"] + + if len(surfaces) != 3: + raise ValueError("surface must have 3 surfaces: wing, strut, and jury") + + surf_wing = surfaces[0] + surf_strut = surfaces[1] + surf_jury = surfaces[2] + + # no fuel and point mass loads for struts and jury + surf_strut["distributed_fuel_weight"] = False + surf_strut.pop("n_point_masses", None) + surf_jury["distributed_fuel_weight"] = False + surf_jury.pop("n_point_masses", None) + + # find spanwise index of jury-wing and jury-strut joints + jury_joint_idx = { + "wing": get_joint_idx_from_y(surf_wing, surf_jury["wing_jury_joint_y"]), + "strut": get_joint_idx_from_y(surf_strut, surf_jury["strut_jury_joint_y"]), + } + + # --- Strut-braced wing FEM --- + # compute RHS force vector for wing and strut + for surf in [surf_wing, surf_strut]: + name = surf["name"] + + # compute jury joint loads + self.add_subsystem(f"external_loads_{name}", JuryJointLoads(surface=surf, joint_idx=jury_joint_idx[name]), promotes_inputs=["jury_nodes", "joint_axias_force"]) + + # see if we have load_factor input + if surf["struct_weight_relief"] or surf["distributed_fuel_weight"] or "n_point_masses" in surf.keys(): + promotes_inputs = ["load_factor"] + else: + promotes_inputs = [] + + # RHS vector + self.add_subsystem(f'forces_{name}', FEMloads(surface=surf), promotes_inputs=promotes_inputs) + self.connect(f"external_loads_{name}.loads_out", f"forces_{name}.loads") + self.connect(f"forces_{name}.forces", f"fem.forces_{name}") + + # FEM of wing-strut system + self.add_subsystem("fem", FEMStrutBraced(surfaces=[surf_wing, surf_strut]), promotes_inputs=["local_stiff_transformed_*"]) + + # reshape displacements and remove Lagrange multipliers from disp_aug + for surf in [surf_wing, surf_strut]: + name = surf["name"] + self.add_subsystem(f"disp_{name}", Disp(surface=surf), promotes_outputs=[("disp", f"disp_{name}")]) + self.connect(f"fem.disp_aug_{name}", f"disp_{name}.disp_aug") + + # --- compute jury axial displacement from strut-braced wing FEM solution --- + self.add_subsystem("jury_disp_from_SBW", JuryDispFromSBW(), promotes_inputs=["jury_nodes"], promotes_outputs=["disp_jury"]) + self.connect("disp_wing", "jury_disp_from_SBW.disp_wing", src_indices=om.slicer[jury_joint_idx["wing"], :3]) + self.connect("disp_strut", "jury_disp_from_SBW.disp_strut", src_indices=om.slicer[jury_joint_idx["strut"], :3]) + + # --- Compute jury axial load from displacement --- + self.add_subsystem("jury_load", JuryStrutRod(surface=surf_jury), promotes_inputs=["jury_nodes", "jury_A", "disp_jury"], promotes_outputs=["joint_axias_force"]) + + self.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, maxiter=20, atol=1e-6, rtol=1e-6, err_on_non_converge=True, iprint=0) + self.linear_solver = om.DirectSolver() + + +class JuryJointLoads(om.ExplicitComponent): + """ + Compute jury joint loads + """ + + def initialize(self): + self.options.declare("surface", types=dict) + self.options.declare("joint_idx", desc="spanwise index of jury joint on this surface") + + def setup(self): + surface = self.options["surface"] + ny = surface["mesh"].shape[1] + + self.add_input("jury_nodes", val=np.zeros((2, 3)), units='m') + self.add_input("joint_axias_force", val=0., units="N", desc="axial force on jury joint. Positive = tension") + self.add_input("loads", val=np.ones((ny, 6)), units="N") + self.add_output("loads_out", val=np.ones((ny, 6)), units="N") + + self.declare_partials("loads_out", ["jury_nodes", "joint_axias_force"], method="cs", step=1e-100) + self.declare_partials("loads_out", "loads", rows=np.arange(ny * 6), cols=np.arange(ny * 6), val=np.ones(ny * 6)) + + def compute(self, inputs, outputs): + surface = self.options["surface"] + joint_idx = self.options["joint_idx"] + + jury_nodes = inputs["jury_nodes"] + + # direction of jury struts = axial force. Jury node is defined from wing to node + if surface["name"] == "wing": + load_dir = jury_nodes[-1, :] - jury_nodes[0, :] + elif surface["name"] == "strut": + load_dir = jury_nodes[0, :] - jury_nodes[-1, :] + else: + raise ValueError(f"Jury joint loads are only supported for wing and strut surfaces, not {surface['name']}") + load_dir = load_dir / np.linalg.norm(load_dir) + + # compute load vector and add to the loads + load_vec = inputs["joint_axias_force"] * load_dir + outputs["loads_out"] = inputs["loads"] * 1.0 + outputs["loads_out"][joint_idx, :3] += load_vec + + +class JuryStrutRod(om.ExplicitComponent): + """ + Analytic rod model for jury strut + """ + + def initialize(self): + self.options.declare("surface", types=dict) + + def setup(self): + self.add_input("jury_nodes", val=np.zeros((2, 3)), units='m') + self.add_input("jury_A", val=1., units="m**2", desc="cross-sectional area of jury strut") + self.add_input("disp_jury", val=0., units='m', desc="axial displacement of jury strut (tension = positive)") + + self.add_output("joint_axias_force", val=0., units="N", desc="axial force on jury joint. Positive = tension") + + self.declare_partials("*", "*", method="cs", step=1e-100) + + def compute(self, inputs, outputs): + surface = self.options["surface"] + E = surface["E"] + + A = inputs["jury_A"] + jury_nodes = inputs["jury_nodes"] + disp = inputs["disp_jury"] + + jury_len = np.linalg.norm(jury_nodes[-1, :] - jury_nodes[0, :]) + outputs["joint_axias_force"] = E * A * disp / jury_len + + +class JuryDispFromSBW(om.ExplicitComponent): + """ + Compute jury strut displacement from the strut-braced wing FEM displacements + """ + def setup(self): + self.add_input("disp_wing", shape=(3,), units='m') + self.add_input("disp_strut", shape=(3,), units='m') + self.add_input("jury_nodes", shape=(2, 3), units='m') + + self.add_output("disp_jury", val=0., units='m', desc="axial displacement of jury strut (tension = positive)") + + self.declare_partials("*", "*", method="cs", step=1e-100) + + def compute(self, inputs, outputs): + jury_nodes = inputs["jury_nodes"] + + # jury strut direction + jury_dir = jury_nodes[-1, :] - jury_nodes[0, :] + jury_dir = jury_dir / np.linalg.norm(jury_dir) + + # compute jury strut tensile displacement + disp_wing_mapped = np.dot(inputs["disp_wing"], jury_dir) + disp_strut_mapped = np.dot(inputs["disp_strut"], jury_dir) + outputs["disp_jury"] = disp_strut_mapped - disp_wing_mapped \ No newline at end of file From c513559095fd0189c466387949a93865f00cfff2 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Mon, 18 Aug 2025 14:57:57 -0400 Subject: [PATCH 16/30] bugfix for conventional --- openaerostruct/functionals/total_performance.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/openaerostruct/functionals/total_performance.py b/openaerostruct/functionals/total_performance.py index d14b0e7e8..2df5ba0cf 100644 --- a/openaerostruct/functionals/total_performance.py +++ b/openaerostruct/functionals/total_performance.py @@ -40,6 +40,9 @@ def setup(self): raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") else: raise ValueError("surfaces must be in order of [wing_surface, strut_surface, (optional) jury_surface]") + else: + surfaces_all = surfaces + surfaces_AS = surfaces if not self.options["user_specified_Sref"]: self.add_subsystem( From 0fe4a718df2d3ccd1ea75d52b76291135ef97db0 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Mon, 18 Aug 2025 17:20:40 -0400 Subject: [PATCH 17/30] minor fix for old python --- openaerostruct/integration/aerostruct_groups.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 95e230bb0..df07afee6 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -488,7 +488,7 @@ def setup(self): promotes_inputs=["nodes", "Iz", "joint_load"], promotes_outputs=["failure_column_buckling"] ) - self.connect("coupled.FEM_SBW.joint_axias_force", f"{surface["name"]}_perf.joint_load") + self.connect("coupled.FEM_SBW.joint_axias_force", surface["name"] + "_perf.joint_load") # Add functionals to evaluate performance of the system. # Note that only the interesting results are promoted here; not all From c609a34d875a33a4553dbad0954217e98dc29080 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Mon, 18 Aug 2025 17:28:54 -0400 Subject: [PATCH 18/30] more very minor fix --- openaerostruct/integration/aerostruct_groups.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index df07afee6..8357b95b7 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -474,14 +474,14 @@ def setup(self): # post-AS-analysis computations for structure-only components # for surface in [s for s in surfaces_all if s not in surfaces_AS]: # self.add_subsystem( - # f"{surface["name"]}_perf", + # surface["name"] + "_perf", # SpatialBeamFunctionals(surface=surface), # ) # column buckling failure for jury strut for surface in [s for s in surfaces_all if s not in surfaces_AS]: ny = surface["mesh"].shape[1] - perf_group = self.add_subsystem(f"{surface["name"]}_perf", om.Group()) + perf_group = self.add_subsystem(surface["name"] + "_perf", om.Group()) perf_group.add_subsystem( "column_buckling", EulerColumnBucklingFailureKS(surface=surface, ny=ny, joint_load_type="scaler"), From 8fcacb29c17c5a07b16b1af334ad9e818767a4d5 Mon Sep 17 00:00:00 2001 From: kanekosh Date: Tue, 26 Aug 2025 14:02:57 -0400 Subject: [PATCH 19/30] added column length factor for Euler buckling --- openaerostruct/integration/aerostruct_groups.py | 6 +++--- openaerostruct/structures/failure_buckling_ks.py | 3 ++- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/openaerostruct/integration/aerostruct_groups.py b/openaerostruct/integration/aerostruct_groups.py index 8357b95b7..d05e3be18 100644 --- a/openaerostruct/integration/aerostruct_groups.py +++ b/openaerostruct/integration/aerostruct_groups.py @@ -445,7 +445,7 @@ def setup(self): # column buckling for outer strut segment. Use strut-wing joint load as a compression load perf_group.add_subsystem( "column_buckling_outer", - EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_outer), + EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_outer, column_length_factor=surface["column_length_factor_outer"]), promotes_inputs=[("joint_load", "joint_load_strut_wing")], promotes_outputs=[("failure_column_buckling", "failure_column_buckling_outer")] ) @@ -455,7 +455,7 @@ def setup(self): # column buckling for inner strut segment. Use strut-root reaction force as a compression load perf_group.add_subsystem( "column_buckling_inner", - EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_inner), + EulerColumnBucklingFailureKS(surface=surface, joint_load_type="vector", ny=ny_inner, column_length_factor=surface["column_length_factor_inner"]), promotes_inputs=[("joint_load", "reaction_load_strut_root")], promotes_outputs=[("failure_column_buckling", "failure_column_buckling_inner")] ) @@ -466,7 +466,7 @@ def setup(self): ny = surface["mesh"].shape[1] perf_group.add_subsystem( "column_buckling", - EulerColumnBucklingFailureKS(surface=surface, ny=ny, joint_load_type="vector"), + EulerColumnBucklingFailureKS(surface=surface, ny=ny, joint_load_type="vector", column_length_factor=surface["column_length_factor"]), promotes_inputs=["nodes", "Iz", ("joint_load", "joint_load_strut_wing")], promotes_outputs=["failure_column_buckling"] ) diff --git a/openaerostruct/structures/failure_buckling_ks.py b/openaerostruct/structures/failure_buckling_ks.py index b2f8dcadf..c0d3688f2 100644 --- a/openaerostruct/structures/failure_buckling_ks.py +++ b/openaerostruct/structures/failure_buckling_ks.py @@ -163,6 +163,7 @@ def initialize(self): self.options.declare("ny", desc="number of nodes in the strut") self.options.declare("rho", types=float, default=100.0, desc="KS aggregation smoothness parameter") self.options.declare("joint_load_type", types=str, default="vector", desc="type of joint load, either `vector` or `scaler`. If scaler, this should be the axial tension force") + self.options.declare("column_length_factor", types=float, default=1.0, desc="factor to multiply the strut length to account for e.g. mount") def setup(self): ny = self.options["ny"] @@ -194,7 +195,7 @@ def compute(self, inputs, outputs): # strut length # The length we compute from the FEM nodes are longer than the actual length because e.g. we don't model fuselage so the strut root is at the symmetry plane. # To account for this, apply a factor < 1 to the strut length. - column_length_factor = self.options["surface"]["column_length_factor"] + column_length_factor = self.options["column_length_factor"] strut_len = np.linalg.norm(nodes[0, :] - nodes[-1, :]) * column_length_factor # compression load From fc1cb4f4f96c71142825f719034806b26a031f84 Mon Sep 17 00:00:00 2001 From: Shugo Kaneko <49300827+kanekosh@users.noreply.github.com> Date: Wed, 17 Sep 2025 23:42:06 -0400 Subject: [PATCH 20/30] Skip a failing test on OM 3.40 (#467) * skip test * set OM upper bound --- setup.py | 3 ++- tests/integration_tests/test_multi_single.py | 3 +++ 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/setup.py b/setup.py index b8e4b45e0..cb2b35cba 100644 --- a/setup.py +++ b/setup.py @@ -54,7 +54,7 @@ package_data={"openaerostruct": ["tests/*.py", "*/tests/*.py", "*/*/tests/*.py"]}, install_requires=[ # Remember to update the oldest versions in docs/installation.rst - "openmdao>=3.35", + "openmdao>=3.35,<=3.39", "numpy>=1.21", "scipy>=1.7", "matplotlib", @@ -68,3 +68,4 @@ plot_wingbox=openaerostruct.utils.plot_wingbox:disp_plot """, ) + diff --git a/tests/integration_tests/test_multi_single.py b/tests/integration_tests/test_multi_single.py index 3e2f805d5..f2cef9d20 100644 --- a/tests/integration_tests/test_multi_single.py +++ b/tests/integration_tests/test_multi_single.py @@ -1,7 +1,10 @@ from openmdao.utils.assert_utils import assert_near_equal import unittest +import openmdao +# skip this test on OpenMDAO 3.40.0 because it fails due to a bug in OM 3.40.0 +@unittest.skipIf(openmdao.__version__ == "3.40.0", "Skipped on OpenMDAO 3.40.0") class Test(unittest.TestCase): def test(self): import numpy as np From 33685eb564f66362d71bdf9c86b1b0a0f29db810 Mon Sep 17 00:00:00 2001 From: Shugo Kaneko <49300827+kanekosh@users.noreply.github.com> Date: Thu, 18 Sep 2025 10:26:02 -0400 Subject: [PATCH 21/30] Switch to Ruff linting and formatting (#466) * migrated local flake8 to ruff.toml * add extend to ruff local configu * whitespace fixes * trigger checks * Update ruff.toml configuration paths * Update gitignore * Automatic formatting changes * Non-automatic fixes * trigger build * Trigger Build * Formatting --------- Co-authored-by: Alasdair Gray --- .flake8 | 8 - .gitignore | 3 + .readthedocs.yaml | 3 +- CITATION.cff | 2 +- openaerostruct/aerodynamics/aero_groups.py | 2 +- .../aerodynamics/compressible_states.py | 1 + openaerostruct/aerodynamics/geometry.py | 2 +- openaerostruct/aerodynamics/lift_drag.py | 2 +- .../aerodynamics/mesh_point_forces.py | 1 + openaerostruct/aerodynamics/total_drag.py | 2 +- openaerostruct/aerodynamics/wave_drag.py | 4 +- .../docs/_n2html/generate_n2_n2.html | 54 +- openaerostruct/docs/_static/copybutton.js | 4 +- .../advanced_features/custom_mesh_example.rst | 2 +- .../customizing_prob_setup.rst | 2 +- .../docs/advanced_features/figs/CRM_esque.svg | 2 +- .../docs/advanced_features/figs/bspline.svg | 2 +- .../docs/advanced_features/figs/chord.svg | 2 +- .../docs/advanced_features/figs/dihedral.svg | 2 +- .../figs/inverted_gull-wing.svg | 2 +- .../advanced_features/figs/mesh-diagram.svg | 2 +- .../docs/advanced_features/figs/radius.svg | 2 +- .../docs/advanced_features/figs/sweep.svg | 2 +- .../docs/advanced_features/figs/taper.svg | 2 +- .../docs/advanced_features/figs/thickness.svg | 2 +- .../docs/advanced_features/figs/twist.svg | 2 +- .../docs/advanced_features/figs/xshear.svg | 2 +- .../docs/advanced_features/figs/zshear.svg | 2 +- .../docs/advanced_features/ground_effect.rst | 2 - .../docs/advanced_features/mphys_coupling.rst | 2 +- .../multi_section_surfaces.rst | 34 +- .../docs/advanced_features/multipoint.rst | 2 +- .../advanced_features/multipoint_parallel.rst | 8 +- .../openconcept_coupling.rst | 2 +- .../advanced_features/scripts/basic_2_sec.py | 5 +- .../scripts/basic_2_sec_AS.py | 3 +- .../scripts/basic_2_sec_construction.py | 5 +- .../scripts/basic_2_sec_visc.py | 5 +- .../advanced_features/scripts/run_vsp_777.py | 1 + openaerostruct/docs/aerostructural_index.rst | 2 +- openaerostruct/docs/index.rst | 2 +- openaerostruct/docs/user_reference.rst | 2 +- .../docs/user_reference/debugging_tips.rst | 2 +- .../docs/user_reference/mesh_surface_dict.rst | 96 +- .../docs/user_reference/v1_v2_conversion.rst | 2 +- openaerostruct/geometry/ffd_component.py | 3 +- openaerostruct/geometry/geometry_mesh.py | 2 +- .../geometry/geometry_mesh_transformations.py | 3 +- openaerostruct/geometry/radius_comp.py | 2 +- .../integration/aerostruct_groups.py | 2 +- .../meshing/section_mesh_generator.py | 2 +- openaerostruct/mphys/mux_surface_forces.py | 2 +- openaerostruct/structures/failure_ks.py | 2 +- openaerostruct/structures/fuel_loads.py | 4 - openaerostruct/structures/fuel_vol.py | 2 +- .../structures/non_intersecting_thickness.py | 2 +- .../structures/section_properties_tube.py | 2 +- .../structures/section_properties_wingbox.py | 20 +- openaerostruct/structures/spar_within_wing.py | 2 +- .../structures/spatial_beam_functionals.py | 3 +- .../structures/spatial_beam_setup.py | 3 +- openaerostruct/structures/structural_cg.py | 2 +- openaerostruct/structures/vonmises_tube.py | 2 +- openaerostruct/structures/weight.py | 2 +- .../structures/wing_weight_loads.py | 2 +- openaerostruct/utils/testing.py | 6 +- ruff.toml | 10 + setup.py | 1 - ...est_aerostruct_wingbox_+weight_analysis.py | 1 - .../test_aerostruct_wingbox_analysis.py | 1 - ...ostruct_wingbox_fuel_vol_constraint_opt.py | 887 ++++++++-------- .../test_aerostruct_wingbox_opt.py | 833 ++++++++-------- ...ct_wingbox_wave_fuel_vol_constraint_opt.py | 881 ++++++++-------- .../integration_tests/test_multi_2_sec_AS.py | 3 +- tests/integration_tests/test_multi_single.py | 1 + .../integration_tests/test_multi_single_AS.py | 3 +- .../test_multipoint_wingbox_aerostruct.py | 1 - .../test_multipoint_wingbox_derivs.py | 943 +++++++++--------- .../test_wingbox_analysis.py | 1 - .../integration_tests/test_wingbox_derivs.py | 921 +++++++++-------- .../test_wingbox_distributed_fuel.py | 1 - 81 files changed, 2423 insertions(+), 2428 deletions(-) delete mode 100644 .flake8 create mode 100644 ruff.toml diff --git a/.flake8 b/.flake8 deleted file mode 100644 index 7dfa6ed49..000000000 --- a/.flake8 +++ /dev/null @@ -1,8 +0,0 @@ -[flake8] -# OpenAeroStruct-specific settings in addition to the MDO Lab standard config file. -extend-exclude = - openaerostruct/docs/conf.py, - openaerostruct/docs/aero_walkthrough/part*.py, - openaerostruct/docs/_utils/docutil.py, - openaerostruct/docs/_utils/patch.py, - openaerostruct/utils/plot_*.py, diff --git a/.gitignore b/.gitignore index bf8c8aedf..8a939a13d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,6 +1,9 @@ # Ignore OpenMDAO reports folders **/reports/ +# pre-commit config file that should be manually copied from the mdolab/.github repo. +.pre-commit-config.yaml + *.pyc *.py~ *.db diff --git a/.readthedocs.yaml b/.readthedocs.yaml index 4d7f9738a..1257c6976 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -11,8 +11,7 @@ build: tools: python: "3.11" -sphinx: - configuration: openaerostruct/docs/conf.py + python: install: diff --git a/CITATION.cff b/CITATION.cff index f696d8c71..aec49a916 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -28,4 +28,4 @@ preferred-citation: number: "4" pages: "1815--1827" date: "2018-04" - doi: "10.1007/s00158-018-1912-8" \ No newline at end of file + doi: "10.1007/s00158-018-1912-8" diff --git a/openaerostruct/aerodynamics/aero_groups.py b/openaerostruct/aerodynamics/aero_groups.py index 0b5981d24..8e0ee2a1b 100644 --- a/openaerostruct/aerodynamics/aero_groups.py +++ b/openaerostruct/aerodynamics/aero_groups.py @@ -23,7 +23,7 @@ def initialize(self): "compressible", types=bool, default=False, - desc="Turns on compressibility correction for moderate Mach number " "flows. Defaults to False.", + desc="Turns on compressibility correction for moderate Mach number flows. Defaults to False.", ) def setup(self): diff --git a/openaerostruct/aerodynamics/compressible_states.py b/openaerostruct/aerodynamics/compressible_states.py index 5e9cf29e6..308cc052c 100644 --- a/openaerostruct/aerodynamics/compressible_states.py +++ b/openaerostruct/aerodynamics/compressible_states.py @@ -1,6 +1,7 @@ """ Class definition for CompressibleVLMStates. """ + import openmdao.api as om from openaerostruct.aerodynamics.get_vectors import GetVectors diff --git a/openaerostruct/aerodynamics/geometry.py b/openaerostruct/aerodynamics/geometry.py index 91ceeae1c..e93a18a22 100644 --- a/openaerostruct/aerodynamics/geometry.py +++ b/openaerostruct/aerodynamics/geometry.py @@ -11,7 +11,7 @@ class VLMGeometry(om.ExplicitComponent): Some of the quantities, like `normals`, are used to compute the RHS of the AIC linear system. - parameters + Parameters ---------- def_mesh[nx, ny, 3] : numpy array Array defining the nodal coordinates of the lifting surface. diff --git a/openaerostruct/aerodynamics/lift_drag.py b/openaerostruct/aerodynamics/lift_drag.py index a96615828..91fbae99c 100644 --- a/openaerostruct/aerodynamics/lift_drag.py +++ b/openaerostruct/aerodynamics/lift_drag.py @@ -8,7 +8,7 @@ class LiftDrag(om.ExplicitComponent): Calculate total lift and drag in force units based on section forces. This is for one given lifting surface. - parameters + Parameters ---------- sec_forces[nx-1, ny-1, 3] : numpy array Contains the sectional forces acting on each panel. diff --git a/openaerostruct/aerodynamics/mesh_point_forces.py b/openaerostruct/aerodynamics/mesh_point_forces.py index 4699d58c7..0188b25e4 100644 --- a/openaerostruct/aerodynamics/mesh_point_forces.py +++ b/openaerostruct/aerodynamics/mesh_point_forces.py @@ -1,6 +1,7 @@ """ Class definition for the MeshPointForces component. """ + import numpy as np import openmdao.api as om diff --git a/openaerostruct/aerodynamics/total_drag.py b/openaerostruct/aerodynamics/total_drag.py index 38fef9103..9b5e3ff63 100644 --- a/openaerostruct/aerodynamics/total_drag.py +++ b/openaerostruct/aerodynamics/total_drag.py @@ -4,7 +4,7 @@ class TotalDrag(om.ExplicitComponent): """Calculate total drag in force units. - parameters + Parameters ---------- CDi : float Induced coefficient of drag (CD) for the lifting surface. diff --git a/openaerostruct/aerodynamics/wave_drag.py b/openaerostruct/aerodynamics/wave_drag.py index 0b58e6b22..2b013d1f3 100644 --- a/openaerostruct/aerodynamics/wave_drag.py +++ b/openaerostruct/aerodynamics/wave_drag.py @@ -126,9 +126,7 @@ def compute_partials(self, inputs, partials): ccos = np.sum(widths * panel_mid_chords) ccos2w = np.sum(panel_mid_chords * widths**2 / lengths_spanwise) - davgdcos = ( - 2 * panel_mid_chords * widths / lengths_spanwise / ccos - panel_mid_chords * ccos2w / ccos**2 - ) + davgdcos = 2 * panel_mid_chords * widths / lengths_spanwise / ccos - panel_mid_chords * ccos2w / ccos**2 dtocdcos = ( panel_mid_chords * t_over_c / ccos - panel_mid_chords * np.sum(panel_mid_chords * widths * t_over_c) / ccos**2 diff --git a/openaerostruct/docs/_n2html/generate_n2_n2.html b/openaerostruct/docs/_n2html/generate_n2_n2.html index 11ce9d1b5..129601531 100644 --- a/openaerostruct/docs/_n2html/generate_n2_n2.html +++ b/openaerostruct/docs/_n2html/generate_n2_n2.html @@ -356,7 +356,7 @@ /* transition: opacity 0.25s; */ box-shadow: 3px 3px 3px 1px rgba(0, 0, 0, 0.2); border-radius: 10px; - + overflow: hidden; } @@ -487,7 +487,7 @@ .window-theme-value-info .window-body table tbody tr td { text-align: left; font-size: 9pt; - margin: 0; + margin: 0; border-top: 0; border-left: 0; border-bottom: 1px solid #a0a0a0; @@ -507,7 +507,7 @@ .show_value_button { background-color: #c1c1c1; } - + .copy_value_button { background-color: #c1c1c1; } @@ -2858,7 +2858,7 @@

Toolbar Help