Airplane Modal Test

This page contains an example demonstrating some usages of SDynPy. In this analysis, we will walk through the entirety of the typical analysis/test/analysis workflow. We will:

1. Load in the finite element model and use it to select instrumentation locations.

2. Simulate a modal test on the test article to collect time data

3. Compute FRFs from the excitation response

4. Fit modes to the frf data

5. Compare modal fits to finite element model

6. Perform finite element expansion using SEREP

7. Create quicklook reports for the test

Imports¶

For this project, we will import the following modules, including the SDynPy module.

Default Plotting Options¶

Since we will be plotting a lot of shapes, we will set up some options to use for all geometry plots. See the documentation for sdpy.Geometry.plot for these options. We will create a dictionary of these options so they can be passed into the various plotting functions using the **kwargs syntax, or into the various plotting functions that accept a plot_kwargs argument.

Load the Finite Element Model¶

For complex modal tests, it can be non-trival to select a good set of sensors. Additionally, for test articles with internal geometry, sensor locations may not be accessible once the unit is built up. Therefore, selecting a good initial sensor set can be crucial to the success of a given test. For this reason, many tests begin by using a finite element model to select a set of degrees of freedom. We will follow this path in this example.

Many analyses performed at Sandia National Laboratories use the Exodus file format to store output data. We can load in an Exodus file into SDynPy using the sdpy.Exodus class. Additionally, we can reduce the size of the finite element model by reducing it to just the outer surfaces using sdpy.Exodus.reduce_to_surfaces, which returns a sdpy.ExodusInMemory object. This is a resonable step to take, as we will not be able to put any sensors on the interior of material volumes.

Convert Finite Element Model to SDynPy Objects¶

At this point, we would like to convert our finite element model into SDynPy objects to make it easier to work with. We will create both a Geometry as well as Shapes.

Creating a Geometry from an Exodus file¶

Our goal is to create a test geometry. We will note that, while finite element models are often defined exclusively in a single, global coordinate system, we often cannot place sensors in the global coordinate system if, for example, there are surfaces oblique to the principal directions of the part. We often place sensors directly on the surfaces, so the sensor coordinate system will generally be oriented such that it is aligned with the local surface normal of the location that the sensor is placed.

We will therefore create local coordinate systems to use for sensor selection. To do this, we use the local=True argument in the sdpy.geometry.from_exodus function. We will specify a preferred direction in the nose-wise direction preferred_local_orientation=[0,0,1]. The secondary preferred direction will be “up” secondary_preferred_local_orientation=[0,1,0].

Exploring the Geometry object¶

Take a moment here to explore the Geometry object. If we simply type geometry into the IPython console after running the previous command, we get out the representation of the Geometry object (truncated here to save space)

We see there are a large number of nodes, coordinate systems, and elements, but no tracelines. We can access these data by accessing the Geometry attributes: geometry.node, geometry.coordinate_system, geometry.element, and geometry.traceline.

Nodes¶

We will start with the NodeArray object revealed by geometry.node. The NodeArray class is a subclass of SdynpyArray, which is itself a subclass of NumPy’s ndarray. All subclasses of SdynpyArray can therefore take advantage of NumPy functions such as intersect1d, unique, or concatenate and also handle indexing and broadcasting identically to the NumPy ndarray.

Subclasses of SdynpyArray store their data internally as a structured array variant of the ndarray. However, as an alternative to accessing the field data using the syntax array['fieldname'], SdynpyArray allows accessing the fields as if they were attributes using the syntax array.fieldname. Many integrated development environments will not recognize these attributes so all SdynpyArray subclasses have a fields attribute that lists the fields stored in the array that can be accessed.

Returning to the geometry.node, we can identify the fields in the object using the command

Here we see the five fields of the NodeArray object. We can obtain even more information about the shape and type of each of these fields using the dtype attribute, which is inherited from NumPy’s ndarray.

