Table of Contents Author Guidelines Submit a Manuscript
Shock and Vibration
Volume 2018, Article ID 1589794, 12 pages
Research Article

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.


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 [46]. 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 [1518], 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 [2527] 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.

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 as , , and , where the superscript is the DOF number, and the subscript is the substructure number. The obtained constraint matrix is

We can also build the displacement compatibility equations as , , , and , and the following constraint matrix is obtained:

We 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 , hence obtaining inferior results.

Table 1 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 indicates that we kept 2 orders for the first component and 1 order for every other three components, respectively. The Exact (see Table 1) refers to the reference result computed from original system using conventional method. The DCB refers to the Dual-Craig-Bampton method, from (9), and the DCB+Guyan refers to the results from a further static condensation, i.e., from the (22), and the imaginary unit is .

Table 1: The mode frequency comparison in different methods /Hz.

We can infer from Table 1 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 different . We tried the identity matrix , mass matrix , and stiffness matrix 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 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 matrix Gs; orders were kept.

For 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 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 the DCB 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 the DCB+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 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 57.

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 Figure 5, 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 Figure 8, 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 different 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.

Figure 9 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 Figure 9, 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 4p+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 Figure 10 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 proposed Gs 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 Figure 12. 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 proposed Gs 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.

Data Availability

