digital-brain/hatch00_vibrat_matlab_ansys.md at dd52c291671497ab975963af23b0e733df7b1daf

Thomas Dehaeze 22fb3361a5 Update Content - 2022-03-15

2022-03-15 16:40:48 +01:00

74 KiB

Raw Blame History

+++ title = "Vibration Simulation using Matlab and ANSYS" author = ["Dehaeze Thomas"] description = "Nice techniques to analyze resonant systems with Ansys and Matlab." keywords = ["Modal Analysis", "FEM"] draft = false +++

Tags: [Finite Element Model]({{< relref "finite_element_model.md" >}})
Reference: (Hatch 2000)
Author(s): Hatch, M. R.
Year: 2000

Matlab Code form the book is available here.

Introduction

The main goal of this book is to show how to take results of large dynamic finite element models and build small Matlab state space dynamic mechanical models for use in control system models.

The diagram in Figure 1 shows the methodology for analyzing a lightly damped structure using normal modes.

The steps are:

deriving the undamped equations of motion in physical coordinates
solving the eigenvalue problem, yielding eigenvalues (natural frequencies) and eigenvectors (mode shapes). These first 2 steps are usually performed using a Finite Element software
transform the model from physical coordinate system to modal or principal coordinate system using the eigenvector matrix. In the modal coordinate system, the original undamped coupled equations of motion are transform to the same number of undamped uncoupled equations representing the motion of a particular mode of vibration of the system.
proportional damping is applied
the uncoupled equations are easily solved (as each equation corresponds to a sdof system). The desired responses are then back-transformed into the physical coordinate system using the eigenvector matrix.

{{< figure src="/ox-hugo/hatch00_modal_analysis_flowchart.png" caption="Figure 1: Modal analysis method flowchart" >}}

Model Size Reduction

Because finite element models usually have a very large number of states, an important step is the reduction of the number of states while still providing correct responses for the forcing function input and desired output points.

Figure 2 shows such process, the steps are:

start with the finite element model
compute the eigenvalues and eigenvectors (as many as dof in the model)
the goal is then the reduce the size of the model while still maintaining the desired input/output relationships. This is done in two steps:
1. reduce the number of dof of the model from the original set to a new set which includes only those dof where forces are applied and where responses are desired
2. reduce the number of modes of vibration used for the solution. For SISO system, the modes are ranked based on the relative importance to the overall response. For MIMO system, controllability and observability gramians of the modes are estimated.

{{< figure src="/ox-hugo/hatch00_model_reduction_chart.png" caption="Figure 2: Model size reduction flowchart" >}}

Notations

Tables 3, 2 and 3 summarize the notations of this document.

Table 1: Notation for the modes and nodes

	Notation
Number of Nodes	\(n\)
Number of Considered node inputs	\(n_i\)
Number of Considered node outputs	\(n_o\)
Reduced number of nodes	\(n_r\)
Number of Modes	\(m\)
Reduced number of modes	\(m_r\)

Table 2: Notation for the dofs, eigenvectors and eigenvalues

	Notation	Size
Physical dofs	\(\bm{z}\)	\(n \times 1\)
Principal dofs	\(\bm{z}_p\)	\(m \times 1\)
Modal Matrix	\(\bm{z}_m\)	\(n \times m\)
Normalized Modal Matrix	\(\bm{z}_n\)	\(n \times m\)
i'th eigenvector	\(\bm{z}_{mi}\)	\(n \times 1\)
i'th eigenvalue	\(\omega_i\)	\(1 \times 1\)

Table 3: Notation for the mass and stiffness matrices

	Notation	Size
Physical Mass, Stiffness, Damping Matrices	\(\bm{m}\), \(\bm{c}\), \(\bm{k}\)	\(n \times n\)
Normalized Mass and Stiffness Matrices	\(\bm{m}_n\), \(\bm{k}_n\)	\(m \times m\)
Principal Mass, Stiffness, Damping Matrices	\(\bm{m}_p\), \(\bm{c}_p\), \(\bm{k}_p\)	\(m \times m\)
Physical Force Vector	\(\bm{F}\)	\(n \times 1\)
Principal Force Vector	\(\bm{F}_p\)	\(m \times 1\)

Zeros in SISO Mechanical Systems

The origin and influence of poles are clear: they represent the resonant frequencies of the system, and for each resonance frequency, a mode shape can be defined to describe the motion at that frequency.

We here which to give an intuitive understanding for when to expect zeros in SISO mechanical systems and how to predict the frequencies at which they will occur.

Figure 3 shows a series arrangement of masses and springs, with a total of \(n\) masses and \(n+1\) springs. The degrees of freedom are numbered from left to right, \(z_1\) through \(z_n\).

{{< figure src="/ox-hugo/hatch00_n_dof_zeros.png" caption="Figure 3: n dof system showing various SISO input/output configurations" >}}

(Miu 1993) shows that the zeros of any particular transfer function are the poles of the constrained system to the left and/or right of the system defined by constraining the one or two dof's defining the transfer function.

The resonances of the "overhanging appendages" of the constrained system create the zeros.

State Space Analysis

Lightly damped structures are typically analyzed with the "normal mode" method described in this section.

The modal method allows one to replace the n-coupled differential equations with n-uncoupled equations, where each uncoupled equation represents the motion of the system for that mode of vibration.

The overall response of the system is then reconstructed as a superposition of the responses of the different modes of the system.

Heavily damped structures or structures which explicit damping elements, such as dashpots, result in complex modes and require state space solution techniques using the original coupled equations of motion. Thus, the present methods only works for lightly damped structures.

Summarizing the modal analysis method of analyzing linear mechanical systems and the benefits derived:

Solve the undamped eigenvalue problem, which identifies the resonant frequencies and mode shapes (eigenvalues and eigenvectors), useful in themselves for understanding basic motions of the system
Use the eigenvectors to uncoupled or diagonalize the original set of coupled equations, allowing the solution of n-uncoupled sdof problems instead of solving a set of n-coupled equations
Calculate the contribution of each mode to the overall response. This also allows one to reduce the size of the problem by eliminating modes that cannot be excited and/or modes that have not output at the desired dof. Also, high frequency modes that have little contribution to the system at lower frequencies can be eliminated or approximately accounted for, further reducing the size of the system to be analyzed
Write the system matrix \(\bm{A}\), by inspection. Assemble the input and output matrices, \(\bm{B}\) and \(\bm{C}\), using appropriate eigenvector terms. Frequency domain and forced transient response problems can be solved at this point. If complete eigenvectors are available, initial condition transient problems can also be solved. For lightly damped systems, proportional damping can be added, while still allowing the equations to be uncoupled.

Eigenvalue Problem

Equation of Motion

Let's consider the model shown in Figure 4 with \(k_1 = k_2 = k\), \(m_1 = m_2 = m_3 = m\) and \(c_1 = c_2 = 0\).

{{< figure src="/ox-hugo/hatch00_undamped_tdof_model.png" caption="Figure 4: Undamped tdof model" >}}

The equations of motions are:

\begin{equation} \begin{bmatrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & m \end{bmatrix} \begin{bmatrix} \ddot{z}_1 \\ \ddot{z}_2 \\ \ddot{z}_3 \end{bmatrix} + \begin{bmatrix} k & -k & 0 \\ -k & 2k & -k \\ 0 & -k & k \end{bmatrix} \begin{bmatrix} z_1 \\ z_2 \\ z_3 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} \label{eq:tdof_eom} \end{equation}

Principal Mode Definition

Since the system is conservative (it has no damping), normal modes of vibration will exist.

Having normal modes means that at certain frequencies all points in the system will vibrate at the same frequency and in phase, i.e., all points in the system will reach their minimum and maximum displacements at the same point in time.

Having normal modes can be expressed as:

\begin{equation} \bm{z}_i = \bm{z}_{mi} \sin(\omega_i t + \phi_i) = \bm{z}_{mi} \text{Im}(e^{j\omega_i t + \phi_i}) \label{eq:principal_mode} \end{equation}

where:

\(\bm{z}_i\): vector of displacement for all dof's at the i'th frequency
\(\bm{z}_{mi}\): i'th eigenvector
\(\omega_i\): i'th eigenvalue
\(\phi_i\): an arbitrary initial phase angle

Eigenvalues / Characteristic Equation