Here we see that the geometry.node.id array, which contains the node ID number, is a 8-byte (64-bit) unsigned integer. The geometry.node.disp_cs and geometry.node.def_cs arrays, which contain references to the coordinate system in which the node is defined and in which the node displaces, respectively, are also this data type. The geometry.node.color array, while still an unsigned integer, is only 2 bytes, or 16 bits. Finally, the geometry.node.coordinate, which contains the 3D position of the node as defined in the geometry.node.def_cs, consists 8-byte (64-bit) floating-point data, and also has a shape of (3,), which signifies there are three values of the coordinate for each entry in the geometry.node array. These extra dimensions of the field arrays are appended at the end of dimension of the SdynpyArray subclass. For example, if we compare the shape of the geometry.node array to the geometry.node.coordinate array, we will see that the shapes are identical except for the appending of the length-3 extra dimension on the latter array. Here the shape attribute is also an attribute inherited from NumPy’s ndarray.

We see that the shape of our geometry.node array is 12686, meaning the geometry we are examining has that many nodes. We then see that the shape of our geometry.node.coordinate array is 12686 x 3, showing that there are three coordinate values for each node.

Coordinate Systems¶

Coordinate systems in the Geometry object are stored in a CoordinateSystemArray object that can be accessed by geometry.coordinate_system. We will again explore the fields of the CoordinateSystemArray using the dtype.

We now see some new types of fields. We still have id and color, which are consistent with the NodeArray object we previously explored. We now have another integer field cs_type which stores the type of coordinate system (0 - cartesian, 1 - cylindrical, 2 - spherical) in a 16-bit unsigned integer field. We also have a name field, which stores a name of the coordinate system in a string of less than 40 characters. Finally, there is the coordinate system’s transformation matrix, stored in the matrix field, which is stored in a 4 x 3 array of 64-bit floating point numbers. Again, recall the shape of the fields are appended to the shape of the base object, so comparing the shape of the CoordinateSystemArray to the shape of its matrix field, we will see that the latter has 2 extra dimensions of length 4 and 3.

Elements¶

Elements in the Geometry are stored in an ElementArray object, which can be accessed using the geometry.element attribute. The fields of this object are

Like NodeArray and CoordinateSystemArray objects, the ElementArray object also has id and color fields. Each element also has a type field, which is an 8-bit unsigned integer representing the element type as defined by the universal file format dataset 2412. Finally, the element connectivity field is stored as an object array, where each entry in the element array is a NumPy ndarray with length equal to the number of nodes in the element. This construction is necessary as each element might have a different number of nodes, so a single array of fixed size is not possible.

Tracelines¶

The final visualization tool in the Geometry object is the TracelineArray, which represents a line connecting nodes in the geometry. The fields of the TracelineArray object are

Similarly to the other geometry objects, TracelineArray objects have id and color, and like to the ElementArray object, it has a connectivity array that specifies the node IDs to connect with the line. The description field stores a name or description of each item in the TracelineArray as a string with less than 40 characters.

Visualizing Degrees of Freedom on a Geometry¶

Now that we have an understanding of our Geometry object, we would like to visualize it. In particular, we would like to visualize the local coordinate systems in the geometry to verify that they were created correctly.

To do this, we will create a CoordinateArray object, which contains a set of nodes and directions. Since we would like to generate a coordinate system triad for each node, we will use the helper function from_nodelist which creates degrees of freedom for all directions for each node passed in as an argument. Here we will use the entire list of node ID numbers from our geometry.

Coordinates¶

Here again is a good place to explore what makes up a CoordinateArray object, which we have just assigned to the variable coordinates. We can examine the data type of the CoordinateArray to see that it contains fields for a 64-bit unsigned integer as the node field and an 8-bit signed integer for the direction field.

CoordinateArray objects store the direction as an integer where rotations are encoded:

• $X+ = 1$

• $X- = -1$

• $Y+ = 2$

• $Y- = -2$

• $Z+ = 3$

• $Z- = -3$

• $RX+ = 4$

• $RX- = -4$

• $RY+ = 5$

• $RY- = -5$

• $RZ+ = 6$

• $RZ- = -6$

• (empty) $= 0$

When we want to read CoordinateArray objects, the integer directions are typically transformed into the more readable direction strings shown in the first column of the above table. For example, if we type coordinates into the console, the representation of the object displays the string array version of the coordinates.

