Parameters: Linear Solver

Model ElementParam_Linear defines the solution control parameters for a linear analysis. These parameters control the types of linear analyses to be done and the output options.

Description

Two types of linear analyses are available in MotionSolve:

Eigenvalue and eigenvector calculations for the system, and,
State matrix calculation (E, A, B, C, and D matrices) for the linearized system, when a set of inputs and outputs are defined.

Linear analysis can be performed at any operating point. This may be done either after a static equilibrium solution or after a dynamic solution. Usually, a steady state configuration is chosen as the operating condition. Not all the MotionSolve modeling elements are supported for a linear analysis. For more details, please see the Comments section.

Format

<Param_Linear
   [       
        anim_scale              = "real"
        balancing               = {"TRUE" | "FALSE" | "AUTO"}            
        disable_damping         = {"TRUE" | "FALSE" }        
        gyroscopic              = {"NONE"|"PARTIAL"|"FULL"}
        eigen_analysis          = {"TRUE"|"FALSE"}
        write_energy_dist       = {"TRUE"|"FALSE"}
        write_energy_html       = {"TRUE"|"FALSE"} 
        write_eig_info          = {"TRUE"|"FALSE"}
        write_simulinkmdl       = {"TRUE"|"FALSE"}
        write_matlabfiles       = {"TRUE"|"FALSE"} 
        write_oml               = {"TRUE"|"FALSE"}
        pinput_id               = "integer"
        poutput_id              = "integer"
        pstate_id               = "integer"
        numerical_Jacobian      = {"TRUE" | "FALSE"}
       {
         ke_modes_exclude       = "string"
         ke_modes_include       = "string"
        }
        {
         se_modes_exclude       = "string"
         se_modes_include       = "string"
        }
        {
         de_modes_exclude       = "string"
         de_modes_include       = "string"
        }
   ]
/>

Attributes

anim_scale: Specifies a scale factor for magnifying the mode shapes during animation.
The default is 1.0.

balancing

Specifies whether the A matrix should be pre-conditioned using diagonal scaling to improve robustness of the eigenvalue solution. Choose from:

TRUE: MotionSolve will always balance the A matrix using diagonal scaling.

FALSE: MotionSolve will not balance the A matrix.

AUTO: MotionSolve determines when to balance the A matrix based on the condition number of the eigenvector matrix. See Comment 8 for more details.

The default for balancing is AUTO.

disable_damping

Specifies whether the linearization solver should disable damping from all force elements for the eigenvalue solution.

The default is FALSE; damping will be considered for the eigenvalue solution.

gyroscopic

Specifies how MotionSolve should compute gyroscopic effects when linearizing a rotating system. Choose from:

NONE: MotionSolve will not consider gyroscopic effects during linear analysis.
PARTIAL: MotionSolve will consider partial gyroscopic effect without accounting for angular momentum.
FULL: MotionSolve will consider the full gyroscopic effect including angular momentum.

The default is PARTIAL.

eigen_analysis

Specifies whether linear analysis should perform or skip the eigen analysis. This can be used in cases when only the A, B, C, and D matrix are needed. If NO is selected, then eigen information and energy distribution information will not be written. Select from TRUE or FALSE.

The default is TRUE.

write_energy_dist

Specifies whether the modal kinetic strain and dissipative energy distribution is written out to the solver log file and the *_linz.mrf output file. Select from TRUE or FALSE.

The default is FALSE.

See Comment 7 for more details.

The default is FALSE.

write_energy_html

Specifies whether the modal kinetic strain and dissipative energy distribution is written out in an HTML format that can be displayed by any modern browsers. Select from TRUE or FALSE.

The default is TRUE.

write_eig_info

Specifies whether the eigenvalue and eigenvector data are written to an .eig file. Select from TRUE or FALSE.

The default is TRUE.

write_simulinkmdl

Specifies whether the A, B, C, and D matrices that are calculated are to be written out in the Simulink MDL format or not. Select from TRUE or FALSE.