Re-injecting normal modes eq:principal_mode into the equation of motion eq:tdof_eom gives the eigenvalue problem:

\begin{equation} (\bm{k} - \omega_i^2 \bm{m}) \bm{z}_{mi} = 0 \end{equation}

Solving this set of equation gives the eigenvalues:

\begin{equation} \omega_1 = 0, \quad \omega_2 = \pm\sqrt{\frac{3k}{m}}, \quad \omega_3 = \pm\sqrt{\frac{k}{m}} \end{equation}

Eigenvectors

To obtain the eigenvectors of the systems (corresponding to the eigenvalues), any one of the dof, say \(z_1\), is selected as a reference. This means that the eigenvectors are only known as ratios of displacements, not as absolute magnitudes. Then, all but one of the equations of motion is written with that value on the right-hand side:

\begin{equation} (\bm{k} - \omega_i^2 \bm{m}) \bm{z}_{mi} = 0 \end{equation}

One then find:

\begin{equation} \bm{z}_1 = \begin{bmatrix} 1 \\ 1 \\ 1 \end{bmatrix}, \quad \bm{z}_2 = \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix}, \quad \bm{z}_3 = \begin{bmatrix} 1 \\ -2 \\ 1 \end{bmatrix} \end{equation}

Virtual interpretation of the eigenvectors are shown in Figures 5, 6 and 7.

{{< figure src="/ox-hugo/hatch00_tdof_mode_1.png" caption="Figure 5: Rigid-Body Mode, 0rad/s" >}}

{{< figure src="/ox-hugo/hatch00_tdof_mode_2.png" caption="Figure 6: Second Model, Middle Mass Stationary, 1rad/s" >}}

{{< figure src="/ox-hugo/hatch00_tdof_mode_3.png" caption="Figure 7: Third Mode, 1.7rad/s" >}}

The modal matrix is an \(n \times m\) matrix with columns corresponding to the \(m\) system eigenvectors as shown in Eq. eq:modal_matrix

\begin{equation} \bm{z}_m = \begin{bmatrix} \bm{z}_1 & \bm{z}_2 & \bm{z}_3 \end{bmatrix} = \begin{bmatrix} z_{m11} & z_{m12} & z_{m13} \\ z_{m21} & z_{m22} & z_{m23} \\ z_{m31} & z_{m32} & z_{m33} \end{bmatrix} \label{eq:modal_matrix} \end{equation}

Uncoupling the Equations of Motion

At this point, the system is well defined in terms of natural frequencies and modes of vibration. If any further information such as transient or frequency response is desired, solving for it would be laborious because the system equations are still coupled.

It is thus useful to transform the n-coupled second order differential equations to n-uncoupled second order differential equations by transforming from the physical coordinate system to a principal coordinate system.

In linear algebra terms, the transformation from physical to principal coordinates is known as a change of basis.

There are many options for change of basis, but we will show that when eigenvectors are used for the transformation, the principal coordinate system has a physical meaning: each of the uncoupled sdof systems represents the motion of a specific mode of vibration.

The n-uncoupled equations in the principal coordinate system can then be solved for the responses in the principal coordinate system using the well known solutions for the single dof systems. The n-responses in the principal coordinate system can then be transformed back to the physical coordinate system to provide the actual response in physical coordinate.

This procedure is schematically shown in Figure 8.

{{< figure src="/ox-hugo/hatch00_schematic_modal_solution.png" caption="Figure 8: Roadmap for Modal Solution" >}}

The condition to guarantee diagonalization is the existence of n-linearly independent eigenvectors, which is always the case if either:

the mass and stiffness matrices are both symmetric
there are m-different eigenvalues (no repeating eigenvalues)

We start with the Homogeneous equation of motion:

\begin{equation} \bm{m} \bm{\ddot{z}} + \bm{k} \bm{z} = 0 \end{equation}

Consider normal modes: \(\bm{z}_i = \bm{z}_{mi} \sin(\omega_i t + \phi_i)\):

\begin{equation} \bm{k} \bm{z}_{mi} = \omega_i^2 \bm{m} \bm{z}_{mi} \end{equation}

Since \(\bm{m}\) and \(\bm{k}\) are symmetrical (\(\bm{m}^T = \bm{m}\), \(\bm{k}^T = \bm{k}\)), we find:

\begin{equation} (\omega_i^2 - \omega_j^2) \bm{z}_{mj}^T \bm{m} \bm{z}_{mi} = 0 \end{equation}

When \(i \neq j\), the term \((\omega_i^2 - \omega_j^2)\) cannot be equal to zero, meaning that:

\begin{equation} \bm{z}_{mj}^T \bm{m} \bm{z}_{mi} = m_{ij} = 0 \end{equation}

where \(m_{ij}\) is an off-diagonal term in the mass matrix of the principal coordinate system. The two eigenvectors \(\bm{z}_{mi}\) and \(\bm{z}_{mj}\) are said to be orthogonal with respect to \(\bm{m}\).

For \(i = j\), \((\omega_i^2 - \omega_j^2) = 0\), and thus the product \(\bm{z}_{mj}^T \bm{m} \bm{z}_{mi}\) can be set equal to any arbitrary constant \(m_{ii}\):

\begin{equation} \bm{z}_{mi}^T \bm{m} \bm{z}_{mi} = m_{ii} \end{equation}

This is where various normalization techniques for eigenvectors come into play. The stiffness matrix \(\bm{k}\), is normalize the same manner.

Normalizing Eigenvectors

Because eigenvectors are only known as ratios of displacements, not as absolute magnitudes, we can choose how to normalize them.

Normalizing with Respect to Unity

One method is to normalize with respect to unity, making the largest element in each eigenvector equal to unity by dividing each column by its largest value.

\begin{equation} \bm{z}_m = \begin{bmatrix} 1 & 1 & 1 \\ 1 & 0 & -2 \\ 1 & -1 & 1 \end{bmatrix} \Longrightarrow \bm{z}_n \begin{bmatrix} 1 & 1 & -0.5 \\ 1 & 0 & 1 \\ 1 & -1 & -0.5 \end{bmatrix} \end{equation}

where the \(n\) in \(\bm{z}_n\) refers to a "normalized" modal matrix.

Transforming the mass and stiffness matrices give:

\begin{equation} \bm{m}_n = \bm{z}_n^T \bm{m} \bm{z}_n = \begin{bmatrix} 3m & 0 & 0 \\ 0 & 2m & 0 \\ 0 & 0 & 1.5m \end{bmatrix}; \quad \bm{k}_n = \bm{z}_n^T \bm{k} \bm{z}_n = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 2k & 0 \\ 0 & 0 & 4.5k \end{bmatrix} \end{equation}

Normalizing with Respect to Mass

Another method is to normalize with respect to mass using:

\begin{equation} \bm{z}_{ni}^T \bm{m} \bm{z}_{ni} = 1 \end{equation}

making each diagonal mass term equal to 1. This is the method used by default in ANSYS.

Each normalized eigenvector is defined as follows:

\begin{equation} \bm{z}_{ni} = \frac{\bm{z}_{mi}}{\sqrt{\bm{z}_{mi}^T \bm{m} \bm{z}_{mi}}} = \frac{\bm{z}_{mi}}{q_i} \end{equation}

And the normalized mass and stiffness matrices are:

\begin{equation} \bm{m}_n = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}; \quad \bm{k}_n = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 3 \end{bmatrix} \frac{k}{m} \end{equation}

Note that the diagonal terms of the stiffness matrix are the squares of the corresponding three eigenvalues. The normalized stiffness matrix is known as the spectral matrix.

Normalizing with respect to mass results in an identify principal mass matrix and squares of the eigenvalues on the diagonal in the principal stiffness matrix, this normalization technique is thus very useful for the following reason.

Since we know the form of the principal matrices when normalizing with respect to mass, no multiplying of modal matrices is actually required: the homogeneous principal equations of motion can be written by inspection knowing only the eigenvalues.

Transforming Initial Conditions and Forces

Let's determine how to transform initial conditions and forces to the principal coordinate system. Then, we can solve for transient and forced responses in the principal coordinate system using the uncoupled equations.

We start with the equation of motion in the physical coordinates:

\begin{equation} \bm{m} \ddot{\bm{z}} + \bm{k} \bm{z} = \bm{F} \end{equation}