The stiffness and mass matrix data in the first example to support the findings of this study is within the article, and the data in MATLAB format in the second example can be available on request, together with all the MATLAB files we used. They can be consulted from the email address: (Mou Xiaolong).

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This research is supported by the Programme of Introducing Talents of Discipline to Universities of China (Grant no. B12022) and the National Nature Science Foundation of China (Grant no. 51006010). The authors would like to thank the sponsors.


  1. J. Philippe, F. Thouverez, L. Blanc, and M. Gruin, “Vibratory behavior prediction of mistuned stator vane clusters: An industrial application,” Computers & Structures, vol. 196, pp. 12–23, 2018. View at Publisher · View at Google Scholar · View at Scopus
  2. R. Craig and M. Bampton, “Coupling of substructures for dynamic analyses,” AIAA Journal, vol. 6, no. 7, pp. 1313–1319, 1968. View at Publisher · View at Google Scholar · View at Scopus
  3. R. R. Craig and A. J. Kurdila, Fundamentals of Structural Dynamics, John Wiley, 2004.
  4. J.-H. Kim, J. Kim, and P.-S. Lee, “Improving the accuracy of the dual Craig-Bampton method,” Computers & Structures, vol. 191, pp. 22–32, 2017. View at Publisher · View at Google Scholar · View at Scopus
  5. S.-H. Boo, J.-H. Kim, and P.-S. Lee, “Towards improving the enhanced Craig-Bampton method,” Computers & Structures, vol. 196, pp. 63–75, 2018. View at Google Scholar · View at Scopus
  6. S.-H. Boo and P.-S. Lee, “An iterative algebraic dynamic condensation method and its performance,” Computers & Structures, vol. 182, pp. 419–429, 2017. View at Publisher · View at Google Scholar · View at Scopus
  7. Y.-S. Kim and H.-C. Eun, “Reanalysis of modified structures by adding or removing substructures,” Advances in Civil Engineering, vol. 2018, Article ID 3084078, pp. 1–9, 2018. View at Publisher · View at Google Scholar
  8. R. R. Craig and C.-J. Chang, “Free-interface methods of substructure coupling for dynamic analysis,” AIAA Journal, vol. 14, no. 11, pp. 1633–1635, 1976. View at Publisher · View at Google Scholar · View at Scopus
  9. R. R. Craig, “Coupling of substructures for dynamic analyses - An overview,” AIAA Journal, pp. 1–12, 2000. View at Publisher · View at Google Scholar
  10. R. H. MacNeal, “A hybrid method of component mode synthesis,” Computers & Structures, vol. 1, no. 4, pp. 581–601, 1971. View at Publisher · View at Google Scholar · View at Scopus
  11. S. Rubin, “Improved component-mode representation for structural dynamic analysis,” AIAA Journal, vol. 13, no. 8, pp. 995–1006, 1975. View at Publisher · View at Google Scholar · View at Scopus
  12. D. R. Martinez and D. L. Gregory, A Comparison of Free Component Mode Synthesis Techniques Using MSC/Nastran.No.SAND83-0025, Sandia National Laboratories, Albuquerque, NM, USA, June 1984.
  13. R. R. Craig Jr and C.-J. Chang, “Substructure coupling for dynamic analysis and testing,” NASA Report, CR-2781, February 1977. View at Google Scholar
  14. N. Bouhaddi and J. P. Lombard, “Improved free-interface substructures representation method,” Computers & Structures, vol. 77, no. 3, pp. 269–283, 2000. View at Publisher · View at Google Scholar · View at Scopus
  15. P. L. C. van der Valk, Model Reduction & Interface Modeling in Dynamic Substructuring-Application to a Multi-Megawatt Wind Turbine, 2010.
  16. D. J. Rixen, “A dual Craig-Bampton method for dynamic substructuring,” Journal of Computational and Applied Mathematics, vol. 168, no. 1-2, pp. 383–391, 2004. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  17. D. De Klerk, D. J. Rixen, and S. N. Voormeeren, “General framework for dynamic substructuring: History, review, and classification of techniques,” AIAA Journal, vol. 46, no. 5, pp. 1169–1181, 2008. View at Publisher · View at Google Scholar · View at Scopus
  18. D. J. Rixen, Dual Craig-Bampton with Enrichment to Avoid Spurious Modes, vol. 27th-IMAC, Florida, USA, 2009.
  19. T. Miyamura, “Incorporation of multipoint constraints into the balancing domain decomposition method and its parallel implementation,” International Journal for Numerical Methods in Engineering, vol. 69, no. 2, pp. 326–346, 2007. View at Publisher · View at Google Scholar · View at Scopus
  20. B. Case, “User reference manual for the mystran general purpose finite element structural analysis computer program 2011,” MYSTRAN Software User Manual, 2015,
  21. O. C. Zienkiewicz, R. L. Taylor, and D. D. Fox, The Finite Element Method for Solids & Structural Mechanics, Elsevier, 7th edition, 2014.
  22. K. J. Bathe, Finite Element Procedures, 2nd edition, 2014.
  23. M. Géradin and A. Cardona, “Numerical integration of second-order differential-algebraic systems in flexible mechanics dynamics,” in Computer-Aided Analysis of Rigid and Flexible Mechanical Systems, M. Pereira and J. Ambrosio, Eds., vol. 268 of NATO ASI Series. Applied Sciences, pp. 233–284, Kluwer, Dordrecht, Netherlands, 1993. View at Google Scholar · View at MathSciNet
  24. P. Fisette and B. Vaneghem, “Numerical integration of multibody system dynamic equations using the coordinate partitioning method in an implicit Newmark scheme,” Computer Methods Applied Mechanics and Engineering, vol. 135, no. 1-2, pp. 85–105, 1996. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  25. P. Lidström and P. Olsson, “On the natural vibrations of linear structures with constraints,” Journal of Sound and Vibration, vol. 301, no. 1-2, pp. 341–354, 2007. View at Publisher · View at Google Scholar · View at MathSciNet
  26. W. D'Ambrogio and A. Fregolent, “The role of interface DoFs in decoupling of substructures based on the dual domain decomposition,” Mechanical Systems and Signal Processing, vol. 24, no. 7, pp. 2035–2048, 2010. View at Publisher · View at Google Scholar · View at Scopus
  27. M. Ainsworth, “Essential boundary conditions and multi-point constraints in finite element analysis,” Computer Methods Applied Mechanics and Engineering, vol. 190, no. 48, pp. 6323–6339, 2001. View at Publisher · View at Google Scholar · View at MathSciNet · View at Scopus
  28. M. Xiao-long and F. Hui-hua, Discussions about Projection Methods in Dealing with Linear Nonhomogeneous Multipoint Constraints, NCVTA, Nanning, China, 2017.