The default is TRUE.

write_matlabfiles

Specifies whether the A, B, C, and D matrices that are calculated are to be written out to a file that can be read by MATLAB. Select from TRUE or FALSE.

The default is TRUE.

write_oml

Specifies whether the A, B, C, and D matrices that are calculated are to be written out to a file that can be read by Compose or any other OML-based math tool, such as GNU Octave. Select from TRUE or FALSE.

The default is TRUE.

pinput_id

Specifies the plant input ID used for the B and D state matrices. Can be optionally used with the write_simulinkmdl, write_matlabfile, and/or write_oml option.

poutput_id

Specifies the plant output ID used for the C and D state matrices. Can be optionally used with the write_simulinkmdl, write_matlabfile, and/or write_oml option.

pstate_id

An integer that specifies the id of the Reference_PlantState that defines some or all of the states MotionSolve used for linearizing the equations of motion. When the pstate_id attribute is not specified, MotionSolve uses its own internal states to linearize the equations.

numerical_Jacobian

Defines the process of the linearization. If FALSE, each element of the system is being linearized first before the entire system is assembled. If TRUE, the entire system is being assembled first and then linearized. The latter option is computationally more expensive and can take significantly more time for larger models.

The default is FALSE.

ke_modes_exclude, ke_modes_include, se_modes_exclude, se_modes_include, de_modes_exclude, de_modes_include

These attributes allow you to exclude or include modes from the following energy distribution tables.

Kinetic energy distribution table -- ke_modes
Strain energy distribution table -- se_modes
Damping "energy" distribution table -- de_modes

You may exclude or include modes using any one of the options below. ke_modes are used for the example, but you can use se_modes or de_modes also with these options.

A single mode

ke_modes_exclude = "1"
ke_modes_include = "1"

Multiple modes

ke_modes_exclude = "1, 2, 3"
ke_modes_include = "1, 2, 3"

Range of modes

ke_modes_exclude = "1:3"
ke_modes_include = "1:3"

A combination of the above

ke_modes_exclude = "1, 2, 3:4, 5:7"
ke_modes_include = "1, 2, 3:4, 5:7"

All modes

ke_modes_exclude = "ALL"
ke_modes_include = "ALL"

None of the modes

ke_modes_exclude = "NONE"
ke_modes_include = "NONE"

The default for ke_modes_exclude is "NONE".

The default for ke_modes_include is "ALL".

The default for se_modes_exclude is "NONE".

The default for se_modes_include is "ALL".

The default for de_modes_exclude is "NONE".

The default for de_modes_include is "ALL".

Comments

Eigenvector and Eigenvalue Analysis:
The eigenvalues and eigenvectors of a nonlinear system can change with time. An eigenvalue analysis is useful for understanding the underlying vibration characteristics of the system at a specific operating point or a series of operating points.

The eigenvalues represent the complex frequencies of vibration of the system and the eigenvectors represent the modes of vibration.
State Matrix Analysis:
Given a set of inputs "u", a set of outputs "y", and a linearized mechanical system represented by a set of dynamical states "x", the transfer function for the system can be represented in the time domain as:
$E \dot{x} = A x + B u$

$y = C x + D u$
The matrices E, A, B, C, and D are called the state matrices for the system. If there are nx states, nu inputs, and ny outputs, the dimensions of the matrices are as follows:
- E: nx x nx
- A: nx x nx
- B: nx x nu
- C: nu x nx
- D: nu x nu
The inputs "u" are defined using the Control_PlantInput modeling element specified in the pinput_id and the outputs "y" are defined using Control_PlantOutput modeling element specified in the poutput_id.

The inputs "u" can be displacements, forces, or moments. If the inputs are forces or moments, then matrix E becomes an identity matrix and the transfer function for the system can be simplified to:
$\begin{array}{l} \dot{x} = A x + B u \\ y = C x + D u \end{array}$

