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

Authors: Hui-hua Feng; Xiao-long Mou; Khaled Teffah
Journal: Shock and Vibration
 (2018)

---

## 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.


---

## Body

## 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 extractinginterface 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 Section1, 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:(1) 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; (2) 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 reduction(1)xs≈Φsq+RsbsTλ=ΦsRsbsTqλwhere xs is the vector of physical degrees of freedom of the s-th substructure, the matrix Φs is the lower order normal modes from generalized eigenvalue computations, it includes the rigid body modes if the stiffness matrix is rank deficit, q is the vector of mode coordinates, the matrix Rs is the residual flexibility matrix composed of residual modes or residual attachment modes, λ is the Lagrange multiplier vector, and bsTλ is the interface force. The signed Boolean matrix bs is applied to the boundary DOFs, respectively, the compound matrix ΦsRsbsT forms the substructure reduction base, and the superscript T denotes the matrix transpose.The displacement compatibility equation(2)BX=0where B=b1b2⋯bn is a signed Boolean matrix when conforming mesh is used, bs (s=1, 2...n) is a submatrix, and X=x1Tx2T⋯xnTT 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 function(3)l=12X˙TMX˙-12XTKX+λTBXand apply it to the following Lagrange equation:(4)ddt∂l∂q˙-∂l∂q=fThe system dynamic equations are obtained, and the matrix form is(5)MX¨+KX-BTλ=Fwhere diag⁡· is to build a diagonal block matrix, M=diag⁡Ms is the diagonal block mass matrix, K=diag⁡Ks is the diagonal block stiffness matrix, Ms and Ks are the mass and stiffness matrix, respectively, and F=∑fs 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(6)M000X¨λ¨+K-BT-B0Xλ=F0The 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:(7)Xλ=x1⋮xnλ=ΨDCBq1⋮qnλwhere the matrix ΨDCB is given to be(8)ΨDCB=Φ1R1b1T⋱⋮ΦnRnbnTIλ×λand I 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:(9)M-q¨-+K-q-=F-where(10)M-=ΨDCBTM000ΨDCBK-=ΨDCBTK-BT-B0ΨDCBF-=ΨDCBTF0q-=q1Tq2T⋯qnTλTTThe 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).( 1 ) 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.( 2 ) 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 α≠0 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.( 3 ) 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.( 4 ) 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 M^=QTM, where Q is a singular matrix. The reader should consult the original paper [27] for detailed definition of the matrix Q. 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 Fabω+Fbaω=0 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, F~abω+F~baω=0, where F~abω=Fabω+Δ and F~baω=Fbaω-Δ 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 as(11)I00M-jjq¨λ¨+ΛK-ijK-jiK-jjqλ=F-iF-jwhere 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(12)M-jj=∑bsRsMsRsbsTK-ij=-b1Φ1-b2Φ2⋯-bnΦnT,K-ji=K-ijT,K-jj=∑bsRsKsRsbsT-2∑bsRsbsT=-∑bsRsbsT=-Bdiag⁡RsBTF-i=f1TΦ1f2TΦ2⋯fnTΦnTF-j=∑bsRsfsWe prove that the reduced stiffness matrix is indefinite provided thatRsbsT≠0.Lemma 1.
IfRs≥0 and RsbsT≠0, then, for an arbitrary nonzero matrix bs, bsRsbsT>0.Proof.
Consider the quadratic form,∏y=yTbsRsbsTy, for ∀y∈Rm and y≠0, by the SVD to obtain Rs=UsSsUsT, where Us is the matrix of singular vectors, and the diagonal matrix Ss≥0, and we have ∏y=yTbsUsSsUsTbsTy=vTSsv=∑iSs,iivi2≥0, with the definition by v≜UsTbsTy.
Further assume thatRsbsT≠0, and we have ∏y>0; hence, we can assert that K-jj=-∑bsRsbsT<0.Lemma 2.
Assume matrixΛm×m≥0 and Cn×n<0, and then, for an arbitrary matrix H, the augmented matrix ΛHTHC is indefinite.Proof.
Again, for∀y1∈Rm,y2∈Rn and subject to y1Ty1+y2Ty2≠0, consider the quadratic form(13)∏y1,y2=y1y2TΛHTHCy1y2
and we just consider the special cases(14)∏0,y2≠0=y1y2TΛHTHCy1y2=y2TCy2<0∏y1≠0,0=y1y2TΛHTHCy1y2=y1TΛy1≥0
Hence, 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λ1≤Rq≤λn, where λ1 is the smallest eigenvalue and λn is the largest eigenvalue, together with the property of quadratic form of the two matrixes, we can easily prove that the minimum eigenvalue λ1 is negative, and the maximum eigenvalue λn 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](15)X=SqIwhere qI is the independent generalized coordinate, and the mode synthesis matrix is(16)S=I-diag⁡GsBTA-1Bdiag⁡Φsand other matrixes are defined by(17)I=diag⁡ΦsTMsΦs=diag⁡IsA=Bdiag⁡GsBTThe matrixGs can be a properly selected matrix so that the matrix A 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 B 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 defineQ=I-P and P=diag⁡GsBTA-1B, then P&Q are projection matrixes, and this can be verified to be(18)P2=diag⁡GsBTA-1Bdiag⁡GsBTA-1B=diag⁡GsBTA-1AA-1B=P
and(19)Q2=I-P2=I-2P+P2=I-P=Q
If the rank ofB is rB=c, i.e., there are c constraints to describe the displacement compatibility, assume the rank of diag⁡Gs is larger thanc, we can have rQ=∑ns-c, and ns is the total number of DOFs of s-th substructure.Remark 5.
MatrixesQ&S lie in the null space of matrix B, i.e., Q,S∈NB=∀u∈R,Bu=0, and we have(20)BQ=BI-diag⁡GsBTA-1B=B-AA-1B=0
and similarly(21)BS=BQdiag⁡Φs=0
Substitute (15) to (5), left multiply by ST, and we can have the reduced matrix equation(22)STMSq¨I+STKSqI=STFwhere we have used the property BS=0. 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., BX=0, and such obtained matrixes are symmetric as well.(23)STMS=I-diag⁡ΦsTBTA-TBdiag⁡GsTMsΦs-diag⁡ΦsTMsGsBTA-1Bdiag⁡Φs+diag⁡ΦsTBTA-TBdiag⁡GsTMsGsBTA-1Bdiag⁡ΦsSTKS=Λ-diag⁡ΦsTBTA-TBdiag⁡GsTKsdiag⁡Φs-diag⁡ΦsTdiag⁡KsGsBTA-1Bdiag⁡Φs+diag⁡ΦsTBTA-TBdiag⁡GsTKsGsBTA-1Bdiag⁡ΦswithΛ = diag ⁡ Φ s T K s Φ s, Bdiag⁡GsTMsGsBT=∑bsGsTMsGsbsTRemark 6.
The null space matrix can be decoupled into two matrixes(24)S=I-diag⁡GsBTA-1Bdiag⁡Φs=diag⁡Φsdiag⁡GsI-BTA-1Bdiag⁡Φs=Φ1G1b1T⋱⋮ΦNGNbNTI-A-1Bdiag⁡Φs=Ψcc∗Tg∗where(25)Ψcc∗=Φ1G1b1T⋱⋮ΦNGNbNT,Tg∗=I-A-1Bdiag⁡Φs
Obviously from this decomposition, using the null space matrix is equivalent to applying the matrixesΨcc∗ and Tg∗ to (5) sequentially in reduction, the matrix Ψcc∗ is closely related to substructure modal reduction, and further reduction with the matrix Tg∗ 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 matrixGs is selected to be the residual flexibility matrix, i.e., Gs=Rs, the matrix Tg∗ becomes the Guyan reduction matrix; due to the orthogonal property, the condensed system matrixes can be simplified(26)STMS=I+diag⁡ΦsTBTA-1Bdiag⁡RsMsRsBTA-1Bdiag⁡ΦsSTKS=Λ+diag⁡ΦsTBTA-1Bdiag⁡Φs
Usually, 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 nonlinearity(27)BX=g
Introduce a new state variable(28)y=x-Pg
and the time derivative gives(29)y¨=x¨
We obtain the new transformed system equations, and the governing equations are(30)By=0(31)My¨+Ky-BTλ=F-KPg
Now the power of the new derivations becomes noticeable, and the above transformed equation can be resolved in the same fashion. Aftery is computed, x can be retrieved as(32)x=Sq+Pg
and the interface forces are(33)λ=BT-1F-Mx¨-Kx+v∀v∈NBTRemark 9.
The condition for a unique solution is not discussed in Craig’s presentation, and now we simply append this condition without proof(34)φ:∀ω,K-ω2Mφ=0∩NB=0
This 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 matrixGsbsT, especially for the singular stiffness matrix. Another fatal weakness of using the residual modes is that the matrix A=Bdiag⁡GsBT 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 matrixB and the substructure reduction. A series of results can be obtained if we modify the substructure reduction; for example, if we choose Gs=Ms, the mode synthesis matrix can also be quickly generated. Different choices of Gs result in different convergence character. The identity matrix may be favorable for no extra computation a priori; i.e., if we set Gs as identity matrix, the mode synthesis matrix becomes(35)S=I-BTBBT-1Bdiag⁡ΦsAt the same time, the numerical instability of inverting the matrixA=BBT 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 become(36)xs=Φsq+0ii-Kii-1KibSYIbbsbsTλwhere 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(37)Gs=0ii-Kii-1KibSYIbbsIf the conforming mesh is utilized, then the matrixGsbsT becomes the constraint modes, and such obtained matrix A 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 Rs, the submatrixes vary with different Gs. While (15)~(24) keep in universal form, they are valid for residual flexibility matrix and the proposed Gs in formula (37).


