Substructuring using the Transmission Simulator Method
######################################################
This example will demonstrate the capabilities in SDynPy for performing more
advanced substructuring calculations using the Transmission Simulator Method.
The transmission simulator method is a useful form of substructuring which
allows the interface of the substructures to be mass-loaded, which provides a
better set of basis shapes to use in the substructuring calculation.
.. contents::
Setting up the Problem
----------------------
We will first use the Cubit software to generate a mesh for the finite element
analysis that will be used in this example problem. We will create a 3-part
model of a box inside of a container. The container is labeled ``A``, the lid is
labeled ``B``, and the internal box is labeled ``C``. The attached input file can
be executed using Cubit's Python interface to create the model.
.. figure:: figures/transmission_simulator_setup.png
:width: 600
:alt: Mesh generated by Cubit
:align: center
:figclass: align-center
Mesh generated by Cubit shown both in whole and in cross section to reveal
the three components.
:download:`Cubit Input File<assets/create_models.py>`
The model contains 3 element blocks, block 1, block 2, and block 3, which are
the container ``A``, the lid ``B``, and the internal box ``C``, respectively. The attached
input file will generate several exodus files that can be used to investigate
substructuring. Each object will be generated independently, ``a.exo``, ``b.exo``,
and ``c.exo``. Additionally, meshes will be generated for pairs of components,
``ab.exo`` and ``bc.exo``. Finally, a full model ``abc.exo`` will be generated which
will provide truth data to compare the substructuring results to.
Input files are then set up for the Sierra/SD code. Each input file will
contain an Eigensolution ``eigen`` which will solve for the component modes
used in the substructuring. The input files for each mesh are attached:
:download:`Sierra/SD Input for ABC<assets/abc.inp>`
:download:`Sierra/SD Input for AB<assets/ab.inp>`
:download:`Sierra/SD Input for BC<assets/bc.inp>`
:download:`Sierra/SD Input for A<assets/a.inp>`
:download:`Sierra/SD Input for B<assets/b.inp>`
:download:`Sierra/SD Input for C<assets/c.inp>`
When run, these will generate output files containing the natural frequencies
and mode shapes of the structures.
To perform the transmission simulator approach, we will add substructure ``AB``
to substructure ``BC`` and then subtract off the extra ``B``. In this case, the
structure ``B`` serves as the transmission simulator because it mass-loads the
interfaces of ``A`` and ``C`` to give a better set of basis shapes to use in the
substructuring.
.. figure:: figures/transmission_simulator_math.png
:width: 600
:alt: Transmission Simulator Operations
:align: center
:figclass: align-center
Substructuring operations performed in the transmission simulator method.
The system ``AB`` is added to ``BC`` and then system ``B`` is subtracted to
produce system ``ABC``.
Loading the results into SDynPy
-------------------------------
We will now load the models into SDynPy in order to perform the analysis. We
will start by importing the required modules and setting up plotting options.
.. code-block:: python
# Import modules
import sdynpy as sdpy # For dynamics
import numpy as np # For numerics
import matplotlib.pyplot as plt # For plotting
from scipy.spatial import QhullError # For creating elements for visualization
# Set up options to use for plotting
plot_options = {'node_size':0,'line_width':1,'show_edges':False,
'view_up':[0,1,0]}
We will now load in the models and plot their mode shapes to make sure
everything looks right. For each model we will reduce to just the surface
nodes and elements for visualization using the
:py:func:`reduce_to_surfaces<sdynpy.fem.sdynpy_exodus.Exodus.reduce_to_surfaces>`
method of the :py:class:`Exodus<sdynpy.fem.sdynpy_exodus.Exodus>` class.
We will then create :py:class:`Geometry<sdynpy.core.sdynpy_geometry.Geometry>`
as well as a :py:class:`Shapes<sdynpy.core.sdynpy_shape.ShapeArray>` from the
finite element data using the respective
:py:func:`sdpy.geometry.from_exodus<sdynpy.core.sdynpy_geometry.from_exodus>` and
:py:func:`sdpy.shape.from_exodus<sdynpy.core.sdynpy_shape.from_exodus>` methods.
.. code-block:: python
# Specify the models to load into SDynPy
models = ['a','b','c','ab','bc','abc']
# For each model, load in the exodus file and reduce it to surfaces
fexos = {model:sdpy.Exodus(model+'-out.exo').reduce_to_surfaces()
for model in models}
# For each model create a SDynPy Geometry from the finite element model
geometries = {model:sdpy.Geometry.from_exodus(fexos[model])
for model in models}
# For each model create a SDynPy ShapeArray from the eigensolution results
shapes = {model:sdpy.shape.from_exodus(fexos[model])
for model in models}
# For each model, plot the mode shapes on the geometry
for model in models:
geometries[model].plot_shape(shapes[model],plot_options,
deformed_opacity=0.5,
undeformed_opacity=0)
.. figure:: figures/transmission_simulator_A_mode_example.gif
:width: 600
:alt: A Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``A``
.. figure:: figures/transmission_simulator_B_mode_example.gif
:width: 600
:alt: B Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``B``
.. figure:: figures/transmission_simulator_C_mode_example.gif
:width: 600
:alt: C Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``C``
.. figure:: figures/transmission_simulator_AB_mode_example.gif
:width: 600
:alt: AB Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``AB``
.. figure:: figures/transmission_simulator_BC_mode_example.gif
:width: 600
:alt: BC Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``BC``
.. figure:: figures/transmission_simulator_ABC_mode_example.gif
:width: 600
:alt: ABC Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``ABC``
Reducing to Test Degrees of Freedom
-----------------------------------
For this example, we will pretend that the finite element results are actually
test data, and for this reason, we will reduce to a more limited set of degrees
of freedom in each model that may have been measured in a test. We will first
define the geometry of the system so we can select the final kept set of nodes
by position.
.. code-block:: python
container_size = [1.0,1.2,1.5]
component_size = [0.4,0.4,0.7]
lid_depth = 0.15
container_thickness = 0.1
lid_position = container_size[-1]/2-lid_depth
component_position = container_size[-1]/2-container_thickness-component_size[-1]/2
grid_size = 3
component_grid_size = 2
grid_inset = 0.1
We will then loop through each face of the model to get sensor positions that
we would like to use in the test. We will use NumPy functions
`linspace <https://numpy.org/doc/stable/reference/generated/numpy.linspace.html>`_
and `meshgrid <https://numpy.org/doc/stable/reference/generated/numpy.meshgrid.html>`_.
to produce grids of sensors on each surface.
.. code-block:: python
# Loop through each face to get sensor positions
positions = []
indices = []
index = 0
# Outer box
for dimension in range(3):
other_dimensions = [v for v in range(3) if not v == dimension]
meshgrid_inputs = [None,None,None]
grid_0 = container_size[dimension]/2
grid_1 = np.linspace(-container_size[other_dimensions[0]]/2+grid_inset,
container_size[other_dimensions[0]]/2-grid_inset,
grid_size)
grid_2 = np.linspace(-container_size[other_dimensions[1]]/2+grid_inset,
container_size[other_dimensions[1]]/2-grid_inset,
grid_size)
meshgrid_inputs[dimension] = grid_0
meshgrid_inputs[other_dimensions[0]] = grid_1
meshgrid_inputs[other_dimensions[1]] = grid_2
this_positions = np.moveaxis(np.array(np.meshgrid(*meshgrid_inputs,indexing='ij')).squeeze(),0,-1).reshape(-1,3)
positions.append(this_positions)
indices.append(np.ones(this_positions.shape[0])*index)
index += 1
# Now do the negative side
this_positions_opposite = this_positions.copy()
this_positions_opposite[...,dimension] *= -1
positions.append(this_positions_opposite)
indices.append(np.ones(this_positions.shape[0])*index)
index += 1
# Component box
for dimension in range(2):
other_dimensions = [v for v in range(3) if not v == dimension]
meshgrid_inputs = [None,None,None]
grid_0 = component_size[dimension]/2
grid_1 = np.linspace(-component_size[other_dimensions[0]]/2+grid_inset,
component_size[other_dimensions[0]]/2-grid_inset,
component_grid_size)
grid_2 = np.linspace(-component_size[other_dimensions[1]]/2+grid_inset+component_position,
component_size[other_dimensions[1]]/2-grid_inset+component_position,
component_grid_size)
meshgrid_inputs[dimension] = grid_0
meshgrid_inputs[other_dimensions[0]] = grid_1
meshgrid_inputs[other_dimensions[1]] = grid_2
this_positions = np.moveaxis(np.array(np.meshgrid(*meshgrid_inputs,indexing='ij')).squeeze(),0,-1).reshape(-1,3)
positions.append(this_positions)
indices.append(np.ones(this_positions.shape[0])*index)
index += 1
# Now do the negative side
this_positions_opposite = this_positions.copy()
this_positions_opposite[...,dimension] *= -1
positions.append(this_positions_opposite)
indices.append(np.ones(this_positions.shape[0])*index)
index += 1
positions = np.concatenate(positions,axis=0)
indices = np.concatenate(indices,axis=0)
We will now reduce each of the geometries down to its test sensors using
the :py:func:`reduce<sdynpy.core.sdynpy_geometry.Geometry.reduce>` method,
using the :py:func:`by_position<sdynpy.core.sdynpy_geometry.NodeArray.by_position>`
method of the :py:class:`NodeArray<sdynpy.core.sdynpy_geometry.NodeArray>`
object to select which nodes to keep. Note that this will simply select the
closest node to the position specified, which isn't exactly what we want. For
example, if a position is specified that corresponds to a node on system ``A``,
but the system we are generating test data for is system ``BC`` (which doesn't
contain system ``A``) it will instead select a node the closest node on ``BC``
to that point, rather than simply not selecting a node at all. We will discard
these points and connect the test nodes with elements and tracelines to ease
visualization.
We will specify a distance threshold to compare the point where we wanted to
select a node and where the closest node was selected, and use that to discard
improperly selected nodes. We will also specify a maximum condition number for
triangles created as elements for visualization. We will also set up a
dictionary of :py:class:`id_map<sdynpy.core.sdynpy_geometry.id_map>` objects
to map the nodes between the test geometry and the finite element model.
.. code-block:: python
distance_threshold = 0.05
max_condition = 10
node_id_maps = {}
We will then loop through all of our test geometries and finish their creation.
.. code-block:: python
# Now create the test geometries
for model in test_geometries:
print('Creating test geometry for {:}'.format(model))
Within this loop, we will set up a node map from the finite element model nodes
to the test nodes, which will be labeled 1 through $n_nodes$ for each system,
where $n_nodes$ is the number of nodes for each test geometry.
.. code-block:: python
geometry = test_geometries[model]
from_ids = geometry.node.id
to_ids = np.arange(geometry.node.size)+1
node_id_maps[model] = sdpy.id_map(from_ids, to_ids)
At this point, we will now discard the nodes that aren't close enough to their
desired positions, and then use the node map to rename the nodes in the model.
.. code-block:: python
# Throw away nodes that aren't in the model
distance = np.linalg.norm(geometry.node.coordinate-positions,axis=-1)
keep_indices = distance < distance_threshold
new_geometry = geometry.reduce(geometry.node.id[keep_indices])
new_geometry.node.id = node_id_maps[model](new_geometry.node.id)
original_indices = indices[keep_indices]
We will then add elements using the
:py:func:`triangulate<sdynpy.core.sdynpy_geometry.NodeArray.triangulate>` method
of the :py:class:`NodeArray<sdynpy.core.sdynpy_geometry.NodeArray>` object. If
the nodes cannot be triangulated, they are instead connected with a traceline.
.. code-block:: python
# Create elements
elements = []
for value in np.unique(original_indices):
face_indices = original_indices == value
face_nodes = new_geometry.node[face_indices]
try:
face_elements = face_nodes.triangulate(new_geometry.coordinate_system,
condition_threshold=max_condition)
elements.append(np.atleast_1d(face_elements))
except QhullError:
new_geometry.add_traceline(face_nodes.id)
new_geometry.element = np.concatenate(elements)
new_geometry.element.id = np.arange(new_geometry.element.size)+1
test_geometries[model] = new_geometry
We can then similarly map the shapes to the new test geometries. We will use
the
:py:func:`transform_coordinate_system<sdynpy.core.sdynpy_shape.ShapeArray.transform_coordinate_system>`
method of the :py:class:`ShapeArray<sdynpy.core.sdynpy_shape.ShapeArray>` object
to simultaneously reduce the shape down to the test geometry and apply the
node map. We will also assign a 1% damping to each of the test shapes. We will
plot the reduced shapes as well as a MAC matrix to ensure the shapes are still
distinguishable in their reduced form.
.. code-block:: python
# Now create the test shapes
test_shapes = {}
for model in shapes:
print('Creating test shapes for {:}'.format(model))
shape = shapes[model]
original_geometry = geometries[model]
new_geometry = test_geometries[model]
id_map = node_id_maps[model]
test_shape = shape.transform_coordinate_system(original_geometry,new_geometry,id_map)
test_shape.damping = 0.01
test_shapes[model] = test_shape
test_geometries[model].plot_shape(test_shapes[model],plot_options,undeformed_opacity=0,deformed_opacity=0.5)
sdpy.correlation.matrix_plot(sdpy.shape.mac(test_shapes[model]))
.. figure:: figures/transmission_simulator_A_mode_reduced_example.gif
:width: 600
:alt: Test A Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``A``
.. figure:: figures/transmission_simulator_B_mode_reduced_example.gif
:width: 600
:alt: Test B Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``B``
.. figure:: figures/transmission_simulator_C_mode_reduced_example.gif
:width: 600
:alt: Test C Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``C``
.. figure:: figures/transmission_simulator_AB_mode_reduced_example.gif
:width: 600
:alt: Test AB Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``AB``
.. figure:: figures/transmission_simulator_BC_mode_reduced_example.gif
:width: 600
:alt: Test BC Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``BC``
.. figure:: figures/transmission_simulator_ABC_mode_reduced_example.gif
:width: 600
:alt: Test ABC Mode Shape
:align: center
:figclass: align-center
Test Mode Shape for System ``ABC``
Substructuring using the Transmission Simulator Method
------------------------------------------------------
Now that we have our models set up, we can proceed with the substructuring.
We will define our bandwidth as 2000 Hz, and select to use 15 transmission
simulator modes to perform the substructuring. We can then reduce down to
just the shapes satisfying these critera, and create modal
:py:class:`System<sdynpy.core.sdynpy_system.System>` from them. This will
result in diagonal mass, stiffness, and damping matrices, with the mode shape
matrix as the transformation to physical coordinates.
Recall, we will be attempting to compute system ``ABC`` from components ``AB``,
``BC``, and ``B``.
.. code-block:: python
bandwidth = 2000
transmission_simulator_modes = 15
system_ab = test_shapes['ab'][test_shapes['ab'].frequency < bandwidth].system()
system_bc = test_shapes['bc'][test_shapes['bc'].frequency < bandwidth].system()
system_b = test_shapes['b'][:transmission_simulator_modes].system()
The first step is to combine the components ``AB``, ``BC``, and ``B`` into one
:py:class:`System<sdynpy.core.sdynpy_system.System>` object. Note that no
constraints have been applied, so the system matrices will be block diagonal.
We will also combine the geometries so we can eventually plot mode shapes of
the combined system. We will combine the geometries first, as it will give us
a node offset that we can use to avoid conflicting nodes in the various
substructures using the
:py:func:`overlay_geometries<sdynpy.core.sdynpy_geometry.Geometry.overlay_geometries>`
static method of the :py:class:`Geometry<sdynpy.core.sdynpy_geometry.Geometry>`
object. We pass it a list of geometries to combine, a set of colors to make
each geometry, as well as tell it to return the node offset used to avoid
conflicts. We will then combine the systems using the
:py:func:`concatenate<sdynpy.core.sdynpy_system.System.concatenate>` static
method of the :py:class:`System<sdynpy.core.sdynpy_system.System>` object
to combine the systems, passing the node offset as an argument to keep the
concatenated systems consistent with the combined geometries. Note that because
we will be subtracting component ``B`` from the assembly, we will specify a
negative ``system_b`` during the concatenation. We can use the
:py:func:`System.spy<sdynpy.core.sdynpy_system.System.spy>` method to visualize
the structure of the created system and ensure that it is block-diagonal.
.. code-block:: python
geoms = (test_geometries['ab'],test_geometries['bc'],test_geometries['b'])
combined_geometry,node_offset = sdpy.Geometry.overlay_geometries(
geoms,
color_override=[1,7,11],
return_node_id_offset=True)
systems = (system_ab,system_bc,-system_b) # Note that system_b is negative
combined_system = sdpy.System.concatenate(
systems,node_offset)
combined_system.spy()
.. figure:: figures/transmission_simulator_concatenated_systems.png
:width: 600
:alt: Concatenated System Structure
:align: center
:figclass: align-center
Structure of the concatenated system, showing block-diagonal structure.
If we were to plot this concatenated system's modes, it would be like we just
overlaid all three systems on top of one another. This is because no constraints
have yet been applied, so, e.g., component ``AB`` has no effect on component
``BC``.
.. code-block:: python
combined_geometry.plot_shape(combined_system.eigensolution(),plot_options)
Note however, that this will not work in its current state, due to the negative
mass on component ``B``. If we instead flip ``B`` to positive and re-run the
previous code block, we can then plot the shapes. If you do this, don't forget
to flip the sign back to negative, otherwise the substructuring won't work!
.. figure:: figures/transmission_simulator_concatenated_B_mode.gif
:width: 600
:alt: Concatenated ``B`` Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``B`` in the concatenated system
.. figure:: figures/transmission_simulator_concatenated_AB_mode.gif
:width: 600
:alt: Concatenated ``AB`` Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``AB`` in the concatenated system
.. figure:: figures/transmission_simulator_concatenated_BC_mode.gif
:width: 600
:alt: Concatenated ``BC`` Mode Shape
:align: center
:figclass: align-center
Mode Shape for System ``BC`` in the concatenated system
Finally, we can apply constraints to the model. We will first find the nodes
that are used in the constraint. These will be the nodes in ``B`` that are also
in ``AB``. We will create another set of
:py:class:`id_map<sdynpy.core.sdynpy_geometry.id_map>` objects to map the
degrees of freedom from the original test geometries to the combined test
geometry, which has the node offset applied.
.. code-block:: python
connection_nodes = np.intersect1d(test_geometries['ab'].node.id,test_geometries['b'].node.id)
connection_map_ab = sdpy.id_map(connection_nodes,connection_nodes + node_offset)
connection_map_bc = sdpy.id_map(connection_nodes,connection_nodes + node_offset*2)
connection_map_b = sdpy.id_map(connection_nodes,connection_nodes + node_offset*3)
Because the transmission simulator method uses softened, least-squares constraints,
we must also construct the basis set of shapes used in the constraint. These
will be the mode shapes of the transmission simulator mapped to the nodes in the
combined system for each substructure. To help verify that they are assembled
correctly, these shapes can be plotted on the combined geometry. If degrees of
freedom have been correctly selected for each constraint, the boundary nodes
of the substructures should be seen to move identically in each shape.
.. code-block:: python
connection_shapes_ab = test_shapes['b'][:transmission_simulator_modes].copy()
connection_shapes_ab.coordinate.node = connection_map_ab(connection_shapes_ab.coordinate.node)
connection_shapes_bc = test_shapes['b'][:transmission_simulator_modes].copy()
connection_shapes_bc.coordinate.node = connection_map_bc(connection_shapes_bc.coordinate.node)
connection_shapes_b = test_shapes['b'][:transmission_simulator_modes].copy()
connection_shapes_b.coordinate.node = connection_map_b(connection_shapes_b.coordinate.node)
# Combine all the shape degrees of freedom into one shape
connection_shapes = sdpy.shape.concatenate_dofs([connection_shapes_ab,connection_shapes_bc,connection_shapes_b])
# Now to make sure this is correct, if we plot these shapes on the combined
# geometry, it should move all of the connection degrees of freedom in the
# shapes of the transmission simulator
combined_geometry.plot_shape(connection_shapes,plot_options,deformed_opacity=0.5,
undeformed_opacity=0)
Once we have these shapes, we can select the coordinates in each system
in each shape and then create the constraint matrix using the
:py:func:`substructure_by_shape<sdynpy.core.sdynpy_system.System.substructure_by_shape>`
method of the :py:class:`System<sdynpy.core.sdynpy_system.System>` object.
To this method get passed the shapes to use as the basis for the constraints
as well as the degrees of freedom to use in each system for the computation of
the least-squares constraint. This method could directly return a constrained
system, but instead we set the argument ``return_constrained_system=False``,
which instead makes the method return the constraint matrix instead. We
do this because we will be applying another set of shape constraints for the
other substructure, and wish to apply the constraints all at once.
.. code-block:: python
connection_dofs_ab = combined_system.coordinate[np.in1d(combined_system.coordinate.node,
connection_map_ab.to_ids)]
connection_dofs_bc = combined_system.coordinate[np.in1d(combined_system.coordinate.node,
connection_map_bc.to_ids)]
connection_dofs_b = combined_system.coordinate[np.in1d(combined_system.coordinate.node,
connection_map_b.to_ids)]
constraint_matrix_ab_b = combined_system.substructure_by_shape(connection_shapes,
connection_dofs_ab,
connection_dofs_b,
return_constrained_system=False)
constraint_matrix_bc_b = combined_system.substructure_by_shape(connection_shapes,
connection_dofs_bc,
connection_dofs_b,
return_constrained_system=False)
Note that all of the bookkeeping is handled automatically by the
:py:func:`substructure_by_shape<sdynpy.core.sdynpy_system.System.substructure_by_shape>`
method, we simply had to specify which degrees of freedom were to be used in the
substructuring, and ensure that those degrees of freedom existed in the shapes
used for the modal constraints.
We can then concatenate the two constraint matrices into a single constraint
matrix and apply it to the system using the
:py:func:`constrain<sdynpy.core.sdynpy_system.System.constrain>`
method of the :py:class:`System<sdynpy.core.sdynpy_system.System>` object.
.. code-block:: python
constraint_matrix = np.concatenate((constraint_matrix_ab_b,constraint_matrix_bc_b),axis=0)
constrained_system = combined_system.constrain(constraint_matrix)
We can then solve for the modes of the constrained system, and compare them
against the truth shapes from the full finite element model.
.. code-block:: python
constrained_shapes = constrained_system.eigensolution()
combined_geometry.plot_shape(constrained_shapes,plot_options,deformed_opacity=0.5,
undeformed_opacity=0.0)
test_geometries['abc'].plot_shape(test_shapes['abc'],plot_options,deformed_opacity=0.5,
undeformed_opacity=0.0)
.. figure:: figures/transmission_simulator_truth_mode_shape.gif
:width: 600
:alt: Example truth mode shape
:align: center
:figclass: align-center
Example truth mode to compare to the substructured modes
.. figure:: figures/transmission_simulator_substructured_mode_shape.gif
:width: 600
:alt: Example substructured mode shape
:align: center
:figclass: align-center
Example substructured mode showing the the constraints have been satisfied.
Note again that the :py:func:`constrain<sdynpy.core.sdynpy_system.System.constrain>`
method automatically updated not only the system mass, stiffness, and damping
matrices to apply the constraint, but also updated the transformation from the
constrained space to the physical degrees of freedom. Therefore when we
compute the eigensolution, the shapes are already set up in the physical degrees
of freedom.
Summary
-------
This example has demonstrated how substructuring using the transmission
simulator method can be performed. A set of finite element models was created
that serve as the example problem. These were then reduced to a simulated test
geometry that was used to perform the substructuring.
With the models set up, the components were reduced to a test bandwidth and
concatenated to form a single system. A set of basis shapes were constructed
to perform the least-squares modal constraints. These were then passed to the
substructuring calculation along with the degrees of freedom to use for each
constraint. The constraint matrices were then used to constrain the system, and
plotting the mode shapes revealed that the model was constrained successfully.