State matrix output is particularly useful for control engineers, who need a linear plant representation to design a controller.

Another common application of state matrices is to calculate the system Frequency Response Functions to a series of time-based inputs. This is the traditional vibration analysis of linearized mechanical systems.
State matrix output:
- When the write_simulinkmdl attribute is chosen, the A, B, C, and D matrices are written in Simulink format in a Simulink MDL file.
- When the write_matlabfiles attribute is chosen, the A, B, C, and D matrices are written in MATLAB format in four separate files with extensions .a, .b, .c and .d, respectively.
- When the write_oml attribute is chosen, the A, B, C, and D matrices are written in OML format in one file with extension .oml.
Eigenvalue analysis output:
- Eigenvalues are written to the *.eig file.
- To visualize the mode shapes, set the linear_animation = “TRUE” attribute in the H3DOutput command. Each mode shape will be written as new subcase in the H3D file.
Significance of eigenvalues and eigenvectors:
Eigenvalues and eigenvectors can be computed from state matrices. Once a linear analysis is initiated, the system is linearized to the form:

$\begin{array}{l} M \ddot{q} + C \dot{q} + E_{1} \dot{z} + E_{2} z = 0 \\ E_{3} \dot{q} + E_{4} q + \dot{z} + E_{5} z = 0 \end{array}$

q represents the linearized displacement coordinates, and z represents non-mechanical states due to other differential equations.

All constraints are eliminated from the formulation. The equations are converted to first-order form as follows:

$[\begin{matrix} M & 0 & E_{1} \\ 0 & I & 0 \\ 0 & E_{3} & I \end{matrix}] {\begin{matrix} \dot{u} \\ \dot{q} \\ \dot{z} \end{matrix}} + [\begin{matrix} C & K & E_{2} \\ - I & 0 & 0 \\ 0 & E_{4} & E_{5} \end{matrix}] {\begin{matrix} u \\ q \\ z \end{matrix}} = {\begin{matrix} 0 \\ 0 \\ 0 \end{matrix}}$

This results in the eigenvalue problem:
$\hat{K} x = λ \hat{M} x$ where $x^{T} = [\begin{matrix} u^{T} & q^{T} & z^{T} \end{matrix}]$ , $\hat{M} = [\begin{matrix} M & 0 & E_{1} \\ 0 & I & 0 \\ 0 & E_{3} & I \end{matrix}]$ , and $\hat{K} = [\begin{matrix} C & K & E_{2} \\ - I & 0 & 0 \\ 0 & E_{4} & E_{5} \end{matrix}]$
- Let an Eigenvalue be represented as: $λ = λ_{r} + j λ_{i}; j = \sqrt{- 1}$
  The undamped natural frequency of the this mode is: $ω_{n} = \sqrt{(λ_{r}^{2} + λ_{i}^{2})}$
  
  The damping ratio for this mode is: $ζ = \frac{λ_{r}}{ω_{n}}$
  
  The damped natural frequency for this mode is: $ω_{d} = ω_{n} \sqrt{(1 + ζ^{2})}$
- The real part of the Eigenvalue corresponds to the system damping. For all stable systems, this must be less than zero. A positive real component for an eigenvalue indicates instability in the system. The system becomes unstable if this mode is ever excited. Such designs are to be avoided.
- By successfully varying design parameters like stiffness and damping in the system and calculating the eigenvalues, a root locus plot of the system at an operating point can be generated.
- The eigenvectors can be used to validate the analytical model with actual physical tests that may have been performed in the frequency domain.
Significance of state matrices:
Take the Laplace transform of the linearized equations of motion:

$L (\begin{array}{l} \dot{x} = A x + B u \\ y = C x + D u \end{array}) = (\begin{array}{l} s X (s) - X (0) = A X (s) + B U (s) \\ Y (s) = C X (s) + D u (s) \end{array})$

Ignoring the initial value X(0), and rearranging terms, you obtain:

