Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adding example for TACS-pyGeo coupling #103

Merged
merged 4 commits into from
Jun 10, 2022
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3,255 changes: 3,255 additions & 0 deletions examples/struct_only/beam_shape_opt/Slender_Beam.bdf

Large diffs are not rendered by default.

14 changes: 14 additions & 0 deletions examples/struct_only/beam_shape_opt/ffd_8_linear.fmt
Original file line number Diff line number Diff line change
@@ -0,0 +1,14 @@
1
8 2 2
-0.001 0.1421428571428571 0.2852857142857142 0.42842857142857127 0.5715714285714284 0.7147142857142855 0.8578571428571425 1.001
-0.001 0.1421428571428571 0.2852857142857142 0.42842857142857127 0.5715714285714284 0.7147142857142855 0.8578571428571425 1.001
-0.001 0.1421428571428571 0.2852857142857142 0.42842857142857127 0.5715714285714284 0.7147142857142855 0.8578571428571425 1.001
-0.001 0.1421428571428571 0.2852857142857142 0.42842857142857127 0.5715714285714284 0.7147142857142855 0.8578571428571425 1.001
-0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004
0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004
-0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004 -0.051000000000000004
0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004 0.051000000000000004
-0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05
-0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05 -0.05
0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
0.05 0.05 0.05 0.05 0.05 0.05 0.05 0.05
140 changes: 140 additions & 0 deletions examples/struct_only/beam_shape_opt/shape_opt.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,140 @@
"""
This example demonstrates TACS structural optimization capabilities.
The beam model that we will be using for this problem is a rectangular beam,
cantilevered, with a shear load applied at the tip. The beam is discretized using
1001 shell elements along it's span and depth.
The optimization problem is as follows:
Minimize the mass of the beam with respect to the depth of the cross-section along the span,
subject to a max stress constraint dictated by the materials yield stress.
In order to change the shape of the FEM we use a free-form deformation (FFD) volume
parmeterization scheme provided by the pyGeo library.
An aproximate analytical solution can be derived from beam theory,
by realizing that the stress at any spanwise cross-section in the beam
can be found independently using:
sigma(x,y) = y*M(x)/I
An analytical solution for this problem can be shown to be:
d(x) = sqrt(6*V*(L-x)/(t*sigma_y))
The optimization is setup using TACS' MPHYS module, which acts as a wrapper
for OpenMDAO.
"""

import numpy as np
import matplotlib.pyplot as plt
import os

import openmdao.api as om
from pygeo.mphys import OM_DVGEOCOMP

from mphys import Multipoint
from mphys.scenario_structural import ScenarioStructural
from tacs.mphys import TacsBuilder
from tacs import elements, constitutive, functions

bdf_file = os.path.join(os.path.dirname(__file__), 'Slender_Beam.bdf')
ffd_file = os.path.join(os.path.dirname(__file__), 'ffd_8_linear.fmt')

# Beam thickness
t = 0.01 # m
# Length of beam
L = 1.0

# Material properties
rho = 2780.0 # kg /m^3
E = 70.0e9
nu = 0.0
ys = 420.0e6

# Shear force applied at tip
V = 2.5E4

# Callback function used to setup TACS element objects and DVs
def element_callback(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwargs):
# Setup (isotropic) property and constitutive objects
prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys)
con = constitutive.IsoShellConstitutive(prop, t=t, tNum=-1)
# TACS shells are sometimes a little overly-rigid in shear
# We can reduce this effect by decreasing the drilling regularization
con.setDrillingRegularization(0.1)
refAxis = np.array([1.0, 0.0, 0.0])
transform = elements.ShellRefAxisTransform(refAxis)
elem = elements.Quad4Shell(transform, con)
return elem

def problem_setup(scenario_name, fea_assembler, problem):
"""
Helper function to add fixed forces and eval functions
to structural problems used in tacs builder
"""

# Add TACS Functions
problem.addFunction('mass', functions.StructuralMass)
problem.addFunction('ks_vmfailure', functions.KSFailure, safetyFactor=1.0,
ksWeight=100.0)

# Add forces to static problem
problem.addLoadToNodes(1112, [0.0, V, 0.0, 0.0, 0.0, 0.0], nastranOrdering=True)


class Top(Multipoint):
def setup(self):
tacs_options = {'element_callback': element_callback,
'problem_setup': problem_setup,
'mesh_file': bdf_file}

# Initialize MPHYS builder for TACS
struct_builder = TacsBuilder(tacs_options, coupled=False)
struct_builder.initialize(self.comm)

# Add mesh component
self.add_subsystem('mesh', struct_builder.get_mesh_coordinate_subsystem())

# add the geometry component, we dont need a builder because we do it here.
self.add_subsystem("geometry", OM_DVGEOCOMP(ffd_file=ffd_file))
self.geometry.nom_add_discipline_coords("struct")

self.mphys_add_scenario('tip_shear', ScenarioStructural(struct_builder=struct_builder))

self.connect("mesh.x_struct0", "geometry.x_struct_in")
self.connect("geometry.x_struct0", "tip_shear.x_struct0")

def configure(self):
# Create reference axis
nRefAxPts = self.geometry.nom_addRefAxis(name="centerline", alignIndex='i', yFraction=0.5)

# Set up global design variables
def depth(val, geo):
for i in range(nRefAxPts):
geo.scale_y["centerline"].coef[i] = val[i]

self.geometry.nom_addGeoDVGlobal(dvName="depth", value=np.ones(nRefAxPts), func=depth)


################################################################################
# OpenMDAO setup
################################################################################
# Instantiate OpenMDAO problem
prob = om.Problem()
prob.model = Top()
model = prob.model

# Declare design variables, objective, and constraint
model.add_design_var('geometry.depth', lower=1e-3, upper=10.0, scaler=20.0)
model.add_objective('tip_shear.mass', scaler=1.0)
model.add_constraint('tip_shear.ks_vmfailure', lower=0.0, upper=1.0, scaler=1.0)

# Configure optimizer
prob.driver = om.ScipyOptimizeDriver(debug_print=['objs', 'nl_cons'], maxiter=1000)
prob.driver.options['optimizer'] = 'SLSQP'

# Setup OpenMDAO problem
prob.setup()

# Output N2 representation of OpenMDAO model
om.n2(prob, show_browser=False, outfile='beam_opt_n2.html')

# Run optimization
prob.run_driver()