Shock and Vibration

Volume 2018, Article ID 1589794, 12 pages

https://doi.org/10.1155/2018/1589794

## Stabilized Solution to Spurious Mode Problem and Ill-Conditioning in Interface Force Based Substructure Coupling Method

School of Mechanical Engineering, Beijing Institute of Technology, Beijing, 100081, China

Correspondence should be addressed to Hui-hua Feng; nc.ude.tib@hhgnef

Received 31 March 2018; Revised 24 May 2018; Accepted 20 June 2018; Published 17 July 2018

Academic Editor: Roberto Nascimbene

Copyright © 2018 Hui-hua Feng et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

There are two major types of substructure mode synthesis methods, i.e., the fixed-interface component mode synthesis and free-interface component mode synthesis. There are two coupling methods, the interface degrees of freedom based coupling method and the interface force based coupling method, the former one is referred to as the primary assembly method, and the latter is referred to as the dual assembly method. However, the dual assembly method is theoretically shown to be unstable in this research, such reduced stiffness matrix is indefinite, this fatal weakness can be conquered by further interface reduction, and the interface compatibility is therefore rigorously enforced. Unfortunately, Craig’s method leads to another numerical instability when inverting a submatrix of residual flexibility on the interface degrees of freedom, this problem is neglectable in small dimensional matrix problems, but it is prominent in large models when the number of interface degrees of freedom is large, this ill-conditioning problem may be circumvented by truncated singular value decomposition technique; here, a more efficient strategy is proposed, the substructure reduction is modified, this modification proves to be numerically stable, and it can be even more accurate than the prevailing Craig-Bampton method; the numerical examples validate the suggestion.

#### 1. Introduction

Modal analysis is very popular in the dynamic community; the experts use this technique in the vibration analysis, dynamic optimization, condition monitoring, damage detection, and many other civil engineering applications. Modal analysis is performed in modal domain, which is similar to the time domain or the frequency domain. The mode shape is simply a description of structural movement, it can be visually seen by adding the mode shape vector to the space coordinates, they are ordered from the lowest frequency all the way to the highest frequency, there are as many orders as the total number of degrees of freedom in the structure, and the ordered modes capture the stiffness and mass distribution of structural design. From a mathematical viewpoint, mode shape just stands for a basis in vector space, so the modal matrix is filled with mode shape vectors, the linear combination of these vectors can be used to approximate the vibration response.

With the development of the finite element technique, and the increase of computer power, many problems can be resolved; however, the need for greater speed and higher efficiency is still a challenge [1] to finish such task, when millions of degrees of freedom are involved, which motivates new computational techniques. The component mode synthesis (CMS) methods [2, 3] are one of the most efficient techniques applicable to modern large dynamic systems [4–6]. First of all, these domain decomposition methods can reduce the overall computational burden and improve the computer efficiency in parallel computations. Secondly, these methods change the way people communicate the design information, engineers can design every component or substructure of a complex system independently, and they can communicate the component mass and stiffness matrix instead of the detailed design information when they are interested in the assembly properties. Last but not least, these techniques are applicable to the fast designs and predictions when the local mesh increases or decreases, the dynamic system can be quickly synthesized and reanalyzed, which makes the CMS methods overwhelm [7] the most efficient structural dynamic modification algorithms.

The Craig-Bampton method [2] is one of the well-known mode synthesis methods, it is accurate and robust, it becomes one of the industrial standards in solving some dynamic problems, the reduced system matrixes are sparse and this property allows it to be implemented in multilayer analysis. However, this method can hardly be validated in experiments, since the fixed-interface condition is difficult to realize. The Craig-Chang method [8, 9] is another CMS method, it deals with free normal modes and provides the potential of experimental validations, and the residual mode effect ensures the first-order accuracy. This method was further improved by introducing the second-order residual mode effect by MacNeal [10]. Later, Rubin [11] advanced MacNeal’s method and improved the component reduction using a consistent reduction base [12], and second-order accuracy may be obtained by solving the frequency dependent nonlinear equations; however, their methods are too complex. Craig and Chang [13] also reformulated this nonlinear problem in a report. N. Bouhaddi [14] reformulated MacNeal’s and Rubin’s presentations in discrete or matrix form, but the nonlinear eigenvalue solution finds it hard to adapt to the forced response analysis. Furthermore, N. Bouhaddi proposed a method to select the master DOFs by a least square solution; the proposed method is excessively complex and adds much more computer time expense. van der Valk [15] simplified Rubin’s method, and the interface force is extracted to correct the lower order normal modes by using the residual modes; however, this simplification encounters the same instability problem of ill-conditioning in large model applications.