System Transfer Function = $\frac{Y (s)}{U (s)} = C {(s I - A)}^{- 1} B + D$
The meaning of the transfer function can change based on its input and output types:
- When the input "u" are displacements and the output "y" forces, the matrices E-A-B-C-D together represent the dynamic compliance at the operating point
- When the input "u" are Forces and the output "y" displacements, the matrices A-B-C-D together (E is identity) represent the dynamic stiffness at the operating point
- When the input "u" are Forces and the output "y" velocities, the matrices A-B-C-D together (E is identity) represent the mechanical impedance at the operating point
- When the input "u" are Forces and the output "y" accelerations, the matrices A-B-C-D together (E is identity) represent the apparent mass at the operating point
Compliance matrices are used in the automotive industry for suspension design. A compliance matrix is a MIMO transfer function that relates a set of forces and torques applied at specific points on the suspension to the displacements measured at other specific points on the suspension.

Traditionally, compliance matrices have been computed about a static equilibrium position and the data used for suspension design. This approach ignores three effects: (a) The effects of inertia and damping on the system compliance, (b) the variation of compliance on during vehicle operation, and (c) The variation of compliance with the frequency of the input.

The use of ABCD matrices to calculate the compliance transfer function about dynamic operating points overcomes all of these issues.

Computing the energy distribution in the system:

If write_energy_dist is set to TRUE, MotionSolve computes the kinetic, strain, and dissipative energy distributions for parts and forces in the model for the modes specified.

If write_energy_dist is set to FALSE, MotionSolve will not compute any energy distribution regardless of the ke_modes_include/ke_modes_exclude, se_modes_include/se_modes_exclude, or de_modes_include/de_modes_exclude attributes and their values.

You can specify the modes that MotionSolve should compute the energy distributions for. This can be useful when there are a large number of modes in your system and you are interested only in the energy distribution for a subset of these modes.

You can specify the modes MotionSolve should include or exclude, separately, for kinetic, strain, or dissipative energy by using the following attributes:

ke_modes_include – for including modes in the kinetic energy distribution
ke_modes_exclude – for excluding modes in the kinetic energy distribution
se_modes_include – for including modes in the strain energy distribution
se_modes_exclude – for excluding modes in the strain energy distribution
de_modes_include – for including certain modes in the dissipative energy distribution
de_modes_exclude – for excluding certain modes in the dissipative energy distribution

Typically, you only need to specify the modes to be included, or the modes to be excluded, but not both. If you specify both - modes to be included and modes to be excluded - the definition that appears later in the Param_Linear statement will be honored.

Kinetic Energy:

MotionSolve computes the modal kinetic energy distribution for all bodies in the MBD model, for the modes specified. These are then written to the log file in tabular format. Each row in this table represents a body. Within a row, each number represents the percentage distribution of modal kinetic energy between the translational (X, Y, Z), rotational (RXX, RYY, RZZ) and cross-rotational, RXY, RXZ, RYZ directions of the part. The total modal kinetic energy distributions for each mode should add up to 100%.

The total kinetic energy for each mode is calculated as follows:

K E = (Φ^{T} M Φ) ω_{n}^{2}

$Φ$ is the eigenvector

$M$ is the mass matrix of size n_DOF*n_Bodies X n_DOF*n_Bodies. n_DOF is the number of degrees of freedom per body = 6 and n_Bodies is the total number of bodies in the system.

$ω_{n}$ is the undamped frequency of the mode.

The modal kinetic energy distribution for each part is calculated for nine coordinates: X, Y, Z, RX, RY, RZ, RXY, RXZ, and RYZ.

For example, the percentage of total KE for the Z direction is calculated as:

% K E_{z} = (\frac{p_{1}^{T} M p_{1}}{p_{1}^{T} M p}) \times 100

$% K E_{z}$ is percentage of modal kinetic energy in the Z direction.

$M$ is the mass matrix as described above.

