Substructuring using Craig-Bampton Superelements from Sierra/SD

This example will demonstrate the capabilities in SDynPy for working with superelements from Sierra/SD. We will generate a small toy problem that we can use to demonstrate the functionality. We will then run Sierra to generate output, load it into SDynPy, and combine the components together.

Setting up the Problem

We will first use the Cubit software to generate a mesh for the finite element analysis. We will create a small spaceship-looking model consisting of 4 rectangular sections. The attached input file can be executed using Cubit’s Python interface to create the model.

Mesh generated by Cubit

Mesh generated by Cubit showing the two element block components and the boundary nodes between them.

Cubit Input File

This model contains two element blocks, block 1 and block 2, which are the ‘nose’ and the ‘thrusters’ of the spaceship, respectively. The model also contains a single nodeset 1 which specifies the boundary nodes that will be retained in the Craig Bampton Substructuring. The attached Cubit input file will generate three exodus files. One containing element block 1 block_1.exo, a second containing element block 2 block_2.exo, and a third containing both element blocks 1 and 2 block_all.exo, which is to be used for truth data.

Input files are then set up for the Sierra/SD code. In each input file, there will be an Eigensolution eigen and a Craig-Bampton Reduction cbr solution. The cbr solution will produce the superelement that will later be used for substructuring, and the eigen solution will compute the modes of the full model, which can be compared to the reduced cbr model to check convergence. A ‘truth’ input file will also be generated to run the full model; this will only contain an eigen solution as no substructuring will be performed with this model. The input files are attached:

Sierra/SD Input for Block 1

Sierra/SD Input for Block 2

Sierra/SD Input for Truth Data

When run, each of the two cbr analyses will output a superelement file block_<i>-cbr.nc4 as well as an Exodus file block_<i>-cbr.exo. The superelement .nc4 file will contain the reduced mass and stiffness matrices (damping also, though this is zero) as well as a map denoting which reduction degrees of freedom correspond to physical degres of freedom and which correspond to fixed base modes. The Exodus file contains the Craig-Bampton transformation shapes, meaning its shapes are the fixed interface mode shapes and the constraint mode shapes. The shapes in this file transform the reduced degrees of freedom back to the physical degrees of freedom. Each analysis will also output an Exodus file called block_<i>-eig.exo. This will contain the eigensolution computed from the full mass and stiffness matrices, which can be used to check convergence of the reduce model.

The ‘truth’ analysis in block_all.inp will output only an Exodus file block_all-eig.exo which contains the eigensolution computed from the full combined system models, which will be treated as truth data for this work.

Analyzing Superelements in SDynPy

With the data generated, we can now load into SDynPy to perform analysis. We will basically load each superelement .nc4 file and transformation .exo file in as a System object. This will populate the mass, stiffness, and damping matrices with the values from the reduced Craig-Bampton model, and will populate the transformation with data from the .exo transformation file. This allows the model to appear to be the full finite element mesh while only relying on a much smaller number of degrees of freedom to form that representation.

As the first step in our analysis, we will import the required packages and define the model numbers we are loading. We also set up lists to hold the different properties of each model that we will load.

import sdynpy as sdpy
import numpy as np

# Get the block numbers
models = [1,2]

# Set up arrays
systems = []
geometries = []
boundary_dofs = []
check_shapes = []
reduced_shapes = []

We will then set up a for loop over each of our element blocks. For each block, we will create a system, geometry, and set of boundary degrees of freedom using the System.from_exodus_superelement method. This will return a System object representing the reduced system, a Geometry object that can be used to plot shapes or displacement from the reduced system, and finally a CoordinateArray object that contains the boundary degrees of freedom, to which forces or substructuring connections could be applied. We will overwrite the color of the geometry to ensure the two different portions are different colors.