One of the most interesting mode synthesis methods is Rixen’s method [15–18], and the dual assembly uses the Lagrange multipliers or interface force as generalized DOFs. Similar matrix sparsity property is preserved, while the interface DOFs are replaced with Lagrange multipliers. Different from MacNeal’s and Rubin’s methods, the range space spanned by the lower order free modes is not corrected by residual modes, and such obtained generalized DOFs are different significantly. Unfortunately, the spurious modes or negative eigenvalues are reported. As a remedy, Rixen proposed [18] to append the first-order nonrigid body Wilson-Ritz vector to the reduction base from inverse iteration. van der Valk [15] discussed the interface forces reduction by extracting* interface force modes*, he called this kind of interface reduction the “*dynamic interface reduction,*” in his study he thought that the displacement compatibility is weakened, and he concluded that “*the Dual-Craig-Bampton and Dual Mixed Craig-Bampton methods can result in spurious modes due to the fact that only a weak interface compatibility is enforced.*” Until now, the mechanism for this spurious mode issue in this method is still undiscovered.

In this investigation, we will discuss this numerical instability problem and show how to avoid the lower order spurious mode. The discussion is organized as follows: In Section 1, we give a quick and brief summary of Rixen’s dual assembly; in Section 2, we discuss the possible reason and solution to the spurious mode issue and how to circumvent the ill-conditioning problem; numerical examples in Section 3 demonstrate the improved accuracy and stability, and finally we end the discussions in Section 4.

The major contribution of this research are as follows: the theoretical proof is provided for the negative eigenvalue issue, and the strong interface compatibility can be satisfied by the projection method or Craig’s method, which shows that a further interface condensation is necessary; the numerical instability of ill-conditioning can be addressed by modifying the substructure reduction, and the residual modes are replaced with the constraint modes, this modification leads to inverting a diagonal matrix with perfect numerical condition, and some numerical examples show that the accuracy is comparable to the Craig-Bampton method.

#### 2. Theoretical Background and Numerical Improvement

##### 2.1. The Dual Assembly and Mode Synthesis

D. J. Rixen showed [16] the mode synthesis method using Lagrange multipliers or interface force degrees of freedom (DOFs). For given s-th substructure, the dynamic property is represented by the reductionwhere is the vector of physical degrees of freedom of the s-th substructure, the matrix is the lower order normal modes from generalized eigenvalue computations, it includes the rigid body modes if the stiffness matrix is rank deficit, is the vector of mode coordinates, the matrix is the residual flexibility matrix composed of residual modes or residual attachment modes, is the Lagrange multiplier vector, and is the interface force. The signed Boolean matrix is applied to the boundary DOFs, respectively, the compound matrix forms the substructure reduction base, and the superscript denotes the matrix transpose.

The displacement compatibility equationwhere is a signed Boolean matrix when conforming mesh is used, (s=1, 2...n) is a submatrix, and is a vector for all physical degrees of freedom of the assembly split into n components.

Introduce the constraints to the dynamic system, to obtain the modified Lagrange functionand apply it to the following Lagrange equation:The system dynamic equations are obtained, and the matrix form iswhere is to build a diagonal block matrix, is the diagonal block mass matrix, is the diagonal block stiffness matrix, and are the mass and stiffness matrix, respectively, and is the summation of all external force vector. The number of equations in (5) is less than the number of unknown variables.