## 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 reduction(1)xs≈Φsq+RsbsTλ=ΦsRsbsTqλwhere xs is the vector of physical degrees of freedom of the s-th substructure, the matrix Φs is the lower order normal modes from generalized eigenvalue computations, it includes the rigid body modes if the stiffness matrix is rank deficit, q is the vector of mode coordinates, the matrix Rs is the residual flexibility matrix composed of residual modes or residual attachment modes, λ is the Lagrange multiplier vector, and bsTλ is the interface force. The signed Boolean matrix bs is applied to the boundary DOFs, respectively, the compound matrix ΦsRsbsT forms the substructure reduction base, and the superscript T denotes the matrix transpose.The displacement compatibility equation(2)BX=0where B=b1b2⋯bn is a signed Boolean matrix when conforming mesh is used, bs (s=1, 2...n) is a submatrix, and X=x1Tx2T⋯xnTT 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 function(3)l=12X˙TMX˙-12XTKX+λTBXand apply it to the following Lagrange equation:(4)ddt∂l∂q˙-∂l∂q=fThe system dynamic equations are obtained, and the matrix form is(5)MX¨+KX-BTλ=Fwhere diag⁡· is to build a diagonal block matrix, M=diag⁡Ms is the diagonal block mass matrix, K=diag⁡Ks is the diagonal block stiffness matrix, Ms and Ks are the mass and stiffness matrix, respectively, and F=∑fs 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(6)M000X¨λ¨+K-BT-B0Xλ=F0The 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:(7)Xλ=x1⋮xnλ=ΨDCBq1⋮qnλwhere the matrix ΨDCB is given to be(8)ΨDCB=Φ1R1b1T⋱⋮ΦnRnbnTIλ×λand I 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:(9)M-q¨-+K-q-=F-where(10)M-=ΨDCBTM000ΨDCBK-=ΨDCBTK-BT-B0ΨDCBF-=ΨDCBTF0q-=q1Tq2T⋯qnTλTTThe 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).( 1 ) 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.( 2 ) 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 α≠0 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.( 3 ) 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.( 4 ) 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 M^=QTM, where Q is a singular matrix. The reader should consult the original paper [27] for detailed definition of the matrix Q. 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 Fabω+Fbaω=0 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, F~abω+F~baω=0, where F~abω=Fabω+Δ and F~baω=Fbaω-Δ 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 as(11)I00M-jjq¨λ¨+ΛK-ijK-jiK-jjqλ=F-iF-jwhere 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(12)M-jj=∑bsRsMsRsbsTK-ij=-b1Φ1-b2Φ2⋯-bnΦnT,K-ji=K-ijT,K-jj=∑bsRsKsRsbsT-2∑bsRsbsT=-∑bsRsbsT=-Bdiag⁡RsBTF-i=f1TΦ1f2TΦ2⋯fnTΦnTF-j=∑bsRsfsWe prove that the reduced stiffness matrix is indefinite provided thatRsbsT≠0.Lemma 1.
IfRs≥0 and RsbsT≠0, then, for an arbitrary nonzero matrix bs, bsRsbsT>0.Proof.
Consider the quadratic form,∏y=yTbsRsbsTy, for ∀y∈Rm and y≠0, by the SVD to obtain Rs=UsSsUsT, where Us is the matrix of singular vectors, and the diagonal matrix Ss≥0, and we have ∏y=yTbsUsSsUsTbsTy=vTSsv=∑iSs,iivi2≥0, with the definition by v≜UsTbsTy.
Further assume thatRsbsT≠0, and we have ∏y>0; hence, we can assert that K-jj=-∑bsRsbsT<0.Lemma 2.
Assume matrixΛm×m≥0 and Cn×n<0, and then, for an arbitrary matrix H, the augmented matrix ΛHTHC is indefinite.Proof.
Again, for∀y1∈Rm,y2∈Rn and subject to y1Ty1+y2Ty2≠0, consider the quadratic form(13)∏y1,y2=y1y2TΛHTHCy1y2
and we just consider the special cases(14)∏0,y2≠0=y1y2TΛHTHCy1y2=y2TCy2<0∏y1≠0,0=y1y2TΛHTHCy1y2=y1TΛy1≥0
Hence, 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λ1≤Rq≤λn, where λ1 is the smallest eigenvalue and λn is the largest eigenvalue, together with the property of quadratic form of the two matrixes, we can easily prove that the minimum eigenvalue λ1 is negative, and the maximum eigenvalue λn 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](15)X=SqIwhere qI is the independent generalized coordinate, and the mode synthesis matrix is(16)S=I-diag⁡GsBTA-1Bdiag⁡Φsand other matrixes are defined by(17)I=diag⁡ΦsTMsΦs=diag⁡IsA=Bdiag⁡GsBTThe matrixGs can be a properly selected matrix so that the matrix A 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 B 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 defineQ=I-P and P=diag⁡GsBTA-1B, then P&Q are projection matrixes, and this can be verified to be(18)P2=diag⁡GsBTA-1Bdiag⁡GsBTA-1B=diag⁡GsBTA-1AA-1B=P
and(19)Q2=I-P2=I-2P+P2=I-P=Q
If the rank ofB is rB=c, i.e., there are c constraints to describe the displacement compatibility, assume the rank of diag⁡Gs is larger thanc, we can have rQ=∑ns-c, and ns is the total number of DOFs of s-th substructure.Remark 5.
MatrixesQ&S lie in the null space of matrix B, i.e., Q,S∈NB=∀u∈R,Bu=0, and we have(20)BQ=BI-diag⁡GsBTA-1B=B-AA-1B=0
and similarly(21)BS=BQdiag⁡Φs=0
Substitute (15) to (5), left multiply by ST, and we can have the reduced matrix equation(22)STMSq¨I+STKSqI=STFwhere we have used the property BS=0. 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., BX=0, and such obtained matrixes are symmetric as well.(23)STMS=I-diag⁡ΦsTBTA-TBdiag⁡GsTMsΦs-diag⁡ΦsTMsGsBTA-1Bdiag⁡Φs+diag⁡ΦsTBTA-TBdiag⁡GsTMsGsBTA-1Bdiag⁡ΦsSTKS=Λ-diag⁡ΦsTBTA-TBdiag⁡GsTKsdiag⁡Φs-diag⁡ΦsTdiag⁡KsGsBTA-1Bdiag⁡Φs+diag⁡ΦsTBTA-TBdiag⁡GsTKsGsBTA-1Bdiag⁡ΦswithΛ = diag ⁡ Φ s T K s Φ s, Bdiag⁡GsTMsGsBT=∑bsGsTMsGsbsTRemark 6.
The null space matrix can be decoupled into two matrixes(24)S=I-diag⁡GsBTA-1Bdiag⁡Φs=diag⁡Φsdiag⁡GsI-BTA-1Bdiag⁡Φs=Φ1G1b1T⋱⋮ΦNGNbNTI-A-1Bdiag⁡Φs=Ψcc∗Tg∗where(25)Ψcc∗=Φ1G1b1T⋱⋮ΦNGNbNT,Tg∗=I-A-1Bdiag⁡Φs
Obviously from this decomposition, using the null space matrix is equivalent to applying the matrixesΨcc∗ and Tg∗ to (5) sequentially in reduction, the matrix Ψcc∗ is closely related to substructure modal reduction, and further reduction with the matrix Tg∗ 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 matrixGs is selected to be the residual flexibility matrix, i.e., Gs=Rs, the matrix Tg∗ becomes the Guyan reduction matrix; due to the orthogonal property, the condensed system matrixes can be simplified(26)STMS=I+diag⁡ΦsTBTA-1Bdiag⁡RsMsRsBTA-1Bdiag⁡ΦsSTKS=Λ+diag⁡ΦsTBTA-1Bdiag⁡Φs
Usually, 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 nonlinearity(27)BX=g
Introduce a new state variable(28)y=x-Pg
and the time derivative gives(29)y¨=x¨
We obtain the new transformed system equations, and the governing equations are(30)By=0(31)My¨+Ky-BTλ=F-KPg
Now the power of the new derivations becomes noticeable, and the above transformed equation can be resolved in the same fashion. Aftery is computed, x can be retrieved as(32)x=Sq+Pg
and the interface forces are(33)λ=BT-1F-Mx¨-Kx+v∀v∈NBTRemark 9.
The condition for a unique solution is not discussed in Craig’s presentation, and now we simply append this condition without proof(34)φ:∀ω,K-ω2Mφ=0∩NB=0
This 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 matrixGsbsT, especially for the singular stiffness matrix. Another fatal weakness of using the residual modes is that the matrix A=Bdiag⁡GsBT 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 matrixB and the substructure reduction. A series of results can be obtained if we modify the substructure reduction; for example, if we choose Gs=Ms, the mode synthesis matrix can also be quickly generated. Different choices of Gs result in different convergence character. The identity matrix may be favorable for no extra computation a priori; i.e., if we set Gs as identity matrix, the mode synthesis matrix becomes(35)S=I-BTBBT-1Bdiag⁡ΦsAt the same time, the numerical instability of inverting the matrixA=BBT 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 become(36)xs=Φsq+0ii-Kii-1KibSYIbbsbsTλwhere 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(37)Gs=0ii-Kii-1KibSYIbbsIf the conforming mesh is utilized, then the matrixGsbsT becomes the constraint modes, and such obtained matrix A 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 Rs, the submatrixes vary with different Gs. While (15)~(24) keep in universal form, they are valid for residual flexibility matrix and the proposed Gs in formula (37).


