COMSOL 6.3 - The Modal Solver Algorithm

The Modal Solver Algorithm

The purpose of the Modal Solver and the Model Reduction study step is to speed up certain simulations by projecting the solution of the original problem onto a space spanned by some basis vectors, so the solution of the underlying system of equations is approximated by a linear combination of the basis vectors. Instead of solving for the original high-dimensional problem, the algorithm solves for the coefficients of the linear combination, which is low-dimensional.

The underlying system can be well represented by some dominant eigenvectors, so making use of the solution to an eigenvalue or eigenfrequency problem to construct a basis becomes a natural choice. Optionally, this basis can be extended with constraint modes. Each constraint mode is a solution to a stationary problem with a nonhomogeneous boundary condition. These modes make it possible to extend the validity of the reduced model. They also make it possible to use reduced-model inputs in constraints for more general use in a reduced-order model produced by the model reduction.

The Modal solver and Modal Reduction solver use the same algorithm to do reduction. However, the Modal solver cannot handle the case when there is reduced-model inputs. It can only be used to for homogeneous problems or nonhomogeneous problems with constant constraints. The Modal Reduction solver can handle the case when there is reduced-model inputs in constraints.

In the following, the algorithm for the general case, which includes inputs in constraints, is illustrated.

Time Domain

For time-dependent problems, the starting point of modal reduction is the matrix form of the FEM discretization, which is formulated in COMSOL as

(20-10)

where U is the dependent variable, μ represents control inputs, NF is the constant constraint force matrix, and

is the vector of Lagrange multipliers.

Linearizing around a solution U0, μ0, and

, and introducing

(20-11)

Let V(t) = U(t) − U0,

, and

. The linearized system can then be written as

(20-12)

where E is the mass matrix, D is the damping matrix, and K is the stiffness matrix. B, Bdot, and Bdotdot are the input matrix, the time derivative input matrix, and the second time derivative input matrix, respectively. Either E or D can be identically zero. N is the constraint Jacobian matrix. l0 is a constant-in-time vector, and l(t) is a scalar load factor. m0 is a constant constraint vector for the linearized problem.

This linear problem has a solution that is a sum of three parts:

(20-13)

where Vh is a solution satisfying homogeneous boundary conditions, Vc is a solution satisfying the nonhomogeneous boundary conditions that are controlled by control inputs, and Vd is a solution satisfying the nonhomogeneous constant boundary conditions that are not controlled by control inputs. The scheme uses basis vectors that can be used to approximate the solutions Vh and Vc. The vector Vd fulfills nonhomogeneous constraints but where υ = 0; that is,

(20-14)

Construction of Eigenmodes

The modal approach approximates the homogeneous solution Vh by a linear combination of eigenmodes,

(20-15)

where qh are the reduced DOFs to be determined, and each column of Φ is an eigenmode.

Construction of Constraint Modes

Here, constraint modes are defined as a special subset of constraint modes as defined in the original version of the Craig–Bampton method (see Ref[7]) from 1968. In this paper, a constraint mode is defined as “the mode shapes of interior freedoms due to successive unit displacement of boundary points, all other boundary points being totally constrained.”

In general, “boundary points” refer to all DOFs on constraint boundaries, and each one of these DOFs defines a constraint mode with the Craig–Bampton method. These modes together with “fixed normal modes,” which corresponds to the homogeneous eigenvalue problem formulated above for the same constraint formulation, constitute one standard way of connecting substructures, where each substructure is a reduced-order model. The advantage with this construction is that the reduced stiffness matrix is block diagonal, with no connection between the boundary and the interior points. The coupling comes only through the mass and damping matrices. This method is also called static condensation for this reason.

The definition of constraint modes here will be slightly different. Because only global model inputs are supported in constraints, it is only necessary to represent constraints that can be described by these global model inputs. Computing a mode for each boundary DOF would therefore only be an extra cost.

When there are reduced-model inputs in constraints, constraint modes can be introduced to control the constraints set by the inputs. Constraint modes are defined from the equations

(20-16)

where the input ν for the mode j is defined by

(20-17)

and δ is the Kronecker delta. Notice that the matrix N is the same for all constraint modes, which means that they fulfill homogeneous constraints for all boundary conditions that do not depend on υj.