Combine (2) and (5) to obtain the augmented system equations

The negative symbol can be absorbed into the Lagrange multipliers; please refer to original equation (7) in [16]. Different from the classical Craig-Bampton method, the Lagrange multipliers or interface forces are used; hence, Rixen calls it the dual assembly.

The dual assembly uses the following matrix reduction:where the matrix is given to beand is the identity matrix, the dimension is of the same order as the number of Lagrange multipliers.

Substitute (7) and (8) to (6), and we obtain the following reduced equation:where

The condensed matrixes preserve some sparsity, there is no coupling between the mode coordinates and the Lagrange multipliers in the condensed mass matrix, and some coupling terms occur in the condensed stiffness matrix. Unfortunately, the Dual-Craig-Bampton method brings spurious modes as negative eigenvalues, which leads to imaginary frequencies, and this misleading phenomenon implies some incorrectness.

The Lagrange multipliers are introduced as interface force, the boundary degrees of freedom doubled in the process of partitioning the domain into substructures, and the augmented system equations contain redundant degrees of freedom. There are many strategies to cope with these redundant degrees of freedom in different applications. Before we go further, we review some existing methods to solve these algebraic-differential equations (ADE).

The master-slave degree cancelation method is accurate and robust; only a subset of the generalized DOFs is condensed and solved. This method is prevailing in commercial FEM analysis with application to linear multipoint constraints [19, 20]; the latter are introduced by some rigid elements (e.g., RBE2), but this method partially destroys the matrix sparsity and profile. The penalty method [21, 22] can also be used without changing the matrix profile and sparsity when dealing with constraints. The two methods, however, are not modal based.

The trapezoidal rule using the Newmark method in the time domain direct integration is preferred in structural dynamic analysis, but it is unconditionally unstable [23] in ADE systems, to recover stability some numerical damping is necessary, the alpha method or the Hibler-Huges-Taylor method with can be used, the Lagrange multipliers are kept in this direct integration method. Besides, [24] discusses how to discretize the equations and eliminate the redundant Lagrange multipliers to ensure stability, so that general Newmark method can be used with quadratic accuracy.

Ritz solution can be obtained as an approximation to such ADE system, this is the central topic in this research, the generalized method is proposed by Craig, see [3, 9], this method tries to pick up the independent generalized coordinates and condense the dependent generalized coordinates, and we will show that it is equivalent to using a null space matrix or projection matrix to eliminate the Lagrange multipliers.

There are other direct methods [25–27] that can provide alternative solutions, Lidström [25] discusses the system with “n” DOFs and “c” linear homogeneous constraints, and when the Lagrange multipliers are eliminated, there are “n-c” modes left to describe the dynamics of the constrained system. However, this method results in a nonsymmetric singular stiffness matrix; this weakness limits its application. There is another method [26] that provides the possibility of eliminating the Lagrange multipliers directly; however, the singular impedance matrix problem and twice matrix inverse operations make this method less attractive in frequency domain applications. Mark [27] discusses a similar problem in nonconforming mesh modeling, and his method can be extended to the dynamic problems by the matrix transforming , where is a singular matrix. The reader should consult the original paper [27] for detailed definition of the matrix . However, the weakness of this directly extended version is that this usually generates a nonsymmetric singular mass matrix. Both Lidström’s method and Mark’s generalized version will encounter singular matrixes and therein the troublesome numerical modes, or the “undefined modes” according to Lidström.

##### 2.2. The Indefinite Stiffness Matrix Leads to Negative Eigenvalue Problem

The interface displacement compatibility principle and the interface force balance principle are two major rules in the development of CMS methods; the two principles can be verified posteriorly by comparing the solution with accurate results. The comparison reveals that the lower frequency results converge faster than the higher frequency results, which means that the two principles are frequency dependent.