## 3. Numerical Verification Examples

### 3.1. A 10-DOF Model

A 10-DOF lumped mass model is shown in Figure1, 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 ∑ns=13 DOFs together with additional c =3 rigid connections.Figure 1
Undamped lumped parameter mass model, only the horizontal dynamics are considered, the global DOFs are labeled in the box, the spring rate K=1E4 N/m, and the unit mass M=1 kg.Figure 2
Mathematically partitioned four substructures; the local DOFs are labeled in the box.Before jumping into details, we show how to build the constraint matrix correctly. There are 4 components sharing one single boundary DOF, therefore, the constraint matrix is nonunique. We can build the displacement compatibility equations asx14=x21, x21=x41, and x33=x41, where the superscript is the DOF number, and the subscript is the substructure number. The obtained constraint matrix is(38)
B
=
b
1
0
0
0
1
0
0
0
0
0
0
0
0
⋮
⋮
⋮
⋮
b
2
-
1
0
0
1
0
0
0
0
0
⋮
⋮
⋮
⋮
b
3
0
0
0
0
0
0
0
0
1
⋮
⋮
⋮
⋮
b
4
0
0
0
-
1
0
0
-
1
0
0We can also build the displacement compatibility equations asx14=x21, x21=x41, x33=x41, and x14=x41, and the following constraint matrix is obtained:(39)
B
=
b
1
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
1
⋮
⋮
⋮
⋮
⋮
b
2
-
1
0
0
1
0
0
0
0
0
0
0
0
⋮
⋮
⋮
⋮
⋮
b
3
0
0
0
0
0
0
0
0
1
0
0
0
⋮
⋮
⋮
⋮
⋮
b
4
0
0
0
-
1
0
0
-
1
0
0
-
1
0
0We have built two constraint matrixes; the first one is legal while the second one is illegal. The constraint matrix given by formula (38) is full rank in rows, and it is a complete description of the system constraints, while formula (39) has a redundant constraint, or it is rank deficit in rows. If we use the singular constraint matrix in formula (39), we have to invert a singular matrix A=Bdiag⁡GsBT, hence obtaining inferior results.Table1 listed the mode analysis results. A vector was used to denote the total number of the lower order kept normal modes (including rigid mode), and the vector [2,1,1,1] indicates that we kept 2 orders for the first component and 1 order for every other three components, respectively. TheExact (see Table1) refers to the reference result computed from original system using conventional method. TheDCB refers to the Dual-Craig-Bampton method, from (9), and theDCB+Guyan refers to the results from a further static condensation, i.e., from the (22), and the imaginary unit is i=-1.Table 1
The mode frequency comparison in different methods /Hz.
Order Exact DCB DCB+Guyan [ 2,1 , 1,1 ] [ 3,2 , 2,2 ] [ 2,1 , 1,1 ] [ 3,2 , 2,2 ] -3  26.4556i 62.1508i -2  22.8070i 55.1870i -1  20.6766i 45.5169i 1 0.0000 0.0000 0.0000 0.0000 0.0000 2 6.5712 6.5984 6.5715 6.5740 6.5712 3 7.0001 7.0626 7.0007 7.0331 7.0001 4 7.6090 7.7488 7.6098 7.6828 7.6090 5 12.8444 - 12.8509 - 12.8445 6 18.7891 19.6764 18.8169 18.1328 18.7894 7 27.0282  28.3200  27.0636 8 32.0634  47.2122  32.0652 9 33.1482  -  33.6621 10 35.4645  55.3266We can infer from Table1 that spurious modes appear in the DCB method; we marked the spurious modes in bold. Some spurious modes result from the negative eigenvalue, and they are in imaginary frequencies. The 47 Hz and 55 Hz modes are far beyond the highest frequency of the system, and obviously they are spurious modes as well. There is no spurious mode after the proposed static condensation, and the mode frequency accuracy even increases significantly, which shows better convergence property.It is interesting to consider the effects of differentGs. We tried the identity matrix Is, mass matrix Ms, and stiffness matrix Ks in the numerical examples, we listed the corresponding mode frequency results in Table 2, and we can find that there is no spurious mode at least. Essentially, selecting different Gs matrix affects the substructure representation, the interface force balance is weakened physically even though the displacement compatibility is enforced mathematically, and, consequently, these results are inferior to the results from the residual modes.Table 2
The mode frequency comparison using different matrixGs; [3,2,2,2] orders were kept.
Order Exact Gs=Is Gs=Ms Gs=Ks Gs=Rs 1 0.0000 0.0000 0.0000 0.0000 0.0000 2 6.5712 6.5713 6.5713 6.5714 6.5712 3 7.0001 7.0022 7.0022 7.0226 7.0001 4 7.6090 7.6131 7.6131 7.6103 7.6090 5 12.8444 12.8487 12.8487 12.8537 12.8445 6 18.7891 18.7935 18.7935 18.8111 18.7894 7 27.0282 27.1237 27.1237 27.2344 27.0636 8 32.0634 32.0658 32.0658 32.1554 32.0652 9 33.1482 33.7795 33.7795 33.5493 33.6621 10 35.4645 Averaged Relative error (%)  0.30095 0.30095 0.348863 0.211159For further validation, the mode shape accuracy was analyzed. We checked the orthogonality by the MAC matrix analysis [28], this criterion is preferred in modal analysis, and we used [3,2,2,2] kept orders for components in subsequent discussions. The mode shapes were expanded by a series of matrix multiplications. The MAC matrix between exact result and theDCB result was shown in Figure 3, obviously there are three spurious modes in the lowest orders, followed by 7 order modes, and the MAC matrix is diagonally dominated (except for the lowest 3 orders and the highest 2 orders). There are 12 order modes to describe the 10-DOF-system dynamics since the system is augmented by 3 Lagrange multipliers, we have (3+2+2+2) +3=12, and, of course, it is irrational to describe the 10-DOF system with 12 order modes in linear analysis.Figure 3
MAC analysis between accurate shapes and shapes from DCB method.We present the MAC matrix between exact result and theDCB+Guyan result in Figure 4, we can find that no spurious mode exists in the lowest frequency range, 8 orders have good accuracy, and the last order is acceptable; this comparison shows that the static condensation successfully improved the mode shape accuracy.Figure 4
MAC analysis between accurate shapes and shapes from DCB+Guyan.To our knowledge, the MAC analysis is simply a crude understanding of mode shape accuracy, it is insufficient to judge the overall accuracy, therefore, a more rigorous check was followed by the frequency response function (FRF) analysis, the acceleration FRFs were computed, and these FRFs are weighted byω2 in relative to displacement FRFs; therefore they are more sensitive in detecting the shape errors in higher frequency range. For illustrating purpose, we assigned the assumed modal damping 0.01 to every order, and we removed the three lowest order spurious modes in the DCB method; in addition, we removed the identified rigid body mode. Finally, we selected three FRFs, as shown in Figures 5–7.Figure 5
Acceleration FRF comparison, H(2,9).Figure 6
Acceleration FRF comparison, H(5,2).Figure 7
Acceleration FRF comparison, H(6,4).From Figure5, we can find that the exact modes dominate the frequency range within 6~36Hz, the frequency response above the 36Hz is contaminated by the two spurious modes in the DCB method, and they are 47Hz and 55Hz, respectively. The results for the lower frequency analysis were shown in Figures 6 and 7, and 10~20dB improvement was achieved in some frequencies. The above results validated the improvements of the proposed static condensation.