From the above, we can see that the coordinate variable we just created contains a degree of freedom for each of the positive X, Y, Z directions at each node.

Many SDynPy objects allow indexing with a CoordinateArray object to automatically handle the bookkeeping aspect of selecting the right data for each coordinate.

Plotting Coordinates¶

At this point, we would like to plot our coordinates on top of our geometry. For this we use the plot_coordinate method of the Geometry object.

Note that due to the density of the mesh, we had to make the arrow_scale much smaller than the default. Here also note that we passed the previously defined plot_options into the plot_kwargs argument. This method should produce a graphic similar to the one shown below.

If we zoom into the coordinate systems on the figure, we see that the Z+ (blue) direction is always oriented perpendicular to the surface, and the X+ (red) direction is always oriented towards the direction specified by the preferred_local_orientation argument to the sdpy.geometry.from_exodus function, which was the global Z+ direction.

Creating Shapes from an Exodus File¶

We just finished extracting geometry with local coordinate systems. We will now extract shapes from the sdpy.ExodusInMemory object. To do this, we will use the sdpy.shape.from_exodus function. Note that shapes gathered from this function will be in the global coordinate system (as defined in the model), so we will assign the variable a name that denotes that these are global shapes.

Shapes¶

At this point, it is useful to explore briefly the ShapeArray object in the IPython console. The data type of the object is:

The data type of ShapeArray objects can change depending on what type of shape and how many degrees of freedom are in the shape. frequency and damping fields are stored as 64-bit floating point numbers with one value per entry in the ShapeArray. modal_mass is also stored in the present ShapeArray, but if the shape is complex, then the modal mass might also be complex. The shape_matrix field holds the underlying shape data. It has one entry for every degree of freedom in the shape, and is represented by a floating point number for normal modes or a complex number for complex modes. Similarly, the coordinate field identifies which degree of freedom belongs to which entry in the shape_matrix field. The coordinate field stores data as CoordinateArray objects, and thus has the same data type as CoordinateArray. Finally, there are five fields available for comments, which store string data up to 80 characters which can be used to store any data the user feels is relevant to the analysis.

One thing to note is that the shape_matrix field, due to the dimension of the field being appended at the end of the array, will be transposed from the typical representation of a mode shape matrix (degrees of freedom as rows and mode indices as columns). The shape_matrix field will instead have the shape of the ShapeArray object itself as its first dimensions, and then the size of the coordinate field as its last dimension.

Correcting Negative Frequencies¶

One thing to note is that sometimes finite element codes that solve for zero frequency modes can have slight errors that result in small negative natural frequencies. These can, however, cause issues due to stiffness instabilities if, for example, the system is integrated over time. It is easy to set all of these negative natural frequencies to positive natural frequencies.

Note carefully here the order of operations in the above statement. We have a ShapeArray object containing our shape data. We then access the frequency field of that ShapeArray object and index that array to select only those frequencies that are less than zero and set them to zero. Note the subtle difference to the following line of code.

In the above line of code, we first index our ShapeArray object to access only the shapes with frequency less than zero, then access the frequency field of just those shapes and assign to zero. This seems like it should be the identical operation, correct? Indeed, both operations return the same negative frequency values.

However, experienced users of NumPy should readily identify the subtle difference between the above two operations. In the second, the shapes_global ShapeArray is indexed by a logical array (a NumPy ndarray consisting of True and False, True where the frequency field is less than 0). This will trigger NumPy’s advanced indexing schemes, which always return a copy of rather than a view into the underlying data. Therefore, in the second case, what we are looking at is the frequency field of a copy of our shapes_global variable rather than the raw data. Similarly, assigning directly to this copy will not change the data of the original shapes_global variable.

Transforming Shape Coordinate Systems¶

As mentioned above, the shapes that come from the sdpy.shape.from_exodus function are defined in the global cordinate system. We can see that if we plot the shapes on the geometry we loaded previously using the plot_shape method of the Geometry object, the shapes do not look right, due to plotting shapes defined using a global coordinate system onto a geometry defined using local coordinate systems. For example, the shape below should be a rigid body mode of the system, but instead appears to be a torsional mode of the system.