$p$ is a column of the eigenvector matrix Φ corresponding to the current mode number. If mode number = 10, then this is the 10^th column vector in the eigenvector matrix. This is a vector of size 6*n_Bodies X 1.

$p_{1}$ is the same vector as $p$ with all entries zeroed-out except for the coordinate of interest, the Z coordinate, in this example.

For cross rotational terms (RXY, RXZ, RYZ), the calculation is slightly different. For example, for the RXZ percentage:

$% K E_{R x z} = (\frac{p_{1}^{T} M p_{2}}{p_{1}^{T} M p}) \times 100$

$% K E_{R x z}$ is the percentage of modal kinetic energy in the Rxz direction.

$M$ is the mass matrix as described above.

$p$ is a column of the eigenvector matrix Φ corresponding to the current mode number.

$p_{1}$ is the same vector as $p$ with all entries zeroed-out except for the Rx coordinate.

$p_{2}$ is the same vector as $p$ with all entries zeroed out except for the Rz coordinate.

Strain and Dissipative Energy:

MotionSolve can also compute the modal strain and dissipative energy distribution for certain force entities in the MBD model, for the modes specified (these are listed in the table below).

The energy distributions are then written to the log file in tabular format. Each row in this table represents a force entity. Within a row, each number represents the percentage distribution of modal strain or dissipative energy between the translational (X, Y, Z) and rotational (RX, RY, RZ) directions.

Model Element	Total	X	Y	Z	RX	RY	RZ
Force_Beam	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Force_Bushing	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Force_Field	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Force_Vector_OneBody, Force_Vector_TwoBody type = ForceAndTorque	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Force_Vector_OneBody, Force_Vector_TwoBody type = ForceOnly	Yes	Yes	Yes	Yes	No	No	No
Force_Vector_OneBody, Force_Vector_TwoBody type = TorqueOnly	Yes	No	No	No	Yes	Yes	Yes
Force_Scalar_TwoBody type = Force	Yes	Yes	Yes	Yes	No	No	No
Force_Scalar_TwoBody type = Torque	Yes	No	No	No	No	No	No
ForceSpringDamper type=TRANSLATION	Yes	Yes	Yes	Yes	Yes	Yes	Yes
Force_SpringDamper type = ROTATION	No	No	No	No	No	No	No

Note: If a force entity does not contribute to the modal strain or dissipative energy, then that force entity’s row will be zero.

The Strain and Dissipative energy calculations are done in a manner similar to the Kinetic energy distributions except that the stiffness and damping matrices are used in place of the mass matrix to calculate the percentage distributions for force elements.

In addition to the energy distribution tables, MotionSolve also displays a header for each mode that includes information about the mode number, damping ratio, undamped natural frequency, and the total value of the modal kinetic energy. An example of this energy distribution that is written to the output log file as shown below:

*************************
Mode number           =     10
Damping ratio         =    5.0000000E-02
Undamped natural freq.=    1.5915494E+00
Kinetic energy        =    1.1427840E-05

                 Percentage distribution of Kinetic energy
              |    X      Y      Z      RXX    RYY    RZZ    RXY    RXZ    RYZ
              +----------------------------------------------------------------
 PART/30301   |   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00   0.00
 PART/30302   |   0.00  98.69   0.00   1.31   0.00   0.00   0.00   0.00   0.00
              +----------------------------------------------------------------
 
               Percentage distribution of Strain energy
              |  Total    X      Y      Z     RXX    RYY    RZZ
              +-------------------------------------------------
 VFOR/30301   |   0.00   0.00   0.00   0.00
 VFOR/30302   |   0.00   0.00   0.00   0.00
 SPDP/303001  |   0.00   0.00   0.00   0.00   0.00   0.00   0.00
 SPDP/303002  |  50.00  50.00   0.00   0.00   0.00   0.00   0.00
 SPDP/303003  |  50.00  50.00   0.00   0.00   0.00   0.00   0.00
              +-------------------------------------------------
 
               Percentage distribution of Dissipative energy
              |  Total    X      Y      Z     RXX    RYY    RZZ
              +-------------------------------------------------
 VFOR/30301   |   0.00   0.00   0.00   0.00
 VFOR/30302   |   0.00   0.00   0.00   0.00
 SPDP/303001  |   0.00   0.00   0.00   0.00   0.00   0.00   0.00
 SPDP/303002  |  50.00  50.00   0.00   0.00   0.00   0.00   0.00
 SPDP/303003  |  50.00  50.00   0.00   0.00   0.00   0.00   0.00
                  +-------------------------------------------------