### 3.2. A Bolted Beam Model

The beam model and its mesh are shown in Figure8, respectively; the assembly is composed of 3 components, upper beam, down beam, and the bolt. One end is fixed ideally, and the other is free. The bolt connects the upper part and the lower part, and all the contact interfaces are treated as linear constraint using conforming mesh; i.e., they are glued together. Now we split the assembly into 4 substructures to perform modal synthesis analysis and generate the reference result from the full model analysis.Figure 8
The suspended beam and its mesh model.This beam model is used here to show the convergence property using differentGs and in this example the identity matrix, substructure mass matrix, substructure stiffness matrix, and the proposed matrix were compared. The 4 matrixes share the similar numerical stability and only differ in the convergence, and additional computation effort is required for the proposed matrix. The residual flexibility matrix and the inertia relief attachment modes are not applicable in this example due to ill-conditioning.Figure9 shows that when 100 normal modes are kept in the dual assembly, the proposed matrix keeps the best convergence among the 4 candidate matrixes.Figure 9
The relative error of mode frequency using 4 methods, and 100 modes were kept for every component.The results of the Craig-Bampton method are not presented in Figure9, it is unfair to compare the two mode synthesis methods using only 100 normal modes for the dual assembly, and the boundary degrees of freedom are kept in the final solution set in the Craig-Bampton method; however, the interface degrees of freedom are condensed in the dual assembly, so how to compare the two methods fairly requires special attention. There are 927 independent degrees of freedom on the boundary, and if p orders fixed-interface normal mode is kept for every substructure, there are 4∗p+927 degrees of freedom in the synthesized system; hence, we must keep p+231 order normal modes in the dual assembly, and the relative error of mode frequency for p=50 is shown in Figure 10.Figure 10
The relative error of mode frequency; 281 norm modes were kept for the dual assembly and 50 norm modes for the primary assembly.It can be inferred from Figure10 that the proposed method preserves comparative accuracy to the Craig-Bampton method. The Craig-Bampton method works well below the 50th order, and the error enlarges quickly over the 50th order, while, for the dual assembly, the accuracy of the first 50 orders is inferior to Craig-Bampton method, and it is superior over 50 orders, as shown in Figure 11.Figure 11
The relative error of mode frequency, convergence comparison between the proposedGs matrix (281 modes each) and the Craig-Bampton method (50 modes for each component and 927 constraint modes).When 281 modes are kept for both methods in every substructure, for Craig-Bampton method, there were 2051 DOFs in the final synthesized model, whereas, for the dual assembly, there were 1124 DOFs in the final synthesized model. It is unfair to implement such comparison; nevertheless, the results are still attached in Figure12. We can find that the lower order modes from the Craig-Bampton method in this case have better accuracy than the dual assembly, and the accuracy is comparable above from the 150th order.Figure 12
The convergence comparison between the proposedGs matrix and the Craig-Bampton method; 281 normal modes were kept for both methods for every component.