Let Ψ be the matrix whose columns are formed by the constraint modes, and qc be the reduced constraint mode DOFs. Define Vc as

(20-18)

By construction, the boundary conditions are fulfilled as long as

(20-19)

The COMSOL Multiphysics software supports two ways to compute constraint modes: (a) via the sensitivity solver and (b) via the parametric solver. The order in which these modes are solved for, in method (a) or (b), dictates the order in which they are appended to the matrix Ψ in the Modal Solver or Modal Reduction solver, and thereby control the order of the reduced constraint mode DOFs qc.

(a) Sensitivity Solver

The constraint solution V in Equation 20-16 is a linear function of the control variables. It is therefore possible to use the forward sensitivity solver, with the reduced-model inputs in the constraints as sensitivity control variables, and directly compute

(20-20)

Each column in Ψ can be a constraint mode, and qc can be the reduced constraint mode DOFs.

(b) Parametric Solver

Each constraint mode can be constructed directly as a solution to the stationary problem (see Equation 20-16); that is, only one reduced-model input is one ("1"), while the others are zero ("0"), according to Equation 20-17. A simple way to generate such solutions is to introduce a parameter for each constraint and then enable an auxiliary sweep for a Stationary study and these parameters and where all other loads and constraints are zero. By using the specified combinations sweep option, you can solve for the constraint modes by setting the parameters to zero or one according to the condition in Equation 20-17. This construction of constraint modes is similar to the classical Craig–Bampton method (see Ref. 46), but where the boundary DOFs are grouped together. The method here is therefore not as general as the original method, but on the other hand it gives a much smaller reduced system.

Computation of Reduced-System Initial Values and Further Expansion of Basis

For some problems, using eigenmodes and constraint modes is not enough to give a good approximation of the unreduced system. It is then possible to also include initial value and initial derivatives to the basis.

To add initial data to the reduced basis, the algorithm

First solves for the initialization problem of the unreduced time-dependent study under some training value of the inputs, and for the initial time t0. This gives an initial value vector V0 (and the first time-derivative vector).

Next, the method removes the constraint mode solution and Vd from the initial values. For the constraint mode solution, the training expressions are used for the model inputs.