Pre-multiplying by \(\bm{z}_n^T\) and inserting \(I = \bm{z}_n \bm{z}_n^{-1}\) gives:

\begin{equation} \bm{z}_n^T \bm{m} \bm{z}_n \bm{z}_n^{-1} \ddot{\bm{z}} + \bm{z}_n^T \bm{k} \bm{z}_n \bm{z}_n^{-1} \bm{z} = \bm{z}_n^T \bm{F} \end{equation}

Which is re-written in the following form:

\begin{equation} \bm{m}_p \ddot{\bm{z}}_p + \bm{k}_p \bm{z}_p = \bm{F}_p \end{equation}

where:

\(\bm{m}_p = \bm{z}_n^T \bm{m} \bm{z}_n\) the diagonal principal mass matrix
\(\bm{k}_p = \bm{z}_n^T \bm{k} \bm{z}_n\) the diagonal principal stiffness matrix
\(\ddot{\bm{z}}_p = \bm{z}_n^{-1} \ddot{\bm{z}}\) acceleration vector in principal coordinates
\(\bm{z}_p = \bm{z}_n^{-1} \bm{z}\) displacement vector in principal coordinates
\(\bm{F}_p = \bm{z}_n^T \bm{F}\) force vector in principal coordinates

The vectors of initial displacements \(\bm{z}_{op}\) and velocities \(\dot{\bm{z}}_{op}\) in the principal coordinate system can be expressed as:

\begin{align} \bm{z}_{op} &= \bm{z}_n^{-1} \bm{z}_0 \\ \dot{\bm{z}}_{op} &= \bm{z}_n^{-1} \dot{\bm{z}}_0 \end{align}

where \(\bm{z}_0\) and \(\dot{\bm{z}}_0\) are the vectors of initial displacements and velocities in the physical coordinate system.

Back-Transforming from Principal to Physical Coordinates

We have now everything required to solve the equations in the principal coordinate system.

The variables in physical coordinates are the positions and velocities of the masses. The variables in principal coordinates are the displacements and velocities of each mode of vibration.

The equations in the principal coordinate system can easily be solved, since the equations are uncoupled. We now need to back transform the results in the principal coordinate system to the physical coordinate system to get the final answer.

We showed previously that the relationship between physical and principal coordinates is:

\begin{equation} \bm{z}_n^{-1} \bm{z} = \bm{z}_p \end{equation}

where:

\(\bm{z}_n\) is the matrix containing the eigenvectors (\(n \times m\))
\(\bm{z}\) is the vector of physical coordinates (\(n \times 1\))
\(\bm{z}_p\) is the vector of the principal coordinates (\(m \times 1\))

And thus, the displacement vector in physical coordinates is obtained by pre-multiplying the vector of displacement in principal coordinates by the normalized modal matrix \(\bm{z}_n\):

\begin{equation} \bm{z} = \bm{z}_n \bm{z}_p \end{equation}

Reducing the model size when only selection degrees of freedom are required

The reduction of the number of degrees of freedom is one of the key steps in order to reduce the size of models derived from large finite element simulations.

This can be done as only portions of the eigenvector matrix are needed when only defined dof's have forces applied and other dof's are needed for output.

Let's first examine the force transformation from physical to principal coordinates:

\begin{equation} \bm{F}_p = \bm{z}_n^T \bm{F} = \begin{bmatrix} z_{n11} & z_{n12} & z_{n13} \\ z_{n21} & z_{n22} & z_{n23} \\ z_{n31} & z_{n32} & z_{n33} \end{bmatrix}^T \begin{bmatrix} F_1 \\ F_2 \\ F_3 \end{bmatrix} \end{equation}

If force is only to be applied on mass 1 (\(F_2 = F_3 = 0\)), then only the first row of the modal matrix is required to transform the force in physical coordinates to the force in principal coordinates.

Let's now examine the displacement transformation from principal to physical coordinates:

\begin{equation} \bm{z} = \bm{z}_n \bm{z}_p = \begin{bmatrix} z_{n11} & z_{n12} & z_{n13} \\ z_{n21} & z_{n22} & z_{n23} \\ z_{n31} & z_{n32} & z_{n33} \end{bmatrix} \begin{bmatrix} z_{p1} \\ z_{p2} \\ z_{p3} \end{bmatrix} \end{equation}

And thus, if we are only interested in the physical displacement of the mass 2 (\(z_2 = z_{n21} z_{p1} + z_{n22} z_{p2} + z_{n23} z_{p3}\)), only the second row of the modal matrix is required to transform the three displacements \(z_{p1}\), \(z_{p2}\), \(z_{p3}\) in principal coordinates to \(z_2\).

Only the rows of the modal matrix that correspond to degrees of freedom to which forces are applied and/or for which displacements are desired are required to complete the model.

Damping in Systems with Principal Modes

If a mechanical system is designed with a specific viscous damping element, for example a dashpot, then that element can be added to the system as a viscous damper. The resulting system is linear, but probably does not exhibit normal modes. This leads to the inability to diagonalize and uncouple the equations of motion.

Damping in typical structures arises from hysteresis losses in the materials as they are strained, in some cases from viscous losses due to structure/fluid interaction but more importantly form relative motion at the interfaces and boundaries where different parts are attached or grounded. Unless a specific damping element is used in a structural design, most structures have damping which varies from mode to mode and will be in the range of 0.05% to 2% of critical damping.

Conditions Necessary for Existence of Principal Modes in Damped System

With a conservative (undamped) system, normal modes of vibration will exist. In order to have normal modes in a damped system, the modes shapes must be the same as for the undamped case, and the various parts of the system must pass through their minimum and maximum positions at the same instant in time.

A sufficient condition for the existence of damped normal modes is that the damping matrix be a linear combination of the mass of stiffness matrices:

\begin{equation} \bm{c} = a \bm{m} + b \bm{k} \end{equation}

The damped equations of motion are:

\begin{equation} \bm{m} \ddot{\bm{z}} + \bm{c} \dot{z} + \bm{k} \bm{z} = \bm{F} \end{equation}

And the damping matrix expressed in the principal coordinates is:

\begin{equation} \bm{c}_p = \bm{z}_n^T \bm{c} \bm{z}_n = a \bm{z}_n^T \bm{m} \bm{z}_n + b \bm{z}_n^T \bm{k} \bm{z}_n = a \bm{I} + b \bm{k}_p \end{equation}

We note:

\begin{equation} c_{pi} = a + b \omega_i^2 = 2 \xi_i \omega_i \end{equation}

where \(\xi_i\) is the percentage of critical damping for the i'th mode:

\begin{equation} \xi_i = \frac{c_i}{2 \sqrt{k_{pi} m_{pi}}} = \frac{c_i}{2 m_{pi} \sqrt{\omega_i^2}} = \frac{a + b \omega_i^2}{2 \omega_i} \end{equation}

And the equation in principal coordinates becomes:

\begin{equation} \ddot{z}_{pi} + 2 \xi_i \omega_i \dot{z}_{pi} + \omega_i^2 z_{pi} = F_{pi} \end{equation}

This type of damping is known as proportional damping, where the damping for each mode is proportional to the critical damping for that mode.

Simple Proportional Damping

Viscous damping in each mode is taken to be an arbitrary percentage \(\xi\) of the critical damping \(c_{cr}\):

\begin{equation} c = \frac{\xi}{2 \sqrt{k m}} = \xi \cdot c_{cr} \end{equation}

Proportional to Stiffness Matrix: "Relative" Damping

Recognizing that the higher modes of vibration damp out quickly, "relative" damping yields damping in proportion to frequencies in normal modes:

\begin{equation} \xi_i = \frac{a + b\omega_i^2}{2 \omega_i} \end{equation}

If a value of \(\xi_1\) for the first mode, is assumed:

\begin{equation} b = \frac{2 \xi_1}{\omega_1} \end{equation}

The value of the damping for any other mode is:

\begin{equation} \xi_i = \xi_1 \frac{\omega_i}{\omega_1} \end{equation}

Proportional to Mass Matrix: "Absolute" Damping

Absolute damping is based on making \(b = 0\), in which case the percentage of critical damping is inversely proportional to the natural frequency of each mode.

The procedure to obtain the frequency response from a modal form is as follow:

use the eigenvalues and eigenvectors to define the equations of motion in principal coordinates and to transform forces to principal coordinates
use Laplace transform to obtain the transfer functions in principal coordinates
back-transform the transfer functions to physical coordinates where the individual mode contributions will be evident