## 3.1. A 10-DOF Model

A 10-DOF lumped mass model is shown in Figure1, 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 ∑ns=13 DOFs together with additional c =3 rigid connections.Figure 1
Undamped lumped parameter mass model, only the horizontal dynamics are considered, the global DOFs are labeled in the box, the spring rate K=1E4 N/m, and the unit mass M=1 kg.Figure 2
Mathematically partitioned four substructures; the local DOFs are labeled in the box.Before jumping into details, we show how to build the constraint matrix correctly. There are 4 components sharing one single boundary DOF, therefore, the constraint matrix is nonunique. We can build the displacement compatibility equations asx14=x21, x21=x41, and x33=x41, where the superscript is the DOF number, and the subscript is the substructure number. The obtained constraint matrix is(38)
B
=
b
1
0
0
0
1
0
0
0
0
0
0
0
0
⋮
⋮
⋮
⋮
b
2
-
1
0
0
1
0
0
0
0
0
⋮
⋮
⋮
⋮
b
3
0
0
0
0
0
0
0
0
1
⋮
⋮
⋮
⋮
b
4
0
0
0
-
1
0
0
-
1
0
0We can also build the displacement compatibility equations asx14=x21, x21=x41, x33=x41, and x14=x41, and the following constraint matrix is obtained:(39)
B
=
b
1
0
0
0
1
0
0
0
0
0
0
0
0
0
0
0
1
⋮
⋮
⋮
⋮
⋮
b
2
-
1
0
0
1
0
0
0
0
0
0
0
0
⋮
⋮
⋮
⋮
⋮
b
3
0
0
0
0
0
0
0
0
1
0
0
0
⋮
⋮
⋮
⋮
⋮
b
4
0
0
0
-
1
0
0
-
1
0
0
-
1
0
0We have built two constraint matrixes; the first one is legal while the second one is illegal. The constraint matrix given by formula (38) is full rank in rows, and it is a complete description of the system constraints, while formula (39) has a redundant constraint, or it is rank deficit in rows. If we use the singular constraint matrix in formula (39), we have to invert a singular matrix A=Bdiag⁡GsBT, hence obtaining inferior results.Table1 listed the mode analysis results. A vector was used to denote the total number of the lower order kept normal modes (including rigid mode), and the vector [2,1,1,1] indicates that we kept 2 orders for the first component and 1 order for every other three components, respectively. TheExact (see Table1) refers to the reference result computed from original system using conventional method. TheDCB refers to the Dual-Craig-Bampton method, from (9), and theDCB+Guyan refers to the results from a further static condensation, i.e., from the (22), and the imaginary unit is i=-1.Table 1
The mode frequency comparison in different methods /Hz.
Order Exact DCB DCB+Guyan [ 2,1 , 1,1 ] [ 3,2 , 2,2 ] [ 2,1 , 1,1 ] [ 3,2 , 2,2 ] -3  26.4556i 62.1508i -2  22.8070i 55.1870i -1  20.6766i 45.5169i 1 0.0000 0.0000 0.0000 0.0000 0.0000 2 6.5712 6.5984 6.5715 6.5740 6.5712 3 7.0001 7.0626 7.0007 7.0331 7.0001 4 7.6090 7.7488 7.6098 7.6828 7.6090 5 12.8444 - 12.8509 - 12.8445 6 18.7891 19.6764 18.8169 18.1328 18.7894 7 27.0282  28.3200  27.0636 8 32.0634  47.2122  32.0652 9 33.1482  -  33.6621 10 35.4645  55.3266We can infer from Table1 that spurious modes appear in the DCB method; we marked the spurious modes in bold. Some spurious modes result from the negative eigenvalue, and they are in imaginary frequencies. The 47 Hz and 55 Hz modes are far beyond the highest frequency of the system, and obviously they are spurious modes as well. There is no spurious mode after the proposed static condensation, and the mode frequency accuracy even increases significantly, which shows better convergence property.It is interesting to consider the effects of differentGs. We tried the identity matrix Is, mass matrix Ms, and stiffness matrix Ks in the numerical examples, we listed the corresponding mode frequency results in Table 2, and we can find that there is no spurious mode at least. Essentially, selecting different Gs matrix affects the substructure representation, the interface force balance is weakened physically even though the displacement compatibility is enforced mathematically, and, consequently, these results are inferior to the results from the residual modes.Table 2
The mode frequency comparison using different matrixGs; [3,2,2,2] orders were kept.
Order Exact Gs=Is Gs=Ms Gs=Ks Gs=Rs 1 0.0000 0.0000 0.0000 0.0000 0.0000 2 6.5712 6.5713 6.5713 6.5714 6.5712 3 7.0001 7.0022 7.0022 7.0226 7.0001 4 7.6090 7.6131 7.6131 7.6103 7.6090 5 12.8444 12.8487 12.8487 12.8537 12.8445 6 18.7891 18.7935 18.7935 18.8111 18.7894 7 27.0282 27.1237 27.1237 27.2344 27.0636 8 32.0634 32.0658 32.0658 32.1554 32.0652 9 33.1482 33.7795 33.7795 33.5493 33.6621 10 35.4645 Averaged Relative error (%)  0.30095 0.30095 0.348863 0.211159For further validation, the mode shape accuracy was analyzed. We checked the orthogonality by the MAC matrix analysis [28], this criterion is preferred in modal analysis, and we used [3,2,2,2] kept orders for components in subsequent discussions. The mode shapes were expanded by a series of matrix multiplications. The MAC matrix between exact result and theDCB result was shown in Figure 3, obviously there are three spurious modes in the lowest orders, followed by 7 order modes, and the MAC matrix is diagonally dominated (except for the lowest 3 orders and the highest 2 orders). There are 12 order modes to describe the 10-DOF-system dynamics since the system is augmented by 3 Lagrange multipliers, we have (3+2+2+2) +3=12, and, of course, it is irrational to describe the 10-DOF system with 12 order modes in linear analysis.Figure 3
MAC analysis between accurate shapes and shapes from DCB method.We present the MAC matrix between exact result and theDCB+Guyan result in Figure 4, we can find that no spurious mode exists in the lowest frequency range, 8 orders have good accuracy, and the last order is acceptable; this comparison shows that the static condensation successfully improved the mode shape accuracy.Figure 4
MAC analysis between accurate shapes and shapes from DCB+Guyan.To our knowledge, the MAC analysis is simply a crude understanding of mode shape accuracy, it is insufficient to judge the overall accuracy, therefore, a more rigorous check was followed by the frequency response function (FRF) analysis, the acceleration FRFs were computed, and these FRFs are weighted byω2 in relative to displacement FRFs; therefore they are more sensitive in detecting the shape errors in higher frequency range. For illustrating purpose, we assigned the assumed modal damping 0.01 to every order, and we removed the three lowest order spurious modes in the DCB method; in addition, we removed the identified rigid body mode. Finally, we selected three FRFs, as shown in Figures 5–7.Figure 5
Acceleration FRF comparison, H(2,9).Figure 6
Acceleration FRF comparison, H(5,2).Figure 7
Acceleration FRF comparison, H(6,4).From Figure5, we can find that the exact modes dominate the frequency range within 6~36Hz, the frequency response above the 36Hz is contaminated by the two spurious modes in the DCB method, and they are 47Hz and 55Hz, respectively. The results for the lower frequency analysis were shown in Figures 6 and 7, and 10~20dB improvement was achieved in some frequencies. The above results validated the improvements of the proposed static condensation.