What we need to do is transform the global shapes into the local shapes defined on the geometry. We can do this using the transform_coordinate_system method of the ShapeArray object. This method takes two arguments; the first is a Geometry object that defines the coordinate systems that we are transforming “from” and the second is a Geometry object that defines the coordinate systems that we are transforming “to”. We already have the “to” geometry in the form of the local geometry we loaded from the Exodus file. We therefore only need to create a global geometry corresponding to the global coordinate system that our shapes are currently defined in to transform “from”. This is easily done using the sdpy.geometry.from_exodus function while not passing in the arguments that specify the local coordinate systems should be created. Plotting the same coordinates CoordinateArray on this new geometry will verify that all node coordinate systems are aligned with the global coordinate system. We can also plot shapes_global with geometry_global to verify that the shapes look right.

At this point, we have the requirements to transform coordinate systems. We can now plot the local shapes with the local geometry to verify that the shapes are indeed correct. These shapes should look identical to the above figure even though the underlying data is represented in a different coordinate system.

Optimizing Instrumentation for Test¶

The reason we went through all the trouble to create shapes in local coordinate systems is that these are the coordinate systems in which test instrumentation is generally placed. We will now want to downselect a set of sensors to use in the actual test from a candidate set of instruments.

For this example, we will use the effective independence algorithm implemented in sdpy.dof.by_effective_independence. This function starts with a candidate set of degrees of freedom and iteratively throws away the degrees of freedom with the smallest contribution to the effective independence of the system. Because there are currently a lot of nodes in our geometry (12,686), it will take a long time to consider each one of these locations as a potential sensor in a test. Instead, we will reduce the candidate set by overlaying a grid of points over the test geometry with a specified spacing using the node_by_global_grid method.

We have reduced our initial candidate set of 12,686 nodes down to 1,217 nodes.

We can visualize this initial grid on our model by creating a CoordinateArray object from these nodes using the sdpy.coordinate_array function and broadcasting the node identification numbers with the three directions. We can then plot these degrees of freedom on our geometry.

Now that we have selected our candidate degrees of freedom, we need to select our target shapes. Here we will consider the bandwidth of the test plus a bit more to capture some effects of out-of-band modes. We are interested in data up to 200 Hz for this test, so we will consider shapes up to 300 Hz for this instrumentation selection. We can use logical indexing to perform this frequency reduction, keeping only shapes where the frequency field is less than the desired bandwidth.

At this point, we will need to create a ShapeArray object containing only the candidate nodes, which will be the starting point of our optimization. We can do this easily using ShapeArray.reduce and passing the array of node identification numbers to keep.

We can now run the sensor selection algorithm using the effective independence objective function. While the sensor selection schemes are contained in sdpy.dof, these capabilities can be accessed more easily through the ShapeArray.optimize_degrees_of_freedom method. We give it our sensors_to_keep and shape_matrix, as well as give the optional flag to return_efi so we can visualize the effective independence of the shape matrix as degrees of freedom are removed. We can plot this function, as well as the kept degrees of freedom by indexing the oridinal set of degrees of freedom with the keep_indices returned from the by_effective_independence function. Intuitively, the sensors seem well-placed, with many positioned at the boundaries of the wing and tail, as well as a few on the nose and fuselage.

Creating a Test Geometry¶

Because the test sensors are a small subset of the nodes on the original finite element model, it can be useful to add visualization features to the test geometry to help visualize responses from the test. We have already explored elements, which our finite element model has in abundance, so for the test geometry we will use tracelines to aid in visualization. We will use the helper function add_traceline.

The first step in this process is to simply reduce the original geometry down to our kept sensor locations. using the reduce method. This method will discard any coordinate systems, tracelines, or elements that do not entirely consist of the nodes that are kept, which in this case means that the geometry will consist of only nodes and their coordinate systems.

At this point, we would like to connect the nodes with tracelines. We can plot the kept sensor degrees of freedom on the reduced test_geometry, and specify label_dofs=True to create labels at each of the coordinates. This will allow us to easily see which nodes should be connected. Alternatively, we can simply plot the geometry and label the nodes.