This will be applied to the model shown in Figure 9.

{{< figure src="/ox-hugo/hatch00_tdof_model.png" caption="Figure 9: tdof undamped model for modal analysis" >}}

Review from Previous Results

Since the problem we are solving is frequency response, or finding the steady state motion of each mass as a function of frequency and of applied forces, initial conditions are not required.

From previous analysis, we know the eigenvalues and eigenvectors normalized with respect to mass:

\begin{equation} \omega_1 = 0, \quad \omega_2 = \pm\sqrt{\frac{3k}{m}}, \quad \omega_3 = \pm\sqrt{\frac{k}{m}} \end{equation}

\begin{equation} \bm{z}_n = \frac{1}{\sqrt{m}} \begin{bmatrix} \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{2}} & \frac{1}{\sqrt{6}} \\ \frac{1}{\sqrt{3}} & 0 & \frac{-2}{\sqrt{6}} \\ \frac{1}{\sqrt{3}} & \frac{-1}{\sqrt{2}} & \frac{1}{\sqrt{6}} \end{bmatrix} \end{equation}

Knowing that in principal coordinates the mass matrix is the identify matrix and the stiffness matrix is a diagonal matrix with the squares of the respect eigenvalues as terms, we can write the matrices by inspection:

\begin{equation} \bm{m}_n = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}, \quad \bm{k}_n = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 3 \end{bmatrix} \frac{k}{m} \end{equation}

The force vector in principal coordinates is:

\begin{equation} \bm{F}_p = \bm{z}_n^T \bm{F} \end{equation}

The equations of motion in principal coordinates are then:

\begin{equation} \bm{m}_n \ddot{\bm{z}}_p + \bm{k}_n \bm{z}_p = \bm{z}_n^T \bm{F} \end{equation}

which give:

\begin{align} \ddot{z}_{p1} &= (F_1 + F_2 + F_3) \frac{1}{\sqrt{3m}} \\ \ddot{z}_{p2} + \frac{k}{m} z_{p2} &= (F_1 - F_3) \frac{1}{\sqrt{2m}} \\ \ddot{z}_{p3} + \frac{3k}{m} z_{p3} &= (F_1 - 2 F_2 + F_3) \frac{1}{\sqrt{6m}} \end{align}

Transfer Functions - Laplace Transforms in Principal Coordinates

Taking the Laplace transform of each equation gives:

\begin{equation} \begin{bmatrix} \frac{z_{p1}}{F_{1}} \\ \frac{z_{p2}}{F_{1}} \\ \frac{z_{p3}}{F_{1}} \end{bmatrix} = \begin{bmatrix} \frac{1}{s^{2}\sqrt{3m}} \\ \frac{1}{(s^{2} + \omega_{2}^{2})\sqrt{2m}} \\ \frac{1}{(s^{2} + \omega_{3}^{2})\sqrt{6m}} \end{bmatrix} = \begin{bmatrix} z_{p11} \\ z_{p21} \\ z_{p31} \end{bmatrix} \end{equation}

\begin{equation} \begin{bmatrix} \frac{z_{p1}}{F_{2}} \\ \frac{z_{p2}}{F_{2}} \\ \frac{z_{p3}}{F_{2}} \end{bmatrix} = \begin{bmatrix} \frac{1}{s^{2}\sqrt{3m}} \\ 0 \\ \frac{-2}{(s^{2} + \omega_{3}^{2})\sqrt{6m}} \end{bmatrix} = \begin{bmatrix} z_{p12} \\ z_{p22} \\ z_{p32} \end{bmatrix} \end{equation}

\begin{equation} \begin{bmatrix} \frac{z_{p1}}{F_{3}} \\ \frac{z_{p2}}{F_{3}} \\ \frac{z_{p3}}{F_{3}} \end{bmatrix} = \begin{bmatrix} \frac{1}{s^{2}\sqrt{3m}} \\ \frac{-1}{(s^{2} + \omega_{2}^{2})\sqrt{2m}} \\ \frac{1}{(s^{2} + \omega_{3}^{2})\sqrt{6m}} \end{bmatrix} = \begin{bmatrix} z_{p13} \\ z_{p23} \\ z_{p33} \end{bmatrix} \end{equation}

Back-Transforming Mode Contributions to Transfer Functions in Physical Coordinates

The transfer functions in principal coordinates can be back-transformed to physical coordinates. This allows one to see the contributions of each mode, where \(z_{ij}\) is the physical displacement at dof i due to a force at dof j:

\begin{equation} \bm{z} = \bm{z}_n \bm{z}_p \end{equation}

with:

\(\bm{z}\) physical displacements (\(n \times 1\))
\(\bm{z}_n\) matrix of normalized eigenvectors (\(n \times m\))
\(\bm{z}_p\) principal displacements (\(m \times 1\))

And the transfer functions \(\frac{z_i}{F_j}\) can be computed. For instance, the contributions to the transfer function \(\frac{z_1}{F_1}\) are:

\begin{align} \frac{z_1}{F_1} &= \underbrace{z_{n11} z_{p11}}_{\text{1st mode}} + \underbrace{z_{n12} z_{p21}}_{\text{2nd mode}} + \underbrace{z_{n13} z_{p31}}_{\text{3rd mode}} \\ & = \frac{\frac{1}{3m}}{s^2} + \frac{\frac{1}{2m}}{s^2 + \omega_2^2} + \frac{\frac{1}{6m}}{s^2 + \omega_3^2} \end{align}

Forcing Function Combinations to Excite Single Mode

It is instructive to see what types of forcing function combinations will excite each of the three modes separately. From the definition of normal modes, we known that if the system is started from initial displacement conditions that match one of the normal modes, the system will respond at that only mode.

An analogous situation exists for combinations of forcing functions.

The forces transform in the principal coordinates using:

\begin{equation} \bm{F}_p = \bm{z}_n^T \bm{F} \end{equation}

