+++ 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](https://in.mathworks.com/matlabcentral/fileexchange/2186-vibration-simulation-using-matlab-and-ansys). ## Introduction {#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. ### Modal Analysis {#modal-analysis} The diagram in Figure [1](#figure--fig:hatch00-modal-analysis-flowchart) shows the methodology for analyzing a lightly damped structure using normal modes.
The steps are: 1. deriving the undamped equations of motion in physical coordinates 2. solving the eigenvalue problem, yielding eigenvalues (natural frequencies) and eigenvectors (mode shapes). These first 2 steps are usually performed using a Finite Element software 3. 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. 4. proportional damping is applied 5. 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 {#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](#figure--fig:hatch00-model-reduction-chart) 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 {#notations} Tables [3](#figure--fig:hatch00-n-dof-zeros), [2](#table--tab:notations-eigen-vectors-values) and [3](#table--tab:notations-stiffness-mass) 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 {#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](#figure--fig:hatch00-n-dof-zeros) 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 {#state-space-analysis} ## Modal Analysis {#modal-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: 1. **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 2. **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 3. **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 4. **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 {#eigenvalue-problem} #### Equation of Motion {#equation-of-motion} Let's consider the model shown in Figure [4](#figure--fig:hatch00-undamped-tdof-model) 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 {#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 {#eigenvalues-characteristic-equation} Re-injecting normal modes into the equation of motion 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 {#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](#figure--fig:hatch00-tdof-mode-1), [6](#figure--fig:hatch00-tdof-mode-2) and [7](#figure--fig:hatch00-tdof-mode-3). {{< 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" >}} #### Modal Matrix {#modal-matrix} The modal matrix is an \\(n \times m\\) matrix with columns corresponding to the \\(m\\) system eigenvectors as shown in Eq. \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 {#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--fig:hatch00-schematic-modal-solution). {{< 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 {#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 {#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 {#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 {#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 {#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 {#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 {#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 {#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 {#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 {#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 {#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. ## Modal Analysis: State Space Form {#modal-analysis-state-space-form} ## Frequency Response: Modal Form {#frequency-response-modal-form} 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--fig:hatch00-tdof-model). {{< figure src="/ox-hugo/hatch00_tdof_model.png" caption="Figure 9: tdof undamped model for modal analysis" >}} ### Review from Previous Results {#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 {#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 {#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 {#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 {#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 and 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](#figure--fig:hatch00-z11-tf-example) 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 {#siso-state-space-matlab-model-from-ansys-model} ### Introduction {#introduction} In this section is developed a SISO state space Matlab model from an ANSYS cantilever beam model as shown in Figure [11](#figure--fig:hatch00-cantilever-beam). 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: 1. define the eigenvector elements for all modes for only the input and output degrees of freedom 2. 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-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 {#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 {#ground-acceleration-matlab-model-from-ansys-model} ### Model Description {#model-description} ### Initial ANSYS Model Comparison {#initial-ansys-model-comparison} ### Matlab State Space Model from Eigenvalue Run {#matlab-state-space-model-from-eigenvalue-run} ## SISO Disk Drive Actuator Model {#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](#figure--fig:hatch00-disk-drive-siso-model)). ### Actuator Description {#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--fig:hatch00-disk-drive-nodes-reduced-model)). {{< 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](#figure--fig:hatch00-disk-drive-nodes-reduced-model) 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 {#ansys-actuator-suspension-model-results} A recommended sequence for analyzing dynamic finite element models is: 1. 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. 2. Plot frequency responses to define which modes couple into the response. 3. 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. 4. Run parameter studies to understand the sensitivity of critical modes to design variables: dimensions, tolerances, material properties, etc. ### Ansys Output Example Listing {#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 {#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 {#building-full-state-space-matrices} From Ansys, we have the eigenvalues \\(\omega\_i\\) and eigenvectors \\(\bm{z}\\). ## Balanced Reduction {#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 {#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 {#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: 1. 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 2. 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. 3. 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: 1. 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. 2. 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. 3. 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 {#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 {#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 {#mimo-two-stage-actuator-model} In this section, a MIMO two-stage actuator model is derived from a finite element model (Figure [14](#figure--fig:hatch00-disk-drive-mimo-schematic)). ### Actuator Description {#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 {#ansys-model-description} In Figure [15](#figure--fig:hatch00-disk-drive-mimo-ansys) 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 {#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-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 {#matlab-simple-example} ### Problem Description {#problem-description} We define the system parameters. ```matlab m = 1; k = 1; ``` We write the mass and stiffness matrices: ```matlab M = diag([m, m, m]); K = [k, -k, 0; -k, 2*k, -k; 0, -k, k]; ``` Compute the eigenvalues and eigenvectors: ```matlab [z, w] = eig(M\K); ``` | | rad/s | |----|-------| | w1 | 0.0 | | w2 | 1.0 | | w3 | 3.0 | Normalization of the eigenvectors: ```matlab 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: ```matlab Mn = zn' * M * zn; Kn = zn' * K * zn; ``` By inspection: ```matlab Mn = eye(3); Kn = w; % Shouldn't this be equal to w.^2 ? ``` We add some simple proportional damping: ```matlab 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}} \\] ```matlab Gp = (tf(Mn)*s^2 + tf(Cn)*s + tf(Kn))\tf(eye(3))*zn'; ``` ```matlab 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}} \\] ```matlab 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 {#matlab-with-ansys} ### Extract values {#extract-values} ```matlab 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 {#define-physical-inputs-and-outputs} First, define the node numbers corresponding to the inputs and outputs ```matlab i_input = 14; i_output = 29; ``` ### Define Damping {#define-damping} We here use uniform damping. ```matlab xi = 0.01; ``` ### Alternative Definition of the matrix A {#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 {#all-modes-included-in-the-model} System Matrix - A ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*n_modes); C(:, 1:2:end) = xn(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_f = ss(A, B, C, D); ``` ### Simple mode truncation {#simple-mode-truncation} Let's plot the frequency of the modes (Figure [18](#figure--fig:hatch00-cant-beam-modes-freq)). {{< 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. ```matlab m_max = 10; xn_t = xn(:, 1:m_max); f0_t = f0(1:m_max); ``` ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*m_max); C(:, 1:2:end) = xn_t(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_t = ss(A, B, C, D); ``` ### Modes sorted by their DC gain {#modes-sorted-by-their-dc-gain} Let's sort the modes by their DC gains and plot their sorted DC gains. ```matlab 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. ```matlab m_max = 10; xn_s = xn(:, index_sort(1:m_max)); f0_s = f0(index_sort(1:m_max)); ``` ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*m_max); C(:, 1:2:end) = xn_s(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_s = ss(A, B, C, D); ``` ### Comparison {#comparison} ```matlab 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 {#effect-of-the-individual-modes} ```matlab 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 {#non-uniform-damping} #### Definition of the Damping {#definition-of-the-damping} If we want to use Rayleigh damping: ```matlab a = 1e-2; b = 1e-6; xi = (a + b * (2*pi*f0).^2)./(2*pi*f0); ``` #### System Creation {#system-creation} System Matrix - A ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*n_modes); C(:, 1:2:end) = xn(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_d = ss(A, B, C, D); ``` #### Comparison with Uniform Damping {#comparison-with-uniform-damping} ```matlab 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 {#modes-sorted-by-their-peak-gain} Let's sort the modes by their peak gains and plot their sorted peak gains. ```matlab 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. ```matlab 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)); ``` ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*m_max); C(:, 1:2:end) = xn_s(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_p = ss(A, B, C, D); ``` ```matlab 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 {#mimo-system} #### Inputs and Outputs {#inputs-and-outputs} Let's choose two inputs and two outputs. ```matlab i_input = [14, 31]; i_output = [14, 31]; ``` #### System Matrices {#system-matrices} System Matrix - A ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*n_modes); C(:, 1:2:end) = xn(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_m = ss(A, B, C, D); ``` #### Balancing Reduction {#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. ```matlab 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). ```matlab [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. ```matlab 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. ```matlab 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. ```matlab G_min = 1e-4; G_br = modred(G_b, G\d+),(?\d+)\]:(?[^\[]+)', '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 ``` ```matlab str = fileread('files/Mdense.txt'); % Remove spaces str = regexprep(str,'\s+',''); % Regex to get the data parts = regexp(str, '\[(?\d+),(?\d+)\]:(?[^\[]+)', '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: ```matlab i_input = 14; i_output = 29; ``` Correspondence with DOF: ```matlab 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 {#read-position-of-nodes} ```matlab a = readmatrix('files/FEA-poutre-noeuds.txt'); pos = a(:, 4:6); node_i = a(:, 7); ``` ```matlab figure; hold on; plot3(pos(:,1),pos(:,3),pos(:,3), 'ko') text(pos(:,1),pos(:,3),pos(:,3), num2cell(node_i)) ``` ### Reduced order model {#reduced-order-model} Define Physical Inputs and Outputs ```matlab i_input = 14; i_output = 29; ``` Damping ```matlab xi = 0.01; ``` ```matlab dc_gain = abs(xn(i_input, :).*xn(i_output, :))./(2*pi*f0).^2; [dc_gain_sort, index_sort] = sort(dc_gain, 'descend'); ``` ```matlab m_max = 13; xn_s = xn(:, index_sort(1:m_max)); f0_s = f0(index_sort(1:m_max)); ``` ```matlab 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 ```matlab 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 ```matlab C = zeros(length(i_output), 2*m_max); C(:, 1:2:end) = xn_s(i_output, :); ``` System Matrix - D ```matlab D = zeros(length(i_output), length(i_input)); ``` State Space Model ```matlab G_s = ss(A, B, C, D); ``` ### Generate Mass and Stiffness matrices {#generate-mass-and-stiffness-matrices} Full Mass and Stiffness matrices in the principal coordinates: ```matlab Mp = eye(length(f0)); Kp = xn'*diag((2*pi*f0).^2)/xn; ``` Reduced Mass and Stiffness matrices in the principal coordinates: ```matlab Mr = Mp() Kr = xn'*diag((2*pi*f0).^2)/xn; ``` Reduced Mass and Stiffness matrices in the physical coordinates: ```matlab ``` ```matlab M = xn_s*eye(m_max)/xn_s; K = xn_s*diag((2*pi*f0_s).^2)/xn_s; ``` ```matlab % 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 {#frames-for-simscape} ```matlab pos_frames = pos([1, i_input, i_output], :); ``` ## Bibliography {#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.