Once these results are loaded, we can compute the eigensolution of the reduced component model. The goal here is to compare the accuracy of the reduced component model to the full solution, so we will load in the component eigensolution results and extract the modes using the ShapeArray.from_exodus method. We can assign comments to each set of shapes to ensure when they are plotted we can keep track of them.

We wish to overlay the shapes computed from the reduced and full systems, so we will compute the alignment between the shapes. This will return either a 1 or a -1 depending on whether or not the shape needs to be flipped to be in phase with the other shape. We can easily compute the alignment using the ShapeArray.shape_alignment static method. We can then use this alignment array to flip the sign on the appropriate reduced shape to make them in phase with the full shape.

Once the shapes are aligned, we can overlay them using the ShapeArray.overlay_shapes static method. This will return a combined Geometry object and a combined ShapeArray object consisting of the overlaid shapes. We can plot them using the Geometry.plot_shape method, specifying undeformed_opacity = 0 and deformed_opacity = 0.5 so we can see any differences between the overlapping shapes.

We also compute and plot a MAC matrix between the sets of shapes to identify the shape similarity. Note that due to the symmetry of the block_1 portion of the model, there are repeated eigenvalues, which will result in mode shapes that can be linear combinations of one another. For this reason there may not be perfect shape agreement per the MAC for this component even though they are technically correct.

At the end of the for loop we store all the pertinent data to the arrays we have set up.

# Loop through each model
for i in models:
    # Load in the super element
    system,geometry,dofs = sdpy.System.from_exodus_superelement('block_{:}-cbr.nc4'.format(i),
                                                                'block_{:}-cbr.exo'.format(i))
    # Set the color of the model so they aren't all blue
    geometry.element.color=1+(i-1)*6

    # Compute modes of the reduced system
    cb_shapes = system.eigensolution()

    # Load in the eigensolution to compare against
    full_shapes = sdpy.shape.from_exodus(sdpy.Exodus('block_{:}-eig.exo'.format(i)))

    # Assign comments so we know which is which
    cb_shapes.comment1 = 'Reduced Block {:}'.format(i)
    full_shapes.comment1 = 'Full Block {:}'.format(i)

    # Compute shape alignment and flip shapes if necessary
    alignment = sdpy.shape.shape_alignment(cb_shapes[:full_shapes.size],full_shapes[:full_shapes.size])
    cb_shapes[:full_shapes.size] = cb_shapes[:full_shapes.size]*alignment[:,np.newaxis]

    # Overlay shapes
    compare_geo, compare_shapes = sdpy.shape.overlay_shapes(
        [geometry,geometry],[cb_shapes[:full_shapes.size],full_shapes],color_override=[1,7])
    compare_shapes.comment1 = 'Blue: CBR, Green: Full'
    compare_geo.plot_shape(compare_shapes,undeformed_opacity=0,deformed_opacity=0.5)

    # Compute a MAC
    mac = sdpy.shape.mac(cb_shapes,full_shapes)
    ax = sdpy.matrix_plot(mac,text_size=6)
    ax.set_ylabel('Reduced Shapes')
    ax.set_xlabel('Full Shapes')
    ax.set_title('Block {:}'.format(i))

    # Append to arrays
    systems.append(system)
    geometries.append(geometry)
    boundary_dofs.append(dofs)
    check_shapes.append(full_shapes)
    reduced_shapes.append(cb_shapes)

We see that the modes show almost perfect agreement at low frequency, and it is only at the very high frequencies that the agreement tends to break down. Note that the rigid body modes and other repeated roots may show poor shape correlation due to the fact that any linear combination of the repeated modes’ shapes are also valid mode shapes.

Error

Block 1 Component

Block 2 Component

Nearly Zero

Block 1 component with nearly zero error Block 2 component with nearly zero error

Small

Block 1 component with small error Block 2 component with small error

Large

Block 1 component with large error Block 2 component with large error
MAC matrix for block 1

MAC matrix between shapes computed from the reduced model and shapes computed from the full model for block 1