These energy distributions are additionally written to the output *_linz.mrf file, which may be used for post-processing.

Note: The *_linz.mrf contains information for all modes regardless of which were specified for inclusion or exclusion.

Matrix Balancing
For certain models, the eigenvalues computed by MotionSolve can be extremely sensitive to small perturbations in the $A$ matrix. This occurs because the eigenvector matrix, $V$ , for such a model is highly ill-conditioned. As a result, even seemingly minor changes for such models may result in large changes in the eigenvalue solution.

$K (V) = {‖ V ‖}^{- 1} ‖ V ‖$ is the condition number of a matrix $V$ . Consider the linear equation $V x = b$ . $K (V)$ measures the rate at which the solution x will change due to a change in b. If $K (V)$ is large, then even a small change in b can cause large changes in the solution x. In such as case, the matrix V is said to be ill-conditioned.

To obtain a robust eigenvalue solution in such cases, MotionSolve must balance the $A$ matrix by finding a diagonal similarity transformation $D A D^{- 1}$ such that the row and column norms of $A$ are numerically close. By doing so, the eigenvector matrix becomes better conditioned and the eigenvalue solution is more robust to numerical noise in $A$ .
The balancing attribute controls whether matrix balancing is done or not.
- Set this to TRUE to always balance the A matrix regardless of whether the eigenvector matrix is ill-conditioned or not.
- Set balancing to FALSE to never balance the A matrix.
- The default for balancing is AUTO, which lets MotionSolve decide if balancing is needed or not.
  Let $R = K {(V)}^{- 1}$ be the reciprocal of the condition number of $V$ .
  
  If $R < 500 * ε_{m a c h i n e}$ , then balancing is performed. Otherwise, it is not performed.
  
  If $ε_{m a c h i n e} = 2.2 E - 16$ , then balancing is performed when $R < 500 * 2.2 E - 16 = 1.11 E - 13$ .
While balancing a matrix is generally a good idea, you need to be judicious in its use. There are some rare situations where balancing will not provide the benefits you expect. Specifically, you should be aware of the following:
- If a matrix contains small elements that are due to round-off error, balancing might make them more significant than they should be. MotionSolve tries to minimize this by zeroing-out very small entries in the Jacobian.
- The condition number of the eigenvector matrix, $K (V)$ , can only be computed after $V$ is computed. When balancing is set to AUTO, the eigenvalue solver may be called twice: first, without balancing and then, if $K (V)$ is small, with balancing. As a result, an increase in run-time may occasionally be seen.
Element Support
MotionSolve does not support all modeling elements for a linear analysis. The following elements are currently not supported:
- Body_Flexible (NLFE bodies only)
- Constraint_General
- Force_Contact
- Force_Penalty
- Reference_2DCluster
- Reference_Variable (Implicit only)
- Subsystem_Planar/Reference 2DCluster
Simulation sequence
The linear analysis step should always be the last step in a sequence of analysis steps. The following sequences are valid:
- 1:Linear
- 1:(Quasi-) Static, 2:Linear
- 1:Transient, 2:Linear
- 1:Static, 2:Transient, 3: Linear
- Other combinations
You can use the Save and Load_Model command statements around a linear analysis step, if further analyses are desired after the linearization. For example:
- 1: Transient
- 2: Save
- 3: Linear
- 4: Reload model from step 2
Cosimulation
The linear analysis does not include elements with which you can perform a cosimulation.