Rixen thought [16] that the interface displacement compatibility is weakened during the assembly process. However, we think the interface force balance is weakened, and this weakening mechanism comes from the approximations applied to substructure reduction, unreasonable Ritz vectors or not enough Ritz vectors lead to a lower interface force resolution, and, as a result, the interface force balance is enforced mathematically but not physically; for example, if we denote the accurate interface balance as at a specific frequency for given DOF, the subscript “ab” refers to the interface force applied to substructure a by substructure b and vice versa. Due to approximations, we have weakened interface balance equation, , where and are the weakened interface forces, and the deviation is . Note that both the accurate interface forces and the weakened interface forces satisfy the interface force balance equation accurately.

The theoretical explanation of “spurious mode” issue is that the reduced system stiffness matrix is indefinite. The indefinite stiffness matrix leads to the negative eigenvalues; hence, the latter can be alleviated neither by more normal modes nor by the load dependent vectors or Wilson vectors.

If the normal modes are mass normalized, from (9), we can rewrite the reduced system aswhere the subscript i refers to the mode coordinates and the subscript j refers to the Lagrange multipliers. The number of Lagrange multipliers equals the total number of independent interface DOFs. The coupling terms are

We prove that the reduced stiffness matrix is indefinite provided that .

Lemma 1. *If and , then, for an arbitrary nonzero matrix , .*

*Proof. *Consider the quadratic form, , for and , by the SVD to obtain , where is the matrix of singular vectors, and the diagonal matrix , and we have , with the definition by .

Further assume that , and we have ; hence, we can assert that .

Lemma 2. *Assume matrix and , and then, for an arbitrary matrix , the augmented matrix is indefinite.*

*Proof. *Again, for and subject to , consider the quadratic formand we just consider the special casesHence, substitute the partitioned block matrixes, and we can find that the reduced stiffness matrix is indefinite. The reduced mass matrix is definite; the proof is similar and omitted here.

Lemma 3. *The system composed of indefinite stiffness matrix and definite mass matrix has both negative eigenvalue and positive eigenvalue.*

*Proof. *Using the maximum and minimum property of Rayleigh quotient (Rq), i.e., we have , where is the smallest eigenvalue and is the largest eigenvalue, together with the property of quadratic form of the two matrixes, we can easily prove that the minimum eigenvalue is negative, and the maximum eigenvalue is positive.

Now, we end the discussion in this section with the conclusion that (11) cannot be solved directly by modal solution strategy, otherwise it leads to negative eigenvalue, and this negative eigenvalue issue is not caused by nonlinear effects; therefore, it is unacceptable. In the following discussion, we will show how to obtain a solution without negative eigenvalue.

##### 2.3. Enforce the Interface Displacement Compatibility by Craig’s Method

Now we start with (5), instead of (6), to show that it is possible to construct a null space matrix directly; the Lagrange multipliers are eliminated directly to obtain the symmetric system matrixes. This theoretical analysis can give more insights into the mode synthesis process, and this derivation can be further generalized to geometrical nonlinearity.

Consider the following mode synthesis transformation [9, 28]where is the independent generalized coordinate, and the mode synthesis matrix isand other matrixes are defined by

The matrix can be a properly selected matrix so that the matrix is invertible; in Craig’s method [9], this matrix is chosen as the residual flexibility matrix or the inertia relief attachment mode; the constraint matrix is assumed to be in full rank, i.e., the rows are linearly independent, and usually this property can be guaranteed. We have the following remarks.

*Remark 4. *If we define and , then are projection matrixes, and this can be verified to beandIf the rank of is , i.e., there are constraints to describe the displacement compatibility, assume the rank of is larger than* c*, we can have , and is the total number of DOFs of s-th substructure.

*Remark 5. *Matrixes lie in the null space of matrix , i.e., , and we haveand similarlySubstitute (15) to (5), left multiply by , and we can have the reduced matrix equationwhere we have used the property . The semidefinite property of the reduced stiffness ensures that the reduced system has no negative eigenvalue, and this is a Ritz process, the projection matrix enforces the solution to satisfy the given linear constraints, i.e., , and such obtained matrixes are symmetric as well.with