MAC matrix for block 2

MAC matrix between shapes computed from the reduced model and shapes computed from the full model for block 2

Combining Component Models

With the systems loaded, we can now combine the models into a single assembly. The model was created such that the boundary nodes are collocated, so we can use the sdpy.system.substructure_by_position function to combine the nodes. The function takes a list of systems to combine and a list of geometries to use to identify collocated nodes. After we perform the substructuring, we can compute the eigensolution and plot the assembled shapes.

# Now let's do some substructuring
combined_system, combined_geometry = sdpy.system.substructure_by_position(systems, geometries)
combined_shapes = combined_system.eigensolution()
combined_geometry.plot_shape(combined_shapes)
Assembled mode shape

First elastic modes of the assembled system.

We can see by the mode shape that the system was assembled successfully as all the nodes at the boundary surface stay together.

We will now compare the substructured modes to the truth modes from the full system eigensolution. We can load in the Exodus file and extract the geometry and shapes.

# Load in the data to compare to
exo = sdpy.Exodus('block_all-eig.exo')
full_geometry = sdpy.geometry.from_exodus(exo)
full_shapes = sdpy.shape.from_exodus(exo)

Note however that we cannot directly compare the shapes. In the sdpy.system.substructure_by_position function, the node IDs are offset so there are no conflicts between substructures that may have the same node numbers. We therefore need to figure out the mapping between the substructured node IDs and the truth node IDs.

Since the nodes coordinates are identical between the substructured and truth models, we can use the node positions to determine the mapping. We will let the original nodes be the nodes from the truth model, and we will select the substructured nodes that correspond to the original nodes using the by_position method of the NodeArray object. We can then construct an id_map object using these nodes, which many correlation functions in SDynPy accept to map nodes that have differing IDs.

# Map the shapes back to the original nodes because the substructuring has
# automatically renamed them to avoid conflicts
original_nodes = full_geometry.node.id
# Get the closest node to the original node
substructured_nodes = combined_geometry.node.by_position(full_geometry.node.coordinate).id
# Create an ID map
id_map = sdpy.id_map(original_nodes,substructured_nodes)

Now that the node map between the substructured model and truth model has been defined, we can use it to perform correlations. We will compute the MAC matrix between the substructured shapes and the truth shapes to identify shape correlations, and then we will compute the frequency differences and print it to a table. Note that the node map is passed as an argument to the mac function, which informs the functions which nodes to compare between the two sets of shapes.

# Compute MAC matrix
mac = sdpy.shape.mac(combined_shapes,full_shapes,id_map)
ax = sdpy.matrix_plot(mac,text_size=6)
ax.set_ylabel('Substructured Shapes')
ax.set_xlabel('Truth Shapes')
ax.set_title('Full Model MAC')

# Print mode table
substructure_indices,truth_indices = np.where(mac[6:,6:] > 0.7)
print(sdpy.shape.shape_comparison_table(combined_shapes[substructure_indices+6],
                                        full_shapes[truth_indices+6],
                                        damping_format = None,node_id_map=id_map,
                                        table_format='markdown'))
MAC matrix between the truth shapes and substructured shapes

MAC matrix between the truth shapes and the substructured shapes.