With the points labeled, it is easy to see which nodes should be connected using tracelines. When we create tracelines, we will also set the colors, using color=1 (blue) for the fuselage, color=7 (green) for the wings, and color=13 (pink) for the tail. After the tracelines are added, it is much easier to visualize the structure of the airplane.

With this new test geometry, we should be able to visualize the reduced test shapes. We should also plot the Modal Assurance Criterion (MAC) matrix with the ShapeArray.plot_mac method to ensure we can distinguish the shapes of interest with our reduced sensor set.

Ooops! An Instrumentation Error!¶

Installing instrumentation is a tedious process that must be well-documented for success of a given test. For large channel-count tests, it is almost inevitable that there will be some kind of sensor or geometry error, so it is useful to develop workflows to verify that sensors are oriented correctly. Here we will purposefully introduce an instrumentation error by setting the displacement coordinate system of one of our sensors to the global coordinate system rather than its correct local coordinate system. We will then investigate a workflow in SDynPy for identifying correcting this error.

Adding Damping¶

Before we start simulating tests, we will add damping to our test shapes. Our finite element model had no damping prescribed, but all real test articles will have damping. To keep it simple, we will add a blanket 2% damping to all modes in the model.

Running a Virtual Experiment: Rigid Body Checkouts¶

One approach to validating the test geometry and channel table is to perform rigid body checks, where the structure is excited at a frequency well below the elastic modes of the system to elicit a rigid body response. These response shapes can be visualized on the geometry and should appear rigid. If, for example, a sensor is moving out of phase with the rest of the sensors, it is fairly certain there is an instrumentation or bookkeeping error. When looking at acceleration response, these rigid shapes should also have very small imaginary parts compared to their real parts. Finally, a projector can be used to project the measured rigid body shapes through analytic rigid body shapes created from the test geometry. This operation will remove any “non-rigid” components of the measured shapes. By subtracting this “rigidized” version from the measured shapes, a residual can be created where large values in the residual point to non-rigid motions.

Creating a System Object to Integrate¶

We will now simulate one of these rigid body tests by integrating equations of motion. SDynPy uses System objects to represent dynamical systems.
System objects contain mass, stiffness, and damping matrices. They also can contain a transformation matrix to return to physical coordinates if the internal state degrees of freedom are represented in a reduced space (for example, a mode shape matrix is the transformation between modal mass, stiffness, and damping matrices and the corresponding physical coordinates). System objects also contain a CoordinateArray to aid in bookkeeping.

We can easily create a modal System from our existing shapes by using the ShapeArray.system method. This will return a state matrix in modal space. We can use the System.spy method to visualize the structure of the created system.

We can see from the above matrix plot that the internal state is diagonal, consisting of 43 modal degrees of freedom, whereas the transformation matrices (which are the mode shape matrices in this case) are full, consisting of 90 physical degrees of freedom.

Exploring the System Object¶

Here is a good place to take time to investigate the System object we have just created. We can look at the mass, stiffness, and damping attributes to see those matrices, and verify that they are indeed diagonal for the modal system we have created.

We can also look at the transformation and coordinate to see that the bookkeeping from internal state (modal) degree of freedom to physical degree of freedom is maintained.

Setting up the Integration and Forcing Function¶

We aim integrate our test_system using the System.time_integrate method, but first we need to set up its required sampling parameters. We will specify our test bandwidth, as well as an integration oversample factor to ensure that the integration is performed accurately.
System objects represent Linear, Time-Invariant systems, so the scipy.signal.lsim function is used for integration within the System.time_integrate function, and this typically only requires a factor of 10x the desired sampling rate to achieve reasonable results. We will also define the desired frequency spacing for the test. The other sampling parameters can be derived from these definitions. We will measure 10 frames so averaging can be performed when computing frequency response functions.

We will now create the forcing function for our rigid body test, which will be a sine wave. The first elastic mode of the test article is near 6 Hz, so we choose a frequency that is much lower than that value. We can then use the sdpy.generator.sine function to produce a sine wave with the proper sampling parameters. However, a more straightforward approach is to use the sdpy.TimeHistoryArray.sine_signal method, which automatically constructs a TimeHistoryArray object. The functions in sdpy.generator only produce NumPy arrays which must then be turned into TimeHistoryArray objects.