Next, the remaining vector is projected onto the space spanned by eigenmodes. This gives the initial values for the reduced vector qh (and

The remaining vector after this step, is in a vector space that is orthogonal to the space spanned by the used eigenmodes. If this vector is nonvanishing, it makes sense to expand the basis with a normalized version of this vector.

In the following, Ni is used to represent the matrix whose columns consist of the basis from the initial values, and qi is used to represent the corresponding reduced DOFs.

Reduction

Assume that P is the projection matrix whose columns consist of basis vectors. In the case that the basis vectors include eigenmodes, constraint modes, and initial value basis, then

and

Substituting V = Pq to the linearized Equation 20-12 and left multiply it by PH, you will get the following reduced equation (an approximation of the original problem):

(20-21)

where

(20-22)

Additional Damping for the Transient Reduced System

The damping matrix D can be present when performing the eigenvalue analysis. It is, however, possible to add additional damping by providing damping ratios per mode (or one ratio for all modes). If λi denotes the i:th eigenvalue and ξi the associated damping ratio, and Dra represents the additional damping matrix, then

(20-23)

For the modal solver, the initial values and initial time derivatives of DOFs are obtained from the initialization process mentioned in Computation of Reduced-System Initial Values and Further Expansion of Basis.

(20-24)

For the modal solver, which supports inputs in constraints, the contribution from the initial value of control inputs needs to be subtracted. The initial condition becomes

(20-25)

where U0 and Udot0 are q0 and qdot0 in Equation 20-24.

Reconstruction

After solving the reduced system, the solution to the original problem is reconstructed as

(20-26)

Linearized Output for Modal Reduction Solver

The Model Reduction study makes it possible to define outputs variables. For linearized output, they can be defined as

(20-27)

where Y0 is the output bias vector, Cr is the reduced output matrix, and F is the input feedback matrix.

The Second-Order Time-Dependent Problem and State-Space Form

Defining x =

, the second-order system in the reduced form can be rewritten as the first-order system below:

(20-28)

Rewrite Equation 20-28 in state-space form as

(20-29)

where

(20-30)

(20-31)

Initial Conditions for Solving the State-Space System

For the modal reduction solver, the initial value to solve the state-space system is

For the modal reduction solver that supports inputs in constraints, the initial condition is

(20-32)

where

(20-33)

Linearized Output for the Modal Reduction Solver in State-Space Form

The linearized output in state-space form in the case of Y0 = 0 can be written as

(20-34)

where

is the reduced output matrix.

Frequency Domain

Equation 20-12 is linear and cannot be used to describe this case. The starting point for the reduction of a frequency-domain problem is a FEM residual vector expansion around the angular frequency ω0, where ω0 = iλ0, λ0 = −2π i f0, and f0 in the first frequency in the frequency list set in the unreduced frequency domain study. The equation can be written as

(20-35)

Expand the unreduced residual vector L in Equation 20-35 around ω0 in a three-term Taylor expansion:

(20-36)

where

(20-37)

Next, considering a linearization point U = U0, μ = μ0 for the residual vector L(ω, U, μ) in Equation 20-36, and substituting the result to Equation 20-35, you obtain

(20-38)

where V = U − U0, ν= μ − μ0 and

(20-39)

(20-40)

Similarly to the time-domain case, for the frequency-domain case, Ll(ω) is introduced as a load vector for the linearization problem (Equation 20-36). Furthermore, a factorization is supported:

where l0(ω) is a frequency-dependent vector and l(ω) is a scalar load factor.

In the COMSOL Multiphysics software, the matrices K, D, and E in Equation 20-39 and the matrices B, Bdot, and Bdotdot are assembled with respect to λ0, so the matrices differ by a factor of (−i )n, where n can be 0, 1, or 2.

Reduction

This process is similar to the transient case. Suppose the projection matrix is P and V = P∗q. Substituting to Equation 20-38 and multiplying it by PH, you get the following reduced equation for the frequency-domain case:

(20-41)

where

(20-42)

and the definitions of Kr, Dr, Er, Br, Brdot, and Brdotdot are the same as in Equation 20-22.

Inhomogeneous Boundary Conditions

For inhomogeneous Dirichlet boundary conditions that are not controlled by inputs, a particular solution is computed by the modal solver. A particular solution Vd is computed from the unreduced equation

(20-43)

Where the matrix K is the one defined in Equation 20-39.

The term

(20-44)

is then subtracted from the right side of the reduced Equation 20-41.

Additional Damping for the Frequency-Domain Reduced System

(20-45)

and the term iωDra is added to the left hand side of Equation 20-41.

Reconstruction and linearized output for the Modal Reduction solver in the frequency domain is the same as in the time domain.

Exporting the Reduced System

•

Reduced matrices that can be exported from both the Modal solver and the Modal Reduction solver for both time and frequency domain: the projection matrix P; each column of P is a basis vector.

•

Reduced matrices that can be exported from the Modal solver for both the time domain and the frequency domain: the stiffness matrix Kr, the damping matrix Dr, the damping ratio matrix Dra, the mass matrix Er, the input matrix Br, the time derivative input matrix Brdot, and the second time derivative input matrix Brdotdot. Those matrices are defined in Equation 20-22. The output matrix Cr and the input feedback matrix F defined in Equation 20-27 can also be exported. For the Modal solver, Cr and F are expected to be zero because there are no output variables defined.

For time-domain problems, the initial value input matrix B0, the initial value time derivative input matrix B0rdot, and the state-space matrices MA, MB, and Mc defined in Equation 20-30 can also be exported, as well as the output matrix C and the input feedback matrix F. State-space matrices can also be evaluated even for first-order systems with what each matrix represents listed in the following table:

Table 20-7: First-order and Second-order Definitions
Name	Description	Definition for First-Order Case	Definition for second-Order Case
Mc	Mass matrix	Dr	[0 Er; Ir 0]
MA	Stiffness matrix	-Kr	[-Kr -Dr; 0 Ir]
MB	Input matrix	-Br	[-Br; 0]
x0	Initial value vector	U0	[U0; Udot0]
C	Output matrix	Cr	[Cr 0]
D	Input feedback matrix	F	F

Note that D in the table is different from the damping matrix D in the unreduced equation.

•

Reduced matrices that can be exported from the Modal Reduction solver for both the time domain and the frequency domain: All matrices mentioned for Modal solver above can be exported. However, they might not be the full size of the reduced matrix. Rather, it is the submatrix of the full-sized reduced matrix defined in Equation 20-22, which corresponds to all unconstrained reduced DOFs according to the current setting in the Constraint Modes section in a reduced-order model node. The counterparts of those submatrices, which correspond to all constrained reduced DOFs, are denoted as, for example, Krc, Drc, Erc, MAc, MBc, and Mcc. They can also be exported.

•

Reduced vectors can be exported from both the Modal solver and the Modal Reduction solver.

For the time domain, the reduced load vector l0r (it is L in the reduced matrices list) and Kud defined in Equation 20-22, the initial value vector U0, and the initial derivative vector Udot0. For second-order problems, the initial vector x0 for the state-space system can also be exported.

For the frequency domain, the reduced load vector l0r(ω) (it is L in the reduced matrices list) defined in Equation 20-42 for the last given angular frequency ω, the mass matrix times the particular solution PHEVd and the damping matrix times the particular solution PHDVd in Equation 20-44. See also the example below.

•

Vectors that can be exported from the Modal Reduction solver only: constraint modes to inputs map CImap, unconstrained DOF indices to the mode index map DofToMode, and constrained DOF index to that mode index map.

•

Vectors that can be exported from the Modal solver only: All load vectors , which export all load vectors (that is, lor(ω0), lor(ω1),... , lor(ωk)), where ωi corresponds to fi in the setting in the l study.

For the frequency domain, to linearize around an angular frequency ω0 is useful, for example, for problems where the underlying problem has a material that depends on frequency. However, it can give some surprising effects compared to the exported matrices from the time domain. For example, consider a simple linear problem:

where L0 is a constant load vector, and KT, DT, and ET are constant time-domain matrices; see Equation 20-11. In the frequency domain, this corresponds to (with a small change of notation):

The definition in Equation 20-39 now gives:

As mentioned above, the expansion is internally done in λ = (−i) ω, which means that the exported matrices are based on the following:

Notice that for ω0 = 0 the matrices coincide with the time domain but not otherwise.

The Differences Between the Modal Solver and the Modal Reduction Solver

Note that although the Modal solver and the Modal Reduction solver use the same algorithm to get a reduced system, they are different in several ways:

•

The Modal solver solves the reduced system directly and returns the reconstructed solution. The Modal Reduction solver does not solve the reduced system; instead, it creates a reduced-order model, which represents the reduced system and can be solved in the same way as any other physics. The reduced-order model is a for time domain, and for frequency domain.

•

As indicated by the Exporting the Reduced System section above, the reduced matrices can be exported by both solvers. For the Modal solver, it is in the section of the settings for the node. The matrices and vectors are of the full size of the reduced system. For the Modal Reduction solver, you need to go to > , create a node, choose a reduced-order model as the solution, and in the section, choose the matrix to output. The matrices and vectors are of eliminated size, which reflects the current setting in the section in the settings for the reduced-order model node. An exported matrix like Kr is the submatrix of the full-sized reduced matrix, corresponding to all unconstrained reduced DOFs. You can also output the constrained part of the stiffness matrix Krc, which corresponds to the constrained reduced DOFs. If you want to get the full-sized reduced matrices and vectors, then clear all checkboxes in the column in the table.

•

The load factor is set by the Modal solver in its solver settings. There is no load factor in the Modal Reduction solver; it is set in the reduced-order model.

•

The Model Reduction study supports defining output variables, while studies using the Modal solver do not support it.

•

Initial conditions for transient case is different in case the second term in the right side of Equation 20-25 is nonzero.

•

The Modal Reduction solver supports inputs in constraints while the Modal solver does not.

If there are inputs in constraints when you run the Modal Reduction solver, the corresponding constraint modes become a constrained DOF in the reduced equation. You can either solve the equation together with the constraints or solve an eliminated reduced equation in the case that the reduced-order model is stateful (that is, the reduced-order model exposes a set of reduced-order equations and state DOFs to the solver). The eliminated system removes the constraints from the equation.

The eliminated system in the time domain is

(20-46)