## 3.2. A Bolted Beam Model

The beam model and its mesh are shown in Figure8, respectively; the assembly is composed of 3 components, upper beam, down beam, and the bolt. One end is fixed ideally, and the other is free. The bolt connects the upper part and the lower part, and all the contact interfaces are treated as linear constraint using conforming mesh; i.e., they are glued together. Now we split the assembly into 4 substructures to perform modal synthesis analysis and generate the reference result from the full model analysis.Figure 8
The suspended beam and its mesh model.This beam model is used here to show the convergence property using differentGs and in this example the identity matrix, substructure mass matrix, substructure stiffness matrix, and the proposed matrix were compared. The 4 matrixes share the similar numerical stability and only differ in the convergence, and additional computation effort is required for the proposed matrix. The residual flexibility matrix and the inertia relief attachment modes are not applicable in this example due to ill-conditioning.Figure9 shows that when 100 normal modes are kept in the dual assembly, the proposed matrix keeps the best convergence among the 4 candidate matrixes.Figure 9
The relative error of mode frequency using 4 methods, and 100 modes were kept for every component.The results of the Craig-Bampton method are not presented in Figure9, it is unfair to compare the two mode synthesis methods using only 100 normal modes for the dual assembly, and the boundary degrees of freedom are kept in the final solution set in the Craig-Bampton method; however, the interface degrees of freedom are condensed in the dual assembly, so how to compare the two methods fairly requires special attention. There are 927 independent degrees of freedom on the boundary, and if p orders fixed-interface normal mode is kept for every substructure, there are 4∗p+927 degrees of freedom in the synthesized system; hence, we must keep p+231 order normal modes in the dual assembly, and the relative error of mode frequency for p=50 is shown in Figure 10.Figure 10
The relative error of mode frequency; 281 norm modes were kept for the dual assembly and 50 norm modes for the primary assembly.It can be inferred from Figure10 that the proposed method preserves comparative accuracy to the Craig-Bampton method. The Craig-Bampton method works well below the 50th order, and the error enlarges quickly over the 50th order, while, for the dual assembly, the accuracy of the first 50 orders is inferior to Craig-Bampton method, and it is superior over 50 orders, as shown in Figure 11.Figure 11
The relative error of mode frequency, convergence comparison between the proposedGs matrix (281 modes each) and the Craig-Bampton method (50 modes for each component and 927 constraint modes).When 281 modes are kept for both methods in every substructure, for Craig-Bampton method, there were 2051 DOFs in the final synthesized model, whereas, for the dual assembly, there were 1124 DOFs in the final synthesized model. It is unfair to implement such comparison; nevertheless, the results are still attached in Figure12. We can find that the lower order modes from the Craig-Bampton method in this case have better accuracy than the dual assembly, and the accuracy is comparable above from the 150th order.Figure 12
The convergence comparison between the proposedGs matrix and the Craig-Bampton method; 281 normal modes were kept for both methods for every component.