,

*Remark 6. *The null space matrix can be decoupled into two matrixeswhereObviously from this decomposition, using the null space matrix is equivalent to applying the matrixes and to (5) sequentially in reduction, the matrix is closely related to substructure modal reduction, and further reduction with the matrix can be considered as a general modal synthesis process; in other words, this modal synthesis process is achieved by a further interface reduction, and this further interface reduction enforces the strong interface displacement compatibility.

*Remark 7. *The residual modes are selected as a special case. When the matrix is selected to be the residual flexibility matrix, i.e., , the matrix becomes the Guyan reduction matrix; due to the orthogonal property, the condensed system matrixes can be simplifiedUsually, the sparsity is lost in the reduced matrixes.

*Remark 8. *The above derivation can be generalized to consider the nonhomogeneous multipoint constraints, and the latter may be introduced by some geometrical nonlinearityIntroduce a new state variableand the time derivative givesWe obtain the new transformed system equations, and the governing equations areNow the power of the new derivations becomes noticeable, and the above transformed equation can be resolved in the same fashion. After is computed, can be retrieved asand the interface forces are

*Remark 9. *The condition for a unique solution is not discussed in Craig’s presentation, and now we simply append this condition without proofThis condition implies that the eigenvectors of original unconstrained system do not lie in the null space of constraint matrix, and it applies both to (5) and to the nonhomogeneous case in (31). If this condition is violated, this implies some erroneous modeling.

##### 2.4. Numerical Instability in Craig’s Method and Final Suggestion by Modifying the Substructure Reduction

Though the residual modes can simplify the reduced matrixes, extra computations prior to the condensation for every substructure are required, this may be the limitation for large scale application where there are numerous boundary DOFs, and it is computationally expensive to obtain the matrix , especially for the singular stiffness matrix. Another fatal weakness of using the residual modes is that the matrix is ill-conditioned sometimes, invert this matrix is unstable, this problem is inherent in both Craig-Chang method [8] and the simplified MacNeal-Rubin method [15], and this ill-conditioning is subtle in small models, but it becomes prominent in large models according to our experience.

The projection matrix depends on the constraint matrix and the substructure reduction. A series of results can be obtained if we modify the substructure reduction; for example, if we choose , the mode synthesis matrix can also be quickly generated. Different choices of result in different convergence character. The identity matrix may be favorable for no extra computation a priori; i.e., if we set as identity matrix, the mode synthesis matrix becomes

At the same time, the numerical instability of inverting the matrix is alleviated completely.

There is a good candidate to balance the numerical stability and convergence, the columns of this matrix on the boundary DOFs are the constraint modes, and the substructure reduction is modified to becomewhere SY designates to keep symmetry, and only the columns related to the boundary degrees of freedom are utilized. For clarity, the proposed matrix is rewritten as below

If the conforming mesh is utilized, then the matrix becomes the constraint modes, and such obtained matrix is a diagonal matrix, and it is perfectly conditioned with positive integers on the diagonal. By the way, the ingredients of (11) are valid only for the residual flexibility matrix , the submatrixes vary with different . While (15)~(24) keep in universal form, they are valid for residual flexibility matrix and the proposed in formula (37).

#### 3. Numerical Verification Examples

##### 3.1. A 10-DOF Model

A 10-DOF lumped mass model is shown in Figure 1, the square box is the lumped mass, the wave line is the spring, the straight line indicates the rigid link, and only the horizontal dynamics are considered. The DOFs are numbered from the 1st to 10th. Partition the assembly into 4 components or substructures, as shown in Figure 2, the 4th DOF is split into 4 DOFs, and the mass is equally distributed to the 4 substructures. The mass for every DOF and spring rate for every connection are labeled in the nearby blank. Note this cutting operation is only mathematically, and they are still physically connected to constitute a whole assembly; therefore, the dynamic properties of the partitioned system keep the same to the original system. The system now has =13 DOFs together with additional =3 rigid connections.