|  Mode|  Freq 1 (Hz)|  Freq 2 (Hz)|  Freq Error|  MAC|
|------|-------------|-------------|------------|-----|
|     1|       463.46|       463.46|        0.0%|  100|
|     2|       585.96|       585.95|        0.0%|  100|
|     3|       750.96|       750.96|        0.0%|  100|
|     4|       834.17|       834.17|        0.0%|  100|
|     5|      1368.52|      1368.49|        0.0%|  100|
|     6|      1434.50|      1434.50|        0.0%|  100|
|     7|      1458.51|      1458.50|        0.0%|  100|
|     8|      1665.30|      1665.14|        0.0%|  100|
|     9|      1727.18|      1727.17|        0.0%|  100|
|    10|      2021.07|      2021.03|        0.0%|  100|
|    11|      2767.00|      2766.90|        0.0%|  100|
|    12|      3140.80|      3138.37|        0.1%|  100|
|    13|      3432.48|      3432.29|        0.0%|  100|
|    14|      3767.55|      3765.91|        0.0%|  100|
|    15|      3871.50|      3870.86|        0.0%|  100|
|    16|      4063.53|      4062.56|        0.0%|  100|
|    17|      4251.80|      4250.85|        0.0%|  100|
|    18|      4312.29|      4312.22|        0.0%|  100|
|    19|      4319.99|      4319.28|        0.0%|  100|
|    20|      4948.01|      4947.56|        0.0%|  100|
|    21|      5195.45|      5193.05|        0.0%|  100|
|    22|      5713.47|      5687.49|        0.5%|  100|
|    23|      6213.86|      6198.35|        0.3%|  100|
|    24|      6215.22|      6210.73|        0.1%|  100|
|    25|      6945.63|      6945.25|        0.0%|  100|
|    26|      7110.41|      7102.79|        0.1%|  100|
|    27|      7247.83|      7246.99|        0.0%|  100|
|    28|      7398.35|      7397.53|        0.0%|  100|
|    29|      7732.32|      7729.57|        0.0%|  100|
|    30|      8052.87|      8049.50|        0.0%|  100|
|    31|      8101.66|      7992.09|        1.4%|   97|
|    32|      8336.22|      8310.52|        0.3%|   99|
|    33|      8877.21|      8860.38|        0.2%|   99|
|    34|      9108.49|      9069.66|        0.4%|   99|
|    35|      9323.28|      9263.32|        0.6%|   99|
|    36|      9627.15|      9591.93|        0.4%|   98|
|    37|      9914.18|      9875.70|        0.4%|   99|
|    38|     10095.29|     10060.16|        0.3%|   98|
|    39|     10308.96|     10307.29|        0.0%|  100|
|    40|     10391.74|     10369.17|        0.2%|   99|
|    41|     12108.51|     12079.93|        0.2%|   99|

The results are very good. MAC values between the truth and substructure shapes are generally above 95, and only one mode has a frequency error of more than 1%, with several having a frequency error of less than 0.1%.

Finally, we can overlay the substructure shapes. We again compute alignments to ensure the shapes are not overlaid out of phase with one another. We can pass the node map to the shape_alignment static method to inform which nodes should be compared between the two sets of shapes. The shapes are then multiplied by the alignment vector to flip their signs if necessary, and finally the shapes can be overlaid and plotted, again using the ShapeArray.overlay_shapes static method. It isn’t until the 31st mode shape that there is a visible difference in the mode shape, so the substructuring can be considered very accurate.

# Align shapes
alignment = sdpy.shape.shape_alignment(combined_shapes[substructure_indices+6],full_shapes[truth_indices+6],id_map)
combined_shapes[substructure_indices+6] = combined_shapes[substructure_indices+6]*alignment[:,np.newaxis]
# Overlay shapes
compare_geo, compare_shapes = sdpy.shape.overlay_shapes(
    [combined_geometry,full_geometry],[combined_shapes[substructure_indices+6],full_shapes[truth_indices+6]],color_override=[1,7])
compare_shapes.comment1 = 'Blue: CB Substructured, Green: Truth'
compare_geo.plot_shape(compare_shapes,undeformed_opacity=0,deformed_opacity=0.5)
Mode comparison between the truth and substructured models

Mode shape comparison for the first shape where differences between the truth and substructured shapes are readily visible.

Summary

This example has demonstrated how Craig-Bampton superelement and transformation files from Sierra/SD can be loaded into SDynPy as a system object and then used for substructuring. The results were compared against truth data from an eigensolution of the full system matrices, and very good agreement was obtained.