## 4. Conclusions

This research investigates the dual assembly technique; it is a free-interface mode synthesis technique. The mechanism of negative eigenvalue is caused by the indefinite stiffness matrix, the redundant Lagrange multipliers should be condensed in the modal domain analysis, or the interface condensation is necessary. We applied the projection method to show that a general matrix can be chosen without weakening the displacement compatibility. The dual assembly with residual flexibility matrix or inertial relief attachment modes becomes Craig-Chang’s method, it has an inherent instability caused by numerical ill-conditioning. Finally, a new substructure reduction is proposed based on the projection method, the residual flexibility matrix is replaced by the constraint modes, this matrix is well conditioned, and the accuracy is comparable to the classical Craig-Bampton technique. The bolted beam model had shown that the Craig-Bampton method is superior to the dual assembly, provided that the kept normal modes are the same for every component. If we restrict the final reduced matrixes are of the same dimension, then the dual assembly with the proposed matrix can have better accuracy in the higher order modes, while the Craig-Bampton method still keeps a better convergence in the lowest orders. This contribution is important in applications where the free-interface is necessary and the residual flexibility matrix fails; in this case, the suggested substructure reduction can be used to ensure the numerical stability and to keep acceptable accuracy.



---

*Source: 1589794-2018-07-17.xml*