As we construct these signals, we must also select the degrees of freedom at which the forces will be applied. We will select degrees of freedom that largely point through the airplane’s center of mass to ensure largely translational response.

The sine_signal method accepts arguments such as amplitude, phase, and various comments; however, we will leave these as the default values for now.

Running the Integration to Generate Synthetic Test Data¶

We will now loop through each of our excitation degrees of freedom and integrate the System’s response to the force located at that position using System.time_integrate. This will give us responses as TimeHistoryArray objects. We will then truncate off the first second of the response in order to remove the startup transients from the data using TimeHistoryArray.extract_elements_by_abscissa. We will then pass the truncated data into the sdpy.TransferFunctionArray.from_time_data function to produce frequency response functions as a TransferFunctionArray object. Once the frequency response functions are computed, we can extract the deflection shapes to get the deflection shape using to_shape_array, which will produce a ShapeArray object. We append each shape to a running list, and then concatenate all the shapes into one ShapeArray object, which we can plot on our erroneous geometry.

When performing time integration, we must specify which degrees of freedom we would like responses to be generated at. Additionally, we need to specify what type of response (displacement, velocity, or acceleration), we would like. We do this using a dictionary syntax where the keys are the displacement derivative (e.g. 2 for acceleration, which is the second derivative of displacement). We would like all acceleration responses, which means we should populate the response dictionary key of 2 and the value as the degrees of freedom we would like to keep.

When we plot the measured rigid shapes against the geometry where we intentionally introduce errors, we see that one of the sensors on the wing seems to be moving opposite the rest of the wing. We can see that by visualizing rigid body response, we are able to identify incorrect test geometry. We will investigate this more thoroughly in a moment.

Exploring Data Objects¶

Before we move on to analyzing the rigid body datasets, we will first take some time to explore some of the data objects that we have created. Namely, we generated some TimeHistoryArray objects and a TransferFunctionArray object for each excitation degree of freedom.

In SDynPy, all data objects inherit from the NDDataArray class, which represents a general, multi-dimensional data container. This class, in turn, inherits from the SdynpyArray, so all the broadcasting and attribute access features are available as well.

We will start by examining responses from the previous block of code. We see that it is a TimeHistoryArray object. We can look at its dtype to see what data it contains.

Being a TimeHistoryArray object, responses consists of real data. It has floating point abscissa and ordinate fields to contain the independent time and dependent response variables, respectively. Note that the abscissa and ordinate fields have the size of the length of the data record; in this case, it the time data consists of 316,000 samples. It also has five 80-character string fields where comments can be stored. Finally, it has a coordinate field that stores the degree of freedom information for each data record. Note the rather peculiar length 1 on the coordinate field. This is to signify that a time history signal can be considered a 1D data array (one degree of freedom per signal), which we will see is contrary to, for example, TransferFunctionArray objects which are 2D, having both a response (output) coordinate and a reference (input) coordinate for each signal. Again, remember that these extra field dimensions are appended to the TimeHistoryArray dimensions.

Consider the previous object in contrast to the frf variable, which is a TransferFunctionArray object.

TransferFunctionArray has all the same fields as TimeHistoryArray, except they are different shapes and types. Because frequency response data is complex, the ordinate field is now a 16-byte complex number rather than an 8-byte floating point number. Similarly, because there are now reference and response coordinates associated with each data record, the length of the coordinate field is now 2.

Identifying Bad Geometry with Rigid Body Checkouts in SDynPy¶

Now that we’ve understood the data objects a bit better, we can return to the task at hand, which is to sort out our geometry and channel table, which has an intentionally incorrect sensor.

SDynPy has some built-in tools for doing rigid body checkouts, plotting the complex plane and the residual discussed previously. These are contained in the sdpy.shape.rigid_body_check function. Simply give the function the geometry and shapes, and it will attempt to figure out which nodes are suspicious and warrant further investigation.

The first three plots generated by sdpy.shape.rigid_body_check are the complex plane. These plots allow us to visualize the phase of the responses compared to the reference signal. For a force-to-acceleration FRF, we should see the deflection shape be primarily real when far from elastic mode shapes, meaning the imaginary part should be small compared to the real part. Indeed, we see that the imaginary part is several orders of magnitude smaller than the real part.