Thus, if \(\bm{F}\) is aligned with \(\bm{z}_{ni}\) (the i'th normalized eigenvector), then \(\bm{F}_p\) will be null except for its i'th term and only the i'th mode will be excited.

How Modes Combine to Create Transfer Functions

Any transfer function derived from the modal analysis is an additive combination of sdof systems.

Each single degree of freedom system has a gain determined by the appropriate eigenvector entries and a resonant frequency given by the appropriate eigenvalue.

It can be shown that for a general system with \(m\) undamped modes:

\begin{equation} \frac{z_j}{F_k} = \sum_{i = 1}^m \frac{z_{nji} z_{nki}}{s^2 + \omega_i^2} \label{eq:general_add_tf} \end{equation}

If modes have some damping:

\begin{equation} \frac{z_j}{F_k} = \sum_{i = 1}^m \frac{z_{nji} z_{nki}}{s^2 + 2 \xi_i \omega_i s + \omega_i^2} \label{eq:general_add_tf_damp} \end{equation}

Equations eq:general_add_tf and eq:general_add_tf_damp shows that in general every transfer function is made up of additive combinations of single degree of freedom systems, with each system having its DC gain determined by the appropriate eigenvector entry product divided by the square of the eigenvalue, \(z_{nji} z_{nki}/\omega_i^2\), and with resonant frequency defined by the eigenvalue \(\omega_i\).

Figure 10 shows the separate contributions of each mode to the total response \(z_1/F_1\).

{{< figure src="/ox-hugo/hatch00_z11_tf.png" caption="Figure 10: Mode contributions to the transfer function from \(F_1\) to \(z_1\)" >}}

The zeros for SISO transfer functions are the roots of the numerator, however, from modal analysis we can see that the zeros arise when modes combine with appropriate phase such that the resulting motion is null.

SISO State Space Matlab Model from ANSYS Model

Introduction

In this section is developed a SISO state space Matlab model from an ANSYS cantilever beam model as shown in Figure 11. A z direction force is applied at the midpoint of the beam and z displacement at the tip is the output. The objective is to provide the smallest Matlab state space model that accurately represents the pertinent dynamics.

{{< figure src="/ox-hugo/hatch00_cantilever_beam.png" caption="Figure 11: Cantilever beam with forcing function at midpoint" >}}

The steps to define the smallest model are:

define the eigenvector elements for all modes for only the input and output degrees of freedom
analyze the modal contributions of all the modes and sort them to define which ones have the greatest contribution

One method for reducing the size of a modal model is to simple truncate the higher frequency modes. However, if truncation is performed without understanding the contributions of each of the modes to the response, several problems could arise. One problem is that deleting a high frequency mode with a significant DC gain could adversely affect the model.

ANSYS Eigenvalue Extraction Methods

ANSYS has a number of different eigenvalue extraction techniques, but for most problems only two methods are commonly used:

The first method, Block Lanczos, is the fastest and calculates all the eigenvalues or eigenvectors in a specific frequency range. This method is usually preferred as most practical models require knowledge of the modes from DC through a specified higher frequency.
The second method, Reduced, performs a Guyan reduction on the model to reduce its size, then calculates all the eigenvalues for the reduced model. Obtaining eigenvector components for the reduced degrees of freedom requires an additional calculation step in ANSYS. This method is still useful to obtain a "super-element".

We can choose to use ANSYS to output only the eigenvectors for nodes of interest or we can output the complete modal matrix and choose the appropriate rows of data within Matlab.

Matlab State Space Model from ANSYS Eigenvalue Run

In order to obtain a State Space model of reasonable order, we need to rank the relative importance of the contributions of each of the individual modes. To do so, we can use a ranking of DC gains.

Once the modes are ranked, the most important can be selected for use, and modes with small DC gains are eliminated from the model. The DC gain contributions of the eliminated modes are not included in the overall DC gain, so there is error in the low frequency gain. In order to eliminate this error, the Matlab function modred is introduced. Using modred is analogous to using Guyan reduction to reduce some less important degrees of freedom.

We will discuss in this section two methods of sorting, one which is applicable for models with the same value of damping for all modes \(\xi_i = \xi\) ("uniform" damping), and another which is applicable for models with different damping values for each mode ("non-uniform" damping).

The general equation for the overall transfer function of undamped and damped systems are:

\begin{align} \frac{z_j}{F_k} &= \sum_{i = 1}^m \frac{z_{nji} z_{nki}}{s^2 + \omega_i^2} \\ \frac{z_j}{F_k} &= \sum_{i = 1}^m \frac{z_{nji} z_{nki}}{s^2 + 2 \xi_i \omega_i s + \omega_i^2} \end{align}

The DC gain of the i'th mode can be obtained by substituting \(s = j\omega = 0\):

\begin{equation} \text{DC gain}_i = \frac{z_{ji}}{F_{ki}} = \frac{z_{nji} z_{nki}}{\omega_i^2} \end{equation}

where \(z_{nji} z_{nki}\) is the product of the j'th (output) row and k'th (force applied) row terms of the i'th eigenvector.

At resonance, the peak gain amplitude of each mode is given by substituting \(s = j \omega_i\):

\begin{equation} \text{Peak gain}_i = \frac{z_{ji}}{F_{ki}} = \frac{-j}{2 \xi_i} \frac{z_{nji} z_{nki}}{\omega_i^2} \end{equation}

And we see that the peak gain for a mode is the DC gain of the same mode divided by \(2 \xi_i\).

If the same value of \(\xi\) is used for all modes, then all the DC gain terms are divided by the same \(2 \xi\) terms and the relative amplitude of the DC gains and peak gains are the same, so there is no difference between sorting a uniform damping model using DC gain or peak gain.

However, if the modes have different damping the relationship between the DC gain and peak gain for all the modes is not a constant value and peak gain must be used to rank modes for importance.

The Matlab command modred can be used for reducing models while retaining the overall system DC gain. The matchdc gain option for the function modred reduces defined states by setting the derivatives of the states to be eliminated to zero, then solving for the remaining states. The method essentially sets up the eliminated states to be "infinitely fast" and is analogous to Guyan reduction in that the low frequency effects of the eliminated states are included in the remaining states. In thus case, modred produces a new state space model with a non-null \(D\) matrix meaning that the reduced system will have a flat response at high frequency which may not be desirable.

The truncate method does not account for the DC gains of the unused modes, which can result in error in the low frequency portion of the frequency response. However, the truncate method has the advantage that it does not exhibit the unusual high frequency direct transmission matrix related behavior of the matchdc method.

If sorting of DC gain values is performed prior to the truncate operation, the system DC gain error may be acceptable while maintaining better high frequency performance.

Ground Acceleration Matlab Model From ANSYS Model

Model Description

Initial ANSYS Model Comparison

Matlab State Space Model from Eigenvalue Run

SISO Disk Drive Actuator Model

In this section we wish to extract a SISO state space model from a Finite Element model representing a Disk Drive Actuator (Figure 12).

Actuator Description

{{< figure src="/ox-hugo/hatch00_disk_drive_siso_model.png" caption="Figure 12: Drawing of Actuator/Suspension system" >}}

The primary motion of the actuator is rotation about the pivot bearing, therefore the final model has the coordinate system transformed from a Cartesian x,y,z coordinate system to a Cylindrical \(r\), \(\theta\) and \(z\) system, with the two origins coincident (Figure 13).

{{< figure src="/ox-hugo/hatch00_disk_drive_nodes_reduced_model.png" caption="Figure 13: Nodes used for reduced Matlab model. Shown with partial finite element mesh at coil" >}}

For reduced models, we only require eigenvector information for dof where forces are applied and where displacements are required. Figure 13 shows the nodes used for the reduced Matlab model. The four nodes 24061, 24066, 24082 and 24087 are located in the center of the coil in the z direction and are used for simulating the VCM force. The arrows at the nodes indicate the direction of forces.

Nodes 22 and 10022 are the nodes for the top and bottom heads and are used for measurement.

Ansys Actuator/Suspension Model Results

A recommended sequence for analyzing dynamic finite element models is:

Plot resonant frequencies versus mode numbers to get a feel for the frequency range. See if there are any significant jumps in frequency between modes which can indicate the system transitioning from one type of characteristic motion to another.
Plot frequency responses to define which modes couple into the response.
Plot and animate the mode shapes that contribute to the response, identifying modes that couple into motions in directions of interest and those that do not. Visually get a sense of how the geometry of the structure affects the modes.
Run parameter studies to understand the sensitivity of critical modes to design variables: dimensions, tolerances, material properties, etc.

Ansys Output Example Listing

A small section of the exported .eig file from ANSYS is shown bellow..

LOAD STEP= 1 SUBSTEP= 1 FREQ= 8.1532 LOAD CASE= 0

THE FOLLOWING DEGREE OF FREEDOM RESULTS ARE IN THE GLOBAL COORDINATE SYSTEM

NODE UX UY UZ ROTX ROTY ROTZ 1 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2 0.0000 1.1981 0.0000 -0.99396E-015-0.14742E-014 0.23368E-001 3 0.31447E-014 0.13758 0.31164E-014 0.33825E-014-0.24915E-015 0.81730E-002 4 -0.51842E-014 0.54120 -0.36031E-014 0.13267E-014 0.30620E-014 0.15962E-001 5 0.0000 4.5604 0.0000 0.50842E-015 0.21315E-015 0.43289E-001

Important information are:

SUBSTEP: mode number
FREQ: eigenvalue in Hz
The eigenvector for the mode in all dofs for all nodes

Inputs and Outputs definition

A unity force is applied at the coil, and evenly distributed among the four nodes. This model is still considered as a Single Input model because the same force is applied to all four coil nodes, requiring only a single column vector for the input matrix \(\bm{B}\).

To compute the DC gain and peak gain in that case, we define a composite forcing function which consists of the force applied to each node times the eigenvector value for that node, \(\bm{F}_p^T \bm{x}_n\). The dimensions of this operation are \((1 \times n) \times (n \times m) = (1 \times m)\), so we have a composite force vector for each mode.

Building Full State Space Matrices

From Ansys, we have the eigenvalues \(\omega_i\) and eigenvectors \(\bm{z}\).

Balanced Reduction

In this chapter another method of reducing models, “balanced reduction”, will be introduced and compared with the DC and peak gain ranking methods.

This method uses the concepts of controllability and observability, commonly referenced in the control community.

One issue with balanced reduction is that we lose the ability to directly identify individual modes in the reduced system model. After balanced reduction one needs to examine the system matrix to identify which modes are included, while the dc and peak gain ranking techniques retain the identities of the individual modes.

Unlike SISO models, which can be easily ranked using simple DC and peak gain techniques, MIMO models will require the balanced reduction method because it easily handles the problem of ranking multiple inputs and outputs.

Reviewing DC gain raking

So far we have used DC or peak gains of the individual modes to rank the importance of including each mode in the reduced system.

For any mode, if the degree of freedom associated with the applied force has a zero value, then the force applied at the degree of freedom cannot excite that mode, so the DC and peak gains will also be zero and the mode be eliminated.

Similarly, if the degree of freedom associated with the output has a zero value, then no matter how much force is applied to that mode, there will be no output, and the mode can be eliminated.

A mode which cannot be excited by the applied force is said to be uncontrollable, and a mode which has no output in the desired direction is said to be unobservable.

Controllability, Observability

For a state space system described by:

\begin{align*} \dot{\bm{x}} &= \bm{A} \bm{x} + \bm{B} u \\ \bm{y} &= \bm{C} \bm{x} \end{align*}

the following definitions of controllability hold:

If there is an input \(u\) that can move the system from some arbitrary state \(\bm{x}_1\) to another arbitrary state \(\bm{x}_2\) in a finite time, then the system is controllable
A controllability matrix \(\bm{\mathcal{C}}\) can be formed as:

\begin{equation} \bm{\mathcal{C}} = \begin{bmatrix} \bm{B} & \bm{A} \bm{B} & \bm{A}^{2} \bm{B} & \dots & \bm{A}^{n-1} \bm{B} \end{bmatrix} \end{equation}

If \(\bm{\mathcal{C}}\) has full (row) rank, the system is controllable.
Another definition of controllability involves the controllability gramian, \(\bm{W}_c\), the solution of the Lyapunov equation:

\begin{equation} \bm{A} \bm{W}_{c} + \bm{W}_{c} \bm{A}^{T} + \bm{B} \bm{B}^{T} = 0 \end{equation}

defined as:

\begin{equation} \bm{W}_{c} = \int_{0}^{\infty} e^{\bm{A}\tau} \bm{B} \bm{B}^{T} e^{{\bm{A}^{T} \tau} d \tau} \end{equation}

If the solution \(\bm{W}_c\) is non-singular, then the system is controllable.

Diagonal elements of the controllability gramian give information about the relative controllability of the different modes and can be used in a manner similar to our use of DC gains to rank the relative controllability of individual modes.

Gramians exists only for system that have all their poles with negative real parts. Thus, if a system has rigid body modes, it should first be partitioned into rigid body dynamics and oscillatory dynamics and then the controllability gramian can be computed on the oscillatory partition.

A similar set of definitions can be made for observability:

If the initial state \(\bm{x}_0\) of a system can be inferred from knowledge of the input \(u\) and output \(\bm{y}\) over a finite period of time \((0, t)\), then the system is said to be observable.
A observability matrix \(\bm{\mathcal{O}}\) can be formed as:

\begin{equation} \bm{\mathcal{O}} = \begin{bmatrix} \bm{C} \\ \bm{C} \bm{A} \\ \bm{C} \bm{A}^{2} \\ \vdots \\ \bm{C} \bm{A}^{n-1} \end{bmatrix} \end{equation}

If \(\bm{\mathcal{O}}\) has full (column) rank, the system is observable.
Another definition of observability involves the observability gramian, \(\bm{W}_o\), the solution of the Lyapunov equation:

\begin{equation} \bm{A}^T \bm{W}_{o} + \bm{W}_{o} \bm{A} + \bm{C}^T \bm{C} = 0 \end{equation}

defined as:

\begin{equation} \bm{W}_{o} = \int_{0}^{\infty} e^{\bm{A}^T \tau} \bm{C}^T \bm{C} e^{{\bm{A} \tau} d \tau} \end{equation}

If the solution \(\bm{W}_o\) is non-singular, then the system is observable.

Diagonal elements of the observability gramian give information about the relative observability of the different modes and can be used in a manner similar to our use of DC gains to rank the relative observability of individual modes.

Ranking Using Controllability/Observability

We could use the controllability curve to rank the states for controllability and eliminate those states with low controllability. Alternately, we could use the observability curve to rank the states for observability and then eliminate states with low observability. The problem with this approach is that the joint controllability/observability is not taken in account.

To take into account both the observability and the controllability, we use the Balance Reduction.

Balanced Reduction

The balanced reduction is used to create a system with identical diagonal controllability and observability Gramians. Since the two gramians are equal, either the diagonal of the controllability gramian or the diagonal of the observability gramian can be used to rank states for elimination.

The diagonal terms of the joint gramian are squares of the Hankel singular values of the system. The Hankel matrix is the product of the controllability and observability gramians.

Because the controllability and observability gramians are identical, there is no ambiguity in deciding whether the most controllable or the most observable states should be chosen. The states to be kept are the states with the largest diagonal terms.

MIMO Two Stage Actuator Model

In this section, a MIMO two-stage actuator model is derived from a finite element model (Figure 14).

Actuator Description

{{< figure src="/ox-hugo/hatch00_disk_drive_mimo_schematic.png" caption="Figure 14: Drawing of actuator/suspension system" >}}

A piezo-actuator is now bounded into one side of each of the arms. The piezo actuator consists of a ceramic element that changes size when a voltage is applied.

Then the fine positioning motion of the piezo is used in conjunction with VCM's coarse positioning motion, higher servo bandwidth is possible.

Instead of applying voltage as the input into the piezo elements, we will assume that we have calculated an equivalent set of forces which can be applied at the ends of the element that will replicate the voltage force function. In this model, we will be applying forces to multiple nodes at the ends of both piezo elements.

Since the same forces are being applied to both piezo elements, they represent the second input to the MIMO system, the first input being the coil force.

Ansys Model Description

In Figure 15 are shown the principal nodes used for the model.

{{< figure src="/ox-hugo/hatch00_disk_drive_mimo_ansys.png" caption="Figure 15: Nodes used for reduced Matlab model, shown with partial mesh at coil and piezo element" >}}

Matlab Model

For a SISO system, we can rank the relative importance of modes using two methods, by using DC or peak gains and by using balancing. However, for a MIMO system, balancing is the only practical solution.

If we were using DC gains to rank the modes, we would have to compute DC gains for the four combinations possible of the two inputs and two inputs.

We would then see that the DC gains of the modes are not ranked to same way for different sets of input/output. Especially, the modes important for the VCM are not the same than for the piezo. If one were to choose a single ranking for the model which would take into account both inputs and both outputs, it is difficult to see how to do it given the DC gain rankings. Thus the necessity of balanced reduction for MIMO models.

Balancing Reduction

Balancing the system involves calculating gramians, which are only defined for negative definite systems. This requires separating the rigid body mode from the oscillatory modes and balancing the oscillatory modes. The system matrices are partitioned and a model of only oscillatory modes is created and balanced. Plotting the diagonal gramian terms reveals the relative important of the states.

The complete system is rebuilt by augmenting the rigid body mode with the reduced oscillatory modes.

Matlab Simple Example

Problem Description

We define the system parameters.

  m = 1;
  k = 1;

We write the mass and stiffness matrices:

  M = diag([m, m, m]);
  K = [k, -k,   0;
      -k, 2*k, -k;
       0, -k,   k];

Compute the eigenvalues and eigenvectors:

  [z, w] = eig(M\K);

	rad/s
w1	0.0
w2	1.0
w3	3.0

Normalization of the eigenvectors:

  zn = z./sqrt(diag(z' * M * z));

zn1	zn2	zn3
-0.6	-0.7	0.4
-0.6	0.0	-0.8
-0.6	0.7	0.4

Non-necessary step:

  Mn = zn' * M * zn;
  Kn = zn' * K * zn;

By inspection:

  Mn = eye(3);
  Kn = w; % Shouldn't this be equal to w.^2 ?

We add some simple proportional damping:

  xi = 0.03;
  Cn = xi/2*sqrt(k*m) * eye(3);

The equations in the principal coordinates are: \[ \bm{m}_n \ddot{\bm{z}}_p + \bm{c}_n \dot{\bm{z}}_p + \bm{k}_n \bm{z}_p = \bm{F}_p = \bm{z}_n^T \bm{F} \]

Let's note \[ \bm{G}_p(s) = \frac{\bm{z}_p}{\bm{F}} \]

  Gp = (tf(Mn)*s^2 + tf(Cn)*s + tf(Kn))\tf(eye(3))*zn';

  bodeFig({Gp(1,1), Gp(2,2), Gp(3,3)})

And we have the Laplace transform in the principal coordinates.

And to convert the modal displacement to the physical displacement, we use: \[ \bm{z} = \bm{z}_n \bm{z}_p \] And we note: \[ \bm{G}(s) = \bm{z}_n\ bm{G}_p(s) = \frac{\bm{z}}{\bm{F}} \]

  G = zn * Gp;

{{< figure src="/ox-hugo/hatch00_z13_tf.png" caption="Figure 16: Mode contributions to the transfer function from \(F_1\) to \(z_3\)" >}}

{{< figure src="/ox-hugo/hatch00_z11_tf.png" caption="Figure 17: Mode contributions to the transfer function from \(F_1\) to \(z_1\)" >}}

Matlab with ANSYS

Extract values

  filename = 'files/cantbeam30bl.eig';

  dir = 3; % UY
  [xn, f0] = readEigFile(filename, dir);

  n_nodes = size(xn, 1);
  n_modes = size(xn, 2);

Define Physical Inputs and Outputs

First, define the node numbers corresponding to the inputs and outputs

  i_input = 14;
  i_output = 29;

Define Damping

We here use uniform damping.

  xi = 0.01;

Alternative Definition of the matrix A

I could define 2x2 sub-matrices each corresponding to a particular mode and then combine them using blkdiag command.

All Modes Included in the Model

System Matrix - A

  Adiag = zeros(2*n_modes,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0);

  Adiagsup = zeros(2*n_modes-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*n_modes-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*n_modes, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*n_modes);
  C(:, 1:2:end) = xn(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_f = ss(A, B, C, D);

Simple mode truncation

Let's plot the frequency of the modes (Figure 18).

{{< figure src="/ox-hugo/hatch00_cant_beam_modes_freq.png" caption="Figure 18: Frequency of the modes" >}}

{{< figure src="/ox-hugo/hatch00_cant_beam_unsorted_dc_gains.png" caption="Figure 19: Unsorted DC Gains" >}}

Let's keep only the first 10 modes.

  m_max = 10;
  xn_t = xn(:, 1:m_max);
  f0_t = f0(1:m_max);

  Adiag = zeros(2*m_max,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0_t);

  Adiagsup = zeros(2*m_max-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*m_max-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0_t).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*m_max, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn_t'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*m_max);
  C(:, 1:2:end) = xn_t(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_t = ss(A, B, C, D);

Modes sorted by their DC gain

Let's sort the modes by their DC gains and plot their sorted DC gains.

  dc_gain = abs(xn(i_input, :).*xn(i_output, :))./(2*pi*f0).^2;

  [dc_gain_sort, index_sort] = sort(dc_gain, 'descend');

{{< figure src="/ox-hugo/hatch00_cant_beam_sorted_dc_gains.png" caption="Figure 20: Sorted DC Gains" >}}

Let's keep only the first 10 sorted modes.

  m_max = 10;

  xn_s = xn(:, index_sort(1:m_max));
  f0_s = f0(index_sort(1:m_max));

  Adiag = zeros(2*m_max,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0_s);

  Adiagsup = zeros(2*m_max-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*m_max-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0_s).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*m_max, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn_s'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*m_max);
  C(:, 1:2:end) = xn_s(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_s = ss(A, B, C, D);

Comparison

  freqs = logspace(0, 5, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_f, freqs, 'Hz'))), 'DisplayName', 'Full');
  plot(freqs, abs(squeeze(freqresp(G_t, freqs, 'Hz'))), 'DisplayName', 'Trun');
  plot(freqs, abs(squeeze(freqresp(G_s, freqs, 'Hz'))), 'DisplayName', 'Sort');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend();

Effect of the Individual Modes

  freqs = logspace(0, 4, 1000);

  figure;
  hold on;
  for mode_i = 1:6
      A = zeros(2);
      A(2,2) = -2*xi.*(2*pi*f0(mode_i));
      A(1,2) = 1;
      A(2,1) = -(2*pi*f0(mode_i)).^2;

      B = [0; xn(i_input, mode_i)'];

      C = [xn(i_output, mode_i), 0];

      D = zeros(length(i_output), length(i_input));

      plot(freqs, abs(squeeze(freqresp(ss(A,B,C,D), freqs, 'Hz'))), ...
           'DisplayName', sprintf('Mode %i', mode_i));
  end
  plot(freqs, abs(squeeze(freqresp(G_f, freqs, 'Hz'))), 'k--', ...
       'DisplayName', 'Full');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend();

Non-Uniform Damping

Definition of the Damping

If we want to use Rayleigh damping:

  a = 1e-2;
  b = 1e-6;
  xi = (a + b * (2*pi*f0).^2)./(2*pi*f0);

System Creation

System Matrix - A

  Adiag = zeros(2*n_modes,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0);

  Adiagsup = zeros(2*n_modes-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*n_modes-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*n_modes, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*n_modes);
  C(:, 1:2:end) = xn(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_d = ss(A, B, C, D);

Comparison with Uniform Damping

  freqs = logspace(0, 5, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_f, freqs, 'Hz'))), 'DisplayName', 'Uniform Damping');
  plot(freqs, abs(squeeze(freqresp(G_d, freqs, 'Hz'))), 'DisplayName', 'Non-Uniform Damping');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend();

Modes sorted by their peak gain

Let's sort the modes by their peak gains and plot their sorted peak gains.

  dc_gain = abs(xn(i_input, :).*xn(i_output, :))./(2*pi*f0).^2;
  peak_gain = dc_gain./xi;

  [peak_gain_sort, index_sort] = sort(peak_gain, 'descend');

Let's keep only the first 10 sorted modes.

  m_max = 10;

  xn_s = xn(:, index_sort(1:m_max));
  f0_s = f0(index_sort(1:m_max));
  xi_x = xi(index_sort(1:m_max));

  Adiag = zeros(2*m_max,1);
  Adiag(2:2:end) = -2*xi_s.*(2*pi*f0_s);

  Adiagsup = zeros(2*m_max-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*m_max-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0_s).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*m_max, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn_s'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*m_max);
  C(:, 1:2:end) = xn_s(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_p = ss(A, B, C, D);

  freqs = logspace(0, 5, 1000);

  figure;
  hold on;
  plot(freqs, abs(squeeze(freqresp(G_f, freqs, 'Hz'))), 'DisplayName', 'Uniform Damping');
  plot(freqs, abs(squeeze(freqresp(G_d, freqs, 'Hz'))), 'DisplayName', 'Non-Uniform Damping');
  plot(freqs, abs(squeeze(freqresp(G_p, freqs, 'Hz'))), 'DisplayName', 'Peak sort');
  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
  ylabel('Amplitude'); xlabel('Frequency [Hz]');
  legend();

MIMO System

Inputs and Outputs

Let's choose two inputs and two outputs.

  i_input = [14, 31];
  i_output = [14, 31];

System Matrices

System Matrix - A

  Adiag = zeros(2*n_modes,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0);

  Adiagsup = zeros(2*n_modes-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*n_modes-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*n_modes, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*n_modes);
  C(:, 1:2:end) = xn(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_m = ss(A, B, C, D);

Balancing Reduction

First, we have to make sure that the rigid body mode is not included in the system (here it is not).

Then, we compute the controllability and observability gramians.

  wc = gram(G_m, 'c');
  wo = gram(G_m, 'o');

And we plot the diagonal terms

{{< figure src="/ox-hugo/hatch00_gramians.png" caption="Figure 21: Observability and Controllability Gramians" >}}

We use balreal to rank oscillatory states.

[SYSB,G] = BALREAL(SYS) computes a balanced state-space realization for the stable portion of the linear system SYS. For stable systems, SYSB is an equivalent realization for which the controllability and observability Gramians are equal and diagonal, their diagonal entries forming the vector G of Hankel singular values. Small entries in G indicate states that can be removed to simplify the model (use MODRED to reduce the model order).

  [G_b, G, T, Ti] = balreal(G_m);

{{< figure src="/ox-hugo/hatch00_cant_beam_gramian_balanced.png" caption="Figure 22: Sorted values of the Gramian of the balanced realization" >}}

Now we can choose the number of states to keep.

  n_states_b = 20;

We now use modred to define reduced order oscillatory system using mathdc or truncate option.

MODRED Model simplification by state elimination.

RSYS = MODRED(SYS,ELIM) simplifies the state-space model SYS by discarding the states specified in the vector ELIM. The full state vector X is partitioned as X = [X1;X2] where Xr=X1 is the reduced state vector and X2 is discarded.

  G_br = modred(G_b, n_states_b+1:size(A,1), 'truncate');

If needed, the rigid body mode should be added to the reduced system.

And other option is to specify the minimum value of the gramians diagonal elements for the modes to keep.

  G_min = 1e-4;
  G_br = modred(G_b, G<G_min, 'truncate');

  figure;
  bodemag(G_br, G_m)

Flexible Bodies with Simscape

Extract values

  filename = 'files/cantbeam30bl.eig';

  [xn, f0] = extractEigs(filename, 'dirs', [0 0 1 0 0 0]);

  n_nodes = size(xn, 1);
  n_modes = size(xn, 2);

Read Mass and Stiffness matrices from Ansys

  str = fileread('files/Kdense.txt');
  % Remove spaces
  str = regexprep(str,'\s+','');
  % Regex to get the data
  parts = regexp(str, '\[(?<row>\d+),(?<col>\d+)\]:(?<val>[^\[]+)', 'names');

  row = cellfun(@str2double, {parts.row}, 'UniformOutput', true);
  col = cellfun(@str2double, {parts.col}, 'UniformOutput', true);
  val = cellfun(@str2double, {parts.val}, 'UniformOutput', true);

  sz = [max(row), max(col)];  % size of output matrix
  MatK = zeros(sz);  % preallocate size
  ix = sub2ind(sz, row, col);  % get matrix positions
  MatK(ix)= val;   % assign data

  str = fileread('files/Mdense.txt');
  % Remove spaces
  str = regexprep(str,'\s+','');
  % Regex to get the data
  parts = regexp(str, '\[(?<row>\d+),(?<col>\d+)\]:(?<val>[^\[]+)', 'names');

  row = cellfun(@str2double, {parts.row}, 'UniformOutput', true);
  col = cellfun(@str2double, {parts.col}, 'UniformOutput', true);
  val = cellfun(@str2double, {parts.val}, 'UniformOutput', true);

  sz = [max(row), max(col)];  % size of output matrix
  MatM = zeros(sz);  % preallocate size
  ix = sub2ind(sz, row, col);  % get matrix positions
  MatM(ix)= val;   % assign data

Find inputs/outputs:

  i_input = 14;
  i_output = 29;

Correspondence with DOF:

  a = readtable('files/Mass_HB.mapping', 'FileType', 'text');
  KM_i = strcmpi(a{:, 3},{'UZ'});

  MatM = MatM(KM_i, KM_i);
  MatK = MatK(KM_i, KM_i);

Read Position of Nodes

  a = readmatrix('files/FEA-poutre-noeuds.txt');
  pos = a(:, 4:6);
  node_i = a(:, 7);

  figure;
  hold on;
  plot3(pos(:,1),pos(:,3),pos(:,3), 'ko')
  text(pos(:,1),pos(:,3),pos(:,3), num2cell(node_i))

Reduced order model

Define Physical Inputs and Outputs

  i_input = 14;
  i_output = 29;

Damping

  xi = 0.01;

  dc_gain = abs(xn(i_input, :).*xn(i_output, :))./(2*pi*f0).^2;

  [dc_gain_sort, index_sort] = sort(dc_gain, 'descend');

  m_max = 13;

  xn_s = xn(:, index_sort(1:m_max));
  f0_s = f0(index_sort(1:m_max));

  Adiag = zeros(2*m_max,1);
  Adiag(2:2:end) = -2*xi.*(2*pi*f0_s);

  Adiagsup = zeros(2*m_max-1,1);
  Adiagsup(1:2:end) = 1;

  Adiaginf = zeros(2*m_max-1,1);
  Adiaginf(1:2:end) = -(2*pi*f0_s).^2;

  A = diag(Adiag) + diag(Adiagsup, 1) + diag(Adiaginf, -1);

System Matrix - B

  B = zeros(2*m_max, length(i_input));

  for i = 1:length(i_input)
      % Physical Coordinates
      Fp = zeros(n_nodes, 1);
      Fp(i_input(i)) = 1;

      B(2:2:end, i) = xn_s'*Fp;
  end

System Matrix - C

  C = zeros(length(i_output), 2*m_max);
  C(:, 1:2:end) = xn_s(i_output, :);

System Matrix - D

  D = zeros(length(i_output), length(i_input));

State Space Model

  G_s = ss(A, B, C, D);

Generate Mass and Stiffness matrices

Full Mass and Stiffness matrices in the principal coordinates:

  Mp = eye(length(f0));
  Kp = xn'*diag((2*pi*f0).^2)/xn;

Reduced Mass and Stiffness matrices in the principal coordinates:

  Mr = Mp()
  Kr = xn'*diag((2*pi*f0).^2)/xn;

Reduced Mass and Stiffness matrices in the physical coordinates:

  M = xn_s*eye(m_max)/xn_s;
  K = xn_s*diag((2*pi*f0_s).^2)/xn_s;

  % M = xn*eye(length(f0))/xn;
  % K = xn*diag((2*pi*f0).^2)/xn;

  M = eye(length(f0));
  K = xn*diag((2*pi*f0).^2)/xn;

Frames for Simscape

  pos_frames = pos([1, i_input, i_output], :);

Bibliography

Hatch, Michael R. 2000. Vibration Simulation Using Matlab and Ansys. CRC Press.

Miu, Denny K. 1993. Mechatronics: Electromechanics and Contromechanics. 1st ed. Mechanical Engineering Series. Springer-Verlag New York.

74 KiB Raw Blame History

Introduction

Modal Analysis

Model Size Reduction

Notations

Zeros in SISO Mechanical Systems

State Space Analysis

Modal Analysis

Eigenvalue Problem

Equation of Motion

Principal Mode Definition

Eigenvalues / Characteristic Equation

Eigenvectors

Modal Matrix

Uncoupling the Equations of Motion

Normalizing Eigenvectors

Normalizing with Respect to Unity

Normalizing with Respect to Mass

Transforming Initial Conditions and Forces

Back-Transforming from Principal to Physical Coordinates

Reducing the model size when only selection degrees of freedom are required

Damping in Systems with Principal Modes

Conditions Necessary for Existence of Principal Modes in Damped System

Simple Proportional Damping

Proportional to Stiffness Matrix: "Relative" Damping

Proportional to Mass Matrix: "Absolute" Damping

Modal Analysis: State Space Form

Frequency Response: Modal Form

Review from Previous Results

Transfer Functions - Laplace Transforms in Principal Coordinates

Back-Transforming Mode Contributions to Transfer Functions in Physical Coordinates

Forcing Function Combinations to Excite Single Mode

How Modes Combine to Create Transfer Functions

SISO State Space Matlab Model from ANSYS Model

Introduction

ANSYS Eigenvalue Extraction Methods

Matlab State Space Model from ANSYS Eigenvalue Run

Ground Acceleration Matlab Model From ANSYS Model

Model Description

Initial ANSYS Model Comparison

Matlab State Space Model from Eigenvalue Run

SISO Disk Drive Actuator Model

Actuator Description

Ansys Actuator/Suspension Model Results

Ansys Output Example Listing

Inputs and Outputs definition

Building Full State Space Matrices

Balanced Reduction

Reviewing DC gain raking

Controllability, Observability

Ranking Using Controllability/Observability

Balanced Reduction

MIMO Two Stage Actuator Model

Actuator Description

Ansys Model Description

Matlab Model

Balancing Reduction

Matlab Simple Example

Problem Description

Matlab with ANSYS

Extract values

Define Physical Inputs and Outputs

Define Damping

Alternative Definition of the matrix A

All Modes Included in the Model

Simple mode truncation

Modes sorted by their DC gain

Comparison

Effect of the Individual Modes

Non-Uniform Damping

Definition of the Damping

System Creation

Comparison with Uniform Damping

Modes sorted by their peak gain

MIMO System

Inputs and Outputs

System Matrices

Balancing Reduction

Flexible Bodies with Simscape

Extract values

74 KiB

Raw Blame History