The last plot generated by sdpy.shape.rigid_body_check is the rigid body residual plot. This is a measure of how “rigid” the deflection shapes look. SDynPy does a best fit of the deflection shapes from the test to a set of rigid body shapes created analytically from the geometry. Sensors that don’t match their analytically created counterparts will show a high residual. This analysis has clearly identified the degrees of freedom at node 6214 as having high residuals, and this is exactly the sensor we applied the error to.

While the example problem shown here only has external sensors, many times internal sensors may be suspicious, and without visual access to them, it can be difficult to investigate the potential incorrect orientations of the sensors. For this reason, SDynPy includes an approach to identify the “best” orientation of the sensor using the sdpy.shape.rigid_body_fix_node_orientation function. If we give this function a geometry, a set of rigid shapes, and a list of suspicious nodes, it will attempt to find the correct orientation of of the sensors at the suspicious nodes.

We see that when we plotted the coordinate systems, SDynPy was able to take the initial erroneous coordinate system (left) and correct it (center) so that it matched the original coordinate system (right) that we had before we introduced the instrumentation error.

We see now that when we plot the shapes on the corrected geometry, it indeed looks like rigid motion.

Running a Virtual Experiment: Modal Testing¶

Now that we have validated test geometry, it is time to perform a modal test. Again, we will simulate this virtually for this test.

The first thing we will do is select drive points to get all of the modes of the system. We will simply select positions here intuitively. We will select wing tips, tail tip, and nose points.

Now we would like to set up a force for our modal test. Here we will use a random excitation, which will enable us to excite the structure with all drive points simultaneously, which results in a Multiple-Input, Multiple-Output modal test.

We can set up a random force using the sdpy.generator.random function or the sdpy.TimeHistoryArray.random_signal method. Again, we will choose the latter because it automatically provides a sdpy.TimeHistoryArray object instead of a np.ndarray. We will provide it the time step size, the total number of samples, the degrees of freedom, and a high-frequency cutoff to ensure that we don’t have aliasing when we downsample from the integration oversampling factor back to our test bandwidth. We will accept default parameters for rms level and other parameters.

We should plot the signal and its frequency content to verify we have set it up correctly.

We can easily compute the FFT of the signal using sdpy.TimeHistoryArray.fft. We should see that we only have frequency content up to the test bandwidth of 200 Hz.

We will again do integration of our test system using the System.time_integrate function, and we will then downsample the resulting data to our test bandwidth. We can again pass the reference and response data into the sdpy.TransferFunctionArray.from_time_data function to compute FRFs. Note that because we now have random data, we will use a window function and apply an overlap.

Since we now have a large number of FRFs, it can be difficult to visualize them all simultaneously, so we will use SDynpy’s GUIPlot to allow us to quickly look through all the functions to ensure they look right.

Fitting Modes using PolyPy¶

Now that we have frequency response functions created, we can fit modes to them. SDynPy has two mode fitters implemented, PolyPy and SMAC. Both curve fitters can be used via graphical user interface or via Python commands if it is desirable to automate the curve fitting. This example will use the sdpy.PolyPy_GUI approach.

Running PolyPy¶

We open the PolyPy GUI by initializing the sdpy.PolyPy_GUI class with our frequency response function dataset frf_sampled. We can also set the geometry used to plot shapes here. Alternatively, we could load the geometry from disk through the GUI.

The initial screen shows mode indicator functions, as well as options for computing the initial stabilization diagram. We can see from the shown Complex Mode Indicator function that there are a few instances of closely-spaced modes. We can drag the frequency region on the figure to select the frequency range of interest, set the polynomial orders, and press the button to compute the stabilization curve.

Once the stabilization plot is computed, stable poles can be selected by clicking on them in the stabilization plot. Once all poles are selected, shapes can be computed.

If the Auto-Resynthesize Checkbox is checked, then PolyPy will compute shapes and resynthesize after each pole is selected. For large datasets, this can be computationally intensive, so the checkbox can be unchecked. If the box is not checked, the Resynthesize button can be clicked to manually trigger a resynthesis.

To visualize the resynthesis quality, we can click on one of the buttons in the Resynthesis portion of the Select Poles tab. For example, clicking the CMIF button will compute CMIFs for the experimental and resynthesized FRFs. In this case, since the data is synthetic, we should expect very good fits.

To visualize fit shapes, the Plot Shapes button can be clicked. If we hadn’t loaded a geometry already we would have had to press the Load Geometry button to load a geometry onto which the shapes could be plotted.

When we are satisfied with our fit shapes, we can save them to disk by clicking the Export Shapes... button and choosing a file name. Here we will chose shapes_polypy.npy as the file name.

Comparing Test and Finite Element Modes¶

Once modes are fit in the GUI and saved to disk, we can load them back into our analysis. We initially compute a Modal Assurance Criterion matrix between the fit shapes and our finite element shapes to compare the results. We can also print a shape comparison table obtained by the sdpy.shape.shape_comparison_table function.

We can see that we achieve very good comparisons to the analytical shapes. Obviously we have not fit the rigid body shapes in the test, so there is correlation between the first six FEM modes and any of the test shapes. Similarly the bandwidth of our test was truncated, so we do not see correlation of test shapes to the higher frequency FEM modes. One thing to note is the increase of damping at low frequency in our fit modes from the simulated test. This is due to the application of the Hann window when computing FRFs; the window function artificially pushes the signals to zero, which artificially increases the damping. This affects lower frequency modes more significantly because they tend to take longer to damp out for an equivalent damping ratio. Note that all window functions disrupt data, they just disrupt data less than the leakage that would otherwise be obtained when processing the random signals!

Another way to compare shapes is to overlay them. This is easily done within SDynPy using the sdpy.shape.overlay_shapes function. The output of this function is a “combined” geometry and “combined” shape. The nodes, coordinate systems, elements, and tracelines are all combined into one geometry, and the id values of each of these elements is offset to avoid conflicts. The shape degrees of freedom are also concatenated and offset similarly to produce the appropriate shape visualization.

Note that our test shapes might be 180 degrees out of phase with the finite element model shapes, so we will first compute the dot product between the two sets of shapes and flip the sign on any shape where the dot product is negative. We can do this easily using the shape_alignment method.

SEREP Expansion¶

Often times, the “stick” test geometry can be insufficient for communicating results to test stakeholders, particularly when there are only uniaxial sensors on the test. Among other things, the System Equivalent Reduction Expansion Process (SEREP) can be useful for expanding data at test sensors out to the full finite element space for better visualization. SDynPy makes it easy to perform SEREP using the ShapeArray.expand method. Particularly in this case where we have kept the test geometry node IDs equivalent to the original finite element node IDs, the expansion bookkeeping is handled automatically.

We first need to create the set of shapes to use to expand the test shapes. These will be the finite element shapes in the test bandwidth. We then call the ShapeArray.expand method, giving it the test geometry, finite element geometry, and expansion basis shapes. Note here that the global shapes and global geometry are used in the expansion. All coordinate transformations between the local test geometry and global finite element geometry are handled automatically by SDynPy.

Summary¶

This document performed a walkthrough of a typical experimental modal analysis. A FEM was loaded in to perform test design studies. Optimized instrumentation locations were selected from this model.

We then created a test geometry to enable us to plot the eventual mode shapes from the test. We noted that instrumentation is tedious and error prone, so we saw how SDynPy gives us the ability to identify incorrectly documented sensors in the channel table or test geometry and a way to fix improperly oriented sensors in the geometry if we don’t have physical access to them to verify their orientation. For these studies we used a “rigid body” test where we excited the structure with low frequency sine waves and measured the rigid response.

After we validated the geometry, we performed a modal test. We did a MIMO Random modal test with four shaker locations. We computed FRFs from the time data and then fit modes. We noted that the PolyPy mode fitter was able to handle closely spaced modes and was able to fit all the modes in the bandwidth.

We then compared the modes from the synthetic test data back to the finite element model, noting that we got exceptionally good fits to the data. We also demonstrated finite element expansion using SEREP.