Preconditioning strategies for linear dependent generalized component modes in 3D flexible multibody dynamics

The trend away from physical towards virtual prototyping as well as increasing industrial demands require advanced simulation tools for dynamical systems. Virtually all engineering systems are assemblies and are associated with stresses, noise and vibrations; therefore, flexible multibody simulations are inevitable for accurate predictions. However, real-world finite element models contain millions of degrees of freedom that cannot be reasonably handled without model reduction techniques. Generalized component mode synthesis is a promising tool for flexible 3D multibody systems, since the generalized component modes not only represent the deformation modes in any possible orientation, but also rigid body motion, which preserves a linear configuration space, yielding a constant mass matrix, a co-rotated but constant stiffness matrix, no quadratic velocity vector and a simple structure of the equations of motion. In this novel framework, the displacement is approximated by a linear combination of generalized component modes generated from vibration eigenmodes, undeformed nodal coordinates and the Cartesian base vectors. The emerging system matrices may be ill-conditioned and may introduce significant numerical errors, because of linearly dependent generalized component modes and due to different orders of magnitude of their Euclidean norm. However, this issue has not received much attention in the open literature despite its importance. Hence, the current contribution sheds light on this problem and derives preprocessing procedures to convert ill-conditioned into well-conditioned problems, which shall improve the formulation’s applicability. The new findings are illustrated by numerical experiments of simple bodies and a crankshaft.


Introduction
The trend away from physical towards virtual prototyping as well as increasing industrial demands on reliability and efficiency of modern dynamical engineering devices require advanced modeling techniques during the design process. At present, virtually all such engineering systems are assemblies and, hence, made out of multiple components, which interact with each other during operation. The forces required to execute desired motions are associated with stresses, noise and vibrations of the bodies in the system. Thus, it is insufficient to model multibody systems as rigid bodies and extract boundary forces to perform subsequent standard finite element (FE) analyses. Flexible multibody simulations, where the system is spatially discretized using a finite number of elements, are, therefore, inevitable. However, most FE models of relevant engineering problems contain a huge number of degrees of freedom (DOFs) that cannot be simulated within a reasonable amount of time without employing model reduction techniques.
There are several ways to model flexible multibody systems [19,24]. Among them the absolute coordinate formulation known as generalized component mode synthesis [16], which is a promising alternative to well-established flexible multibody formalisms, such as the floating frame of reference formulation (FFRF) [15,20,21], since the so-called generalized component modes not only describe rigid body motion, but can also represent the deformation modes in any possible orientation, which preserves a linear relationship between the global displacement field and the DOFs of the considered domain; this yields a constant mass matrix, a constant but co-rotated stiffness matrix with one co-rotational frame for each system body only, a trivial quadratic velocity vector and a simple structure of the governing equations. The mass and stiffness matrix are simply the standard linear FE system matrices pre-and post-multiplied with a constant reduction matrix. Hence, the FE system matrices are extracted only once during preprocessing and remain constant during the whole simulation, which is why, the algorithm is fully decoupled from any FE package and easily applicableit requires just a few lines of code -to any multibody system subjected to large reference motion, but small deformations of the individual components.
Both the FFRF and the generalized component mode synthesis approximate the flexible deformation by a linear combination of component modes. In case of the FFRF, the component modes are, e.g., the eigenmodes of vibration limited to the frequency range of interest, or a combination of eigenmodes and static modes, see, e.g., the pioneering work of [1,9,10,12,17]. If the deformation is approximated by vibration modes, the reduction matrix containing column-wise the eigenmodes is well-conditioned, since the eigenvectors are linearly independent [2, p. 158] even for repeated eigenvalues -the degeneracy theorem [4, p. 72] allows the generalization of the orthogonality of the eigenvectors in the metrics of the FE mass and stiffness matrices in the case of repeated eigenvalues. Hence, the reduction matrix of the FFRF-based component mode synthesis does not introduce numerical errors; in fact, the condition number is close to one if the eigenvectors are displacement normalized, e.g., smaller than 1.4 for all hereinafter analyzed models, whereas the generalized component mode reduction matrix is in many cases ill-conditioned, e.g., on the order of 10 6 to 10 17 for the hereinafter analyzed beam-like models, because of linear dependencies and Euclidean norms that deviate by orders of magnitude from each other between the generalized component modes, which may arise due to the special structure of the reduction basis Φ, see Sect. 2.1. It is well known that linear dependencies and high condition numbers may lead to large errors or even unsolvable problems, since they preclude the factorization of the system Jacobian impossible.
There is past work concerned with a continuum-mechanics-based mathematical derivation of the generalized component modes and equations of motion [5,16], whereas [6] shows how the idea may be employed without model reduction. It is also known from the literature that the inherent properties of the formulation facilitate the construction of energymomentum-conserving time integration schemes [7]. Also, the formulation has been successfully applied to engineering problems, such as fluid-structure interaction [18] and machine parts with large rotations about one axis only [26], and was extended by means of a nullspace projection approach [25] as well as the global modal parametrization [3,8,13]. However, the problem of linear dependencies and ill-conditioned system matrices has not received much attention despite its importance. The issue was marginally reported in [5,16], but has not been addressed in the available literature. Furthermore, it has been believed that these issues can only arise for symmetric problems, which is, as shown in the present paper, in general not true. Hence, the main objective of this contribution is to analyze this problem associated with the generalized component mode synthesis and to present preprocessing strategies to obtain well-conditioned problems.
The remaining part of the paper is organized as follows: Sect. 2 is devoted to the configuration space -the generalized component modes are derived on a nodal-based level and the linear relationship between the global nodal displacements and DOFs is illustrated. Furthermore, the section contains a traceable and more intuitive derivation of the governing equations of motion with a Lagrangian formulation for a general spatially discretized mechanical system, in contrast to the original continuum-mechanics-based derivation reported in the literature. This novel presentation shall also clarify the idea behind the method to some extent. Section 3 presents a short note on the mathematics of linear dependencies and how to handle them, which is required to analyze and eliminate linear dependencies in the arising reduction matrix; followed by an illustrative presentation of the generalized component modes and the problem of linear dependencies in Sect. 4. Sections 5 and 6 investigate the linear dependencies of the flexible part of the reduction matrix and between the flexible and rigid body motion parts, respectively, with the help of simple FE beamlike 3D models. Section 7 applies the presented theory to an FE model of a crankshaft of a two-cylinder reciprocal combustion engine, and Sect. 8 gives step-by-step preprocessing strategies for handling the inherent problem of the formulation to convert ill-conditioned into well-conditioned problems; followed by a conclusion in Sect. 9.
The present paper is a revised and extended version of the conference paper [27] presented at the Fifth Joint International Conference on Multibody System Dynamics (IMSD) 2018.

Reduction basis
The idea behind the generalized component mode formalism is (i) to define a reduction basis that can represent large rigid body translation and rotation, as well as flexible deformation; and (ii) to obtain a linear relationship between the global FE nodal displacements and the DOFs. As already mentioned, the formulation exploits a modal superposition reduction method, where the flexible deformation is approximated by a linear combination of vibration modes, to reduce the system size from a large number of DOFs to a significantly smaller one. This, of course, restricts its applicability to problems where the deformation of the components remains linear within each body frame.
The so-called generalized component modes not only represent the deformation modes in any possible orientation, but also account for large rigid body motion (translation & rotation), yielding a linear configuration space; see Eq. (22). The generalized component modes Fig. 1 Spatially discretized domain with global F and body fixed F frame, related by the rotation matrix A ∈ R 3×3 . The position vector r (i) = x (i) + c (i) defines the current position of node i, with corresponding undeformed nodal coordinates x (i) and nodal displacement c (i) are generated from the set of Cartesian unit vectors, the undeformed FE nodal coordinates and the original eigenmodes of vibration in the frequency range of interest. However, the linear relationship between the global displacements and the DOFs, i.e., the generalized coordinates q, is obtained at the expense of a nine-fold increase in the number of flexible modal coordinates m ζ , i.e., nine generalized flexible coordinates per original natural mode of vibration m ψ of the bodies in the system, as shown in the following paragraphs.
The formal rules to generate the generalized translational, rotational and flexible component modes are explicitly stated in Eqs. (8), (13) and (20), respectively, according to [16]. It is shown in [16] that these modes may be obtained by rearranging the equations during the derivation of the FFRF-based component mode synthesis and by introducing new coordinates. However, a more intuitive and concise derivation of the generalized component modes is shown in this contribution.
To obtain the desired reduction basis, the global FE nodal displacement c is split into its translational c t , rotational c r and flexible c f part, i.e., The column matrix c contains the nodal displacements of all FE nodes and is related to the nodal coordinates in the reference configuration x and to the current positions of the nodes r via (Fig. 1) where r = ⎡ ⎢ ⎣ r (1) . . .
contain the corresponding n n nodal vectors. 1 All nodes share the same displacement for a translational rigid body motion. Hence, any rigid body translation of the FE nodes may be represented as where e l ∈ R 3 (l = 1, 2, 3) denote the orthonormal set of Cartesian base vectors, I ∈ R 3×3 the identity matrix, Φ t ∈ R 3nn×3 the translational reduction matrix containing column-wise the generalized translational component modes and q t = [q t1 q t2 q t3 ] T the corresponding generalized translational coordinates. Likewise, any rigid body rotation of the FE nodes of a body in the system may be represented as where A bd = diag(A, . . . , A) ∈ R 3nn×3nn and I bd = diag(I , . . . , I ) ∈ R 3nn×3nn denote a block-diagonal matrix with the rotation matrix A ∈ R 3×3 (relating the global F and local F frame, see Fig. 1) and the identity matrix on its diagonal, respectively, x (i) k the ith nodal coordinate of the kth coordinate direction of the undeformed FE model, Φ r ∈ R 3nn×9 the rotational reduction matrix containing column-wise the generalized rotational component modes and q r the corresponding generalized rotational coordinates. Equation (10) is nothing but the standard matrix multiplication, where the involved components are rearranged in new matrices; this restructuring will be also employed in the next paragraph to derive the flexible reduction matrix, see Eq. (16). Finally, the flexible deformation is, as already mentioned, approximated by a linear combination of vibration eigenmodes, i.e., calculated in each body frame F (Fig. 1). Hence, the flexible nodal displacement in the global frame F reads where m ψ (i) k denotes the displacement of the ith node in the kth coordinate direction of the mth natural eigenmode of vibration, Φ f ∈ R 3nn×9nm the flexible reduction matrix containing column-wise the generalized flexible component modes q f the corresponding generalized flexible coordinates and n m the number of vibration eigenmodes included in the reduction basis. Hence, combining Eqs. (1), (7), (12) and (19) yields which gives the desired linear relationship between the full set of global FE nodal displacements c and the reduced set of generalized coordinates q, i.e., where dim(. . . ) denotes the dimension of the quantity in parentheses; the reduction matrix Φ obviously contains the translational Φ t ∈ R 3nn×3 , the rotational Φ r ∈ R 3nn×9 and the flexible Φ f ∈ R 3nn×9nm generalized component modes, i.e., where the order of appearance of the generalized translational, Eq. (8), rotational, Eq. (13), and flexible, Eq. (20), component modes in Φ t , Φ r or Φ f , respectively, is given by if i represents the column index of Φ t , Φ r or Φ f . The generalized coordinates q are, of course, also partitioned in a similar manner as the reduction matrix Φ, see Eq. (23), i.e., To illustrate the generalized component mode generation procedure, consider a simple hypothetical example of one eigenmode of an FE model with only one node and three DOFs, Eqs. (28)

Equations of motion
In this section the governing equations of motion, originally derived on the basis of the FFRF via a continuum mechanics approach, are derived in a compact and traceable way, starting in the already spatially discretized domain. The semi-discretized equations of motion are derived via Lagrange's equation, i.e., d dt with the modified Lagrangian L [11] for a general mechanical system defined as where W denotes the work done by applied nodal forces, λ the vector of Lagrange multipliers and g(c) = 0 the vector of holonomic constraint equations; the product λ T g represents the potential energy of the forces that maintain kinematical constraints [11]. The strain energy V , kinetic energy T and W written on a semi-discrete level read where c flx denotes the flexible FE nodal displacements in the body frame, K and M the linear FE stiffness and mass matrix, respectively, and f the applied FE nodal forces. Note that Eq. (31) is stated in terms of c flx , i.e., only the flexible part of q contributes to V -rigid body motion cannot give rise to strains, however, it contributes, of course, to the kinetic energy and work done by applied forces. The generalized coordinates q can be split, i.e., into two terms associated with the rigid body displacement q rig and the flexible deformation q flx , given by respectively, where q * r denotes proper rotational coordinates obtained from q r by an appropriate orthonormalization procedure and is therefore a function of q r , i.e., q * r = q * r (q r ) and no additional DOFs are introduced here. Note that the rotational q r , Eq. (12), and flexible q f , Eq. (19), coordinates do not correspond directly to rotational and flexible displacements only, since q r may also contain stretch and shear deformations, see Appendix A. The displacements calculated by Eq. (12), with the rotational generalized component modes Φ r , exclusively represent rigid body rotation, only if the rotational coordinates q r are extracted from an orthogonal matrix, i.e., if, see Eqs. (10) to (12), where If the columns of Q r and Q * r are (implicitly) defined as, see also Eqs. (9) to (12) and Eq. (37), respectively, the proper (orthonormal) rotational DOFs q * r are obtained from the orthonormal 3 } through an appropriate orthonormalization procedure, see, e.g., [16].
The original generalized component mode synthesis formulation does not enforce the orthogonality condition, Eq. (36), during simulation, i.e., Q r + I is in general not orthogonal, and q r may also contain stretch and shear deformations, see Appendix A. Hence, the proper global flexible nodal displacement c flx , with contributions from q r and q f , reads and is related to the local flexible FE nodal displacements via as may be seen from Fig. 1. Therefore, from Eq. (40) and Eq. (41), since The matrices in Eq. (43) commute, 2 since the reduction matrix Φ consists of blocks that are multiples of the 3 × 3 identity matrix, see Sect. 2.1. Substituting Eqs. (22), (34) and (42) into Eqs. (31) to (33) yields the strain and kinetic energy, as well as the work done by applied nodal forces in terms of the generalized coordinates, i.e., The equations of motion are then obtained by combining Eqs. (22), (29), (30), (35), as well as Eqs. (44) to (46), and carrying out differentiation with respect to q,q and t ; the individual contributions read d dt with the reduced system matrices where G is the standard constraint Jacobian obtained by differentiating the constraints g, defined in terms of the physical FE nodal DOFs c, with respect to c. The third term on the right hand side of Eq. (48) vanishes, since the virtual elastic work done by rigid body displacements is zero -rigid body displacements cannot give rise to elastic forces, i.e., which must hold for any virtual displacement; hence, which is also verified by numerical experiments. 3 Summing up the contributions of Eqs. (47) to (50) according to Eqs. (29) and (30), yields the governing equations of motion, namely where the nonlinear term due to the strain energy is given by 3 A short note on linear dependence

Definition
As already pointed out in Sect. 1, a straightforward generation of Φ according to Eqs. (8), (13), (20) and Eq. (23) may lead to an ill-conditioned system due to linearly dependent generalized component modes, which is why this section is devoted to a short note on the mathematics of linear dependence. The generalized component modes ϕ i , i.e., the columns of Φ, are said to be linearly dependent if if the scalars a pi are not all zero; this may be written in matrix notation as 4 where Hence, the generalized component modes are linearly independent, if the nullspace of the reduction matrix Φ contains solely the zero vector, i.e., if null(Φ) = {0}.

Singular-value decomposition
To reveal information about the nature of the linear dependencies, the singular value decomposition (SVD) of the real-valued matrix Φ, reading [23, p. 25] may be employed to determine the nullspace of Φ. In Eq. (63), U ∈ R 3nn×3nn and V ∈ R (12+9nm)×(12+9nm) are orthogonal 5 matrices whose columns are the left u i and right 3 Note that ∂q rig /∂q was calculated analytically, see Eq. (35), where q * r was obtained by Gram-Schmidt orthonormalization from q r ; A bd is then given by Eq. (38). Equation (57) was then verified for the beam-like models of Sect. 5 at different times during simulation and for arbitrary q's. 4 Note that the subscript p should indicate that the vector a is not unique. v i singular vectors, respectively, and Σ ∈ R 3nn×(12+9nm) is a rectangular diagonal matrix with the singular values σ i on its diagonal.
Right multiplying Eq. (63) with V yields with 12 + 9n m ≤ 3n n . Equation (66) shows that the right singular vectors v p corresponding to vanishing singular values σ p are elements of the nullspace of Φ. Hence, the components of v p are the coefficients a pi of Eq. (60), i.e., v p = a p , compare Eq. (61) with Eq. (66), and we have identified the nature of the dependencies.

Matrix condition number
The SVD is not only useful to determine the nullspace of a matrix, but may be also used to determine the condition number of a matrix, cond(. . . ), which defines how accurately a linear system of equations can be solved. The condition number may be considered as an amplification of error and one has to expect to lose log 10 [cond(J )], see Eq. (69), digits of accuracy on top of all the other errors, such as finite machine precision [23, p. 95]. If the condition number is small, the matrix is said to be well-conditioned, otherwise ill-conditioned.
For a rectangular matrix, such as Φ ∈ R 3nn×(12+9nm) with 12 + 9n m ≤ 3n n , the condition number is in general defined as [23, p. 95] where Φ + denotes the pseudoinverse; if . . . = . . . 2 , Eq. (67) may be written as the ratio of the largest σ max to the smallest σ min singular value, i.e., which links the SVD and the matrix condition number; and we have obtained a quality measure and an indicator of linear dependence of Φ.
If the governing equations of motion, Eq. (58), are integrated with the Newmark method [14], the inverse of the Jacobian reads [5] with the Newmark parameters β = 1 4 and γ = 1 2 according to [5] and the time increment τ . The condition of the reduction matrix Φ is related to the condition of the Jacobian J , which can be estimated by means of matrix norm relations [23, p. 17] and Eq. (67), but is also evident from numerical experiments, i.e., For example, the condition number of the part of the Jacobian that needs to be factorized during time integration, see Eq. (69), is between four to eight orders of magnitude higher than the condition number of the corresponding ill-conditioned reduction matrices for the hereinafter analyzed models of Sects. 5 and 6, which is why it is essential to ensure that cond(Φ) is sufficiently small to prevent errors.

Cosine similarity
A special case of linear dependence appears whenever two columns of Φ are directly proportional to each other, i.e., if ϕ i = bϕ j for some nonzero scalar b. Directly proportional columns fulfill the condition where the row i and column j indices of Φ correspond to the indices k, l and m of Eqs. (8), (13) and (20) in same manner as before, see Eq. (24), and is a measure of similarity between two vectors and known as the Cosine similarity. A value of ±1 indicates linear dependence, whereas 0 shows orthogonality.

Orthogonal reduction basis
Given a set of linearly dependent generalized component modes {ϕ 1 , . . . , ϕ n } and applying a "shortened" version of the classical Gram-Schmidt process, i.e., where the normalization steps are omitted, generates a set orthogonal generalized component modes {ϕ 1 , . . . , ϕ ν } and a set of zero vectors {0 1 , . . . , 0 (n−ν) }; the algorithm outputs the zero vector in the κth step if for some nonzero scalars b λ , since substituting Eq. (74) into the inner product in the numerator of Eq. (73) and noting that 6 yields The number of nonzero vectors generated by the algorithm is equal to the dimension of the space spanned by the original set of linearly dependent generalized component modes; hence, all zero vectors must be discarded from the reduction basis to obtain the desired independent set of generalized component modes.

Remark
It should be noted that the condition number of the translational reduction matrix Φ t is always equal to one, since the translational modes are always orthogonal to each other, Eq. (8), and that it is neither possible to remove nor alter the rotational part Φ r without changing the formulation. Hence, if linear dependencies arise, only generalized flexible modes may be removed or altered without the need to adopt the formulation. However, it is known from the literature [22] that scaling, i.e., adopting the Euclidean norm of modes to similar values, has a positive effect on the condition of a reduction basis, which is possible for both translational and flexible modes, and will be also employed in here, see Sect. 5.2.
Also, as evident from Sect. 2.1, the reduction basis Φ ∈ R 3nn×(12+9nm) contains predominately generalized flexible component modes, which is why the main part of the rest of the present paper is devoted to the investigation of the flexible part of the reduction matrix, since Φ f ∈ R 3nn×9nm accounts for the majority of linear dependencies and the high condition number in the first place.

An illustrative example
In this section, a simple example of a 60 mm square-sectioned 900 mm long steel beam should illustrate the generalized component modes, as well as the inherent problem of linear dependencies, to gain a deeper understanding of the formulation and to show the significance of the current contribution, respectively.  Finally, Fig. 4 depicts the first two bending eigenmodes of vibration with their nine corresponding generalized flexible component modes. It is verified numerically (see Fig. 5), but may be also seen in the figure that the first three generalized component modes of the first and second bending eigenmodes are identical, i.e., Moreover, the fourth to sixth generalized component modes of these two eigenmodes are related by the factor of negative one, i.e., and are, therefore, directly proportional. Likewise, the set of the first three generalized component modes within each bending eigenmode is directly proportional to the set of the second three generalized component modes, i.e., m φ f1l ∝ m φ f2l for m = 1, 2 and l = 1, 2, 3, which shows that it is also possible for linear dependencies to arise within one generalized component mode set. Equations (77) to (79) also imply that These linear dependencies, Eqs. (77) to (81), are attributed to the symmetric crosssection of the analyzed beam and manifest themselves in a high condition number, cond(Φ) = 9.51 × 10 17 . Note that the first two eigenmodes of vibration form a repeated mode pair with equal eigenfrequencies, where the displacements in the x-and y-directions are simply "exchanged" for the first and second bending mode, since any axis through the centroid of a square cross-section is a principal bending axis. Consequently, only nine out of the full set of 18 generalized flexible component modes are, in this case, linearly independent and the system of equations with the initially 9n m generalized flexible coordinates may be further reduced significantly. The degree of linear dependence of the generalized flexible component modes is dictated by the value of cos θ ij , see Eq. (71) and Sect. 3.4. The Cosine similarity may be represented as a two-dimensional array of numbers and visualized with a color array, see, e.g., the plot of | cos θ ij | of the 18 generalized component modes (Fig. 4) visualized in Fig. 5. The Cosine similarity matrix is symmetric, which is why it is represented with a lower triangular array, hereinafter referred to as heatmap. Every entry of the heatmap represents the absolute value of the Cosine similarity between two generalized flexible component modes; the diagonal entries represent the similarity between a mode with itself and are, therefore, of course, one. The first row and the first column correspond to 1 φ f11 ; increasing the row or column index of the heatmap corresponds to an increase of the indices of m φ fkl in the following manner, see also Eq. This example shows the importance of the mode selection process. The systematic investigation of the generalized component modes, which has not been addressed in the available literature, is not only required to obtain an accurate solution or a solvable system of equa-  Fig. 4; displayed as a lower triangular matrix due to symmetry (Color figure online) tions, but may also enable a further reduction of the generalized coordinates and, therefore, a gain in efficiency.

Detailed analysis of a square-sectioned extruded body
It is clear from the introductory example of Sect. 4 that some linear dependencies between generalized flexible component modes are obvious, for example, it is intuitive that linear dependencies between the generalized component modes obtained from the repeated mode pair (Fig. 4, (a) and (k)) exist, since the displacements of the nodes in the x-and y-directions of 1 ψ and 2 ψ are simply "exchanged". However, this is in general, of course, not the case. This section exhibits that linear dependencies of generalized component modes cannot be identified a priori. To this end, a flexible reduction basis is generated from a selection of eight eigenmodes of a 100 mm square-sectioned beam; all the beam-like models subsequently discussed share the following properties [25]: Young's modulus, E = 1500 MPa Poisson's ratio, μ = 0.3 Density, ρ = 1000 kg m −3 Length, l = 2 m Twenty-noded hexahedrals and 15-noded wedge elements (if necessary) were used for the FE discretization. The flexible reduction basis chosen for the analysis includes the first four bending modes, i.e., the first two bending mode shapes in two perpendicular (⊥) directions (hereinafter distinguished with an asterisk*), as well as the first two torsion and longitudinal eigenmodes, see Fig. 6, which leads to 8 × 9 = 72 generalized flexible component modes contained in Φ f with cond(Φ f ) = 2.97 × 10 6 . In contrast, the condition number of the matrix containing the eight original displacement normalized eigenmodes of vibration (Fig. 6) obtained by an FE solver is equal to 1.00. The high condition number of Φ f indicates that the generalized flexible component modes are not independent, which may be also seen in the Cosine similarity matrix visualized in Fig. 7; the figure depicts the absolute value of the Cosine similarity, Eq. (71), of two generalized component modes in turn. The last 18 columns, corresponding to the repeated mode pair of the first bending eigenmode, show, of course, the same pattern and therefore the same dependencies as Fig. 5, i.e., Eqs. (77) to (81).  Fig. 6 of a square-sectioned beam. The first nine columns and rows correspond to the eigenmode depicted in Fig. 6(a), the following nine to the eigenmode of Fig. 6(b), and so forth, see Eq. (24) (Color figure online) Not surprisingly, the repeated mode pair of the second bending mode shape, Fig. 6(e)-(f), shows the same pattern as the first one, i.e., compare columns 37 to 54 with 55 to 72 in Fig. 7. Figure 7 also reveals that it is possible for linear dependencies between different kinds of modes to arise. Columns 28 to 36 indicate that there is a strong dependence (| cos θ ij | = After the linear dependencies are identified, it is important to deal with them appropriately. One way to improve the condition number of the flexible reduction basis is to exclude one of the generalized component modes involved in each correlation. This has the advantage that the set of generalized coordinates is reduced by the number of excluded modes, which leads to a further improvement of the efficiency of the formulation. The only drawback is that one has to determine a threshold value of the Cosine similarity c th , due to the fact that we are dealing with numerical vectors only and a value of exactly one cannot be expected. The generalized component modes are considered to be linearly dependent if | cos θ ij | ≥ c th . This may be an iterative process that continues as long as the condition number is as low as desired, which depends on the problem to be analyzed and on the system matrices, or until removing correlated columns does not improve the condition number further (Fig. 8). Figure 8(a) shows the matrix condition number of the flexible reduction matrix cond(Φ f ) generated from the eight eigenmodes of vibration displayed in Fig. 6 over the threshold value c th of the Cosine similarity. Figure 8(b) depicts the corresponding number of removed columns of Φ f over c th . The figures indicate that a threshold value relatively close to one (c th = 0.999) suffices to reduce the condition number by three orders of magnitude; such a sharp drop can be generally expected if columns are directly proportional, as evident from all analyzed beam-like models discussed in this paper. Further reducing the threshold (c th = 0.995) leads to a minor improvement by a factor of approximately two. After that, the condition number remains constant within the considered range and removing more columns (from 24 to 27 at c th = 0.980) does not lower the condition number further. Hence, if the Cosine similarity is employed to eliminate linear dependencies, "convergence" of cond(Φ f ) should always be checked prior time integration to avoid numerical errors. Similar numerical experiments on the influence of weak "dependencies" on the reduction matrix's condition number were conducted for all herein analyzed models, see Sect. 5.2, and indicated that a Cosine similarity (absolute) value below 0.950 may be considered as extremely low. In fact, a value of c th = 0.993 yielded sound results throughout the investigations presented in this paper, 7 which is why only entries with | cos θ ij | ≥ 0.993 in the Cosine similarity matrix are further discussed in here.
The initial threshold value may be chosen to be one, or as the Cosine similarity (absolute) value calculated from a pair of generalized component modes where linear dependence is a priori known, which is the case between the generalized component modes obtained from B1 ψ and B1 * ψ (Fig. 6) for the here considered example. Another and a very elegant possibility to eliminate the linear dependencies within Φ f is to apply the "shortened" version of the Gram-Schmidt process, introduced in Sect. 3.5, to the dependent set of modes. This provides an (virtually) orthogonal flexible reduction basis if the arising zero vectors are omitted, as already described in Sect. 3.5. This has the advantage that no threshold value needs to be determined, since the norm of the arising (virtually) zero vectors is by several orders of magnitude smaller than the average norm of the full set of the arising vectors. Another advantage is that the algorithm is only required once, i.e., no iterations, and no user interaction is required at all to obtain a perfectly well-conditioned flexible reduction basis. This algorithm is also faster, not only because no iterations are required, but also since no SVD needs to be performed to estimate the condition number, since the condition of the generated "nonzero set" is a priori sufficiently small. The only drawback is that the arising reduction basis is in general larger than that obtained with a neat Cosine similarity preconditioning. Hence, there is a trade-off between preprocessing and "online" cost, however, the "shortened" Gram-Schmidt process (Sect. 3.5) should be, in general, preferred, due to dispensable user intervention.

Extruded bodies with different cross-sections
The analysis described in Sect. 5.1 was also conducted for bodies with the same properties, as outlined in the first paragraph of the section, and with the same eigenmodes included in the reduction basis, see Fig. 6, but with different cross-sections, see Fig. 9. The main findings are summarized in the following paragraphs and in Table 1.
The condition number of the flexible reduction matrices of the analyzed beams ranges from the order of 10 5 to 10 11 , where the "best" condition number belongs to the arbitrary beam, Fig. 9(e), and the "worst" one to the circular beam, Fig. 9(c). This clearly emphasizes the importance of the present investigations, since one has to expect to lose at least eleven digits of accuracy on top of all the other errors, if the linear dependencies are not handled appropriately.
Removing linearly dependent modes with a Cosine similarity equal to or higher than c th = 0.993 lowers the condition number between two to four orders of magnitude, whereas the "shortened" Gram-Schmidt process (with excluded zero vectors) between two and almost nine orders of magnitude. A subsequent scaling of the columns of Φ f , i.e., converting  Table 1 Comparison of the effect of the Cosine similarity (c th = 0.993) and the Gram-Schmidt preconditioning on the condition number cond(Φ f ) of the flexible reduction matrix generated from the eigenmode selection, Fig. 6, of the beam-like models, Fig. 9. The table shows the number of excluded modes n x , the condition number of the reduced reduction matrix cond x (Φ f ) and of the reduced reduction matrix with scaled columns cond xs (Φ f ), for both methods. The modes are scaled after the Cosine similarity or Gram-Schmidt preconditioning

Models
cond(Φ f ) Cosine similarity precond. Gram-Schmidt precond. all flexible modes to the same Euclidean norm (length), lowers cond(Φ f ) of all beam-like models down to the order of 10 1 to 10 2 after the Cosine similarity preprocessing, or down to virtually one after the "shortened" Gram-Schmidt process. Scaling is a crucial preprocessing step for the formulation, as can be seen especially for the circular-sectioned beam, where scaling lowers the condition number by seven orders of magnitude! The exclusion of linearly dependent modes with | cos θ ij | ≥ 0.993, reduces the initially 72 generalized flexible coordinates by 33% in the best case and by at least 25%, even though only eight eigenmodes of vibration are included in the flexible reduction basis here. Hence, the gain in efficiency may be essential especially for larger problems. This shows that a neat preprocessing is not only important for accuracy, but may also for efficiency, due to the significant reduction of the DOFs. However, as already mentioned, the Gram-Schmidt preprocessing yields in general to a smaller reduction in the number of flexible DOFs, which is also reflected here, i.e., only between 4% and 16%, see Table 1.
The FE meshes for the beam-like models were chosen in order to preserve the lines of symmetry of the cross-sections, Fig. 9, except for the circle, where the infinite number of symmetry lines is reduced to four, due to the discretization. The arbitrary-shaped beam's cross-section, Fig. 9(e), was generated such that the geometry itself and the FE mesh do not exhibit any symmetries. Nevertheless, the flexible reduction matrices of the circular, the square and the arbitrary-shaped beam show the same linear dependencies, except that the correlations between the first torsion mode and the first bending mode pair, see Eq. (87), are slightly stronger for the square beam. Also, the Cosine similarity matrix of the arbitraryshaped beam show additionally some very weak correlations, yet with negligible values. Therefore, it seems that "small" deviations from symmetric geometries and the design of FE Fig. 10 Singular values σ i of the unmodified reduction matrix Φ of the circular beam-like model, Fig. 9(c), in descending order meshes have a minor influence on the dependencies of the flexible reduction matrices, and that not the shape but the number of symmetry lines dictates the dependencies. Finally, the rectangular and equilateral triangular beam exhibit the same correlations as the square beam, except for the dependencies within each generalized mode set generated from the bending modes, see Eq. (79).

Dependencies between the flexible and the rigid body motion part 6.1 General results
It is possible that the condition number of the reduction matrix cond(Φ) is still unacceptably high after preconditioning its flexible part. This happens whenever linear dependencies between the flexible Φ f and rigid body motion part [Φ t Φ r ] of the reduction basis arise. Such dependencies usually appear as linear combinations, i.e., if Eq. (60) is fulfilled with more than two vectors involved. The Cosine similarity fails to identify such linear combinations -also we cannot apply the "shortened" Gram-Schmidt process to the full reduction matrix, since the generalized translational and rotational component modes must remain unchanged (Φ t may be scaled). Which is why we have to resort to the SVD, briefly introduced in Sect. 3.2, to determine the nullspace of Φ. The nullspace identifies the generalized component modes involved in linear combinations. Therefore, we can iteratively (i) determine null(Φ), (ii) remove one of the generalized flexible component modes involved, and (iii) start with (i) again, as long as the condition number is as low as desired. 8 As mentioned in Sect. 3.2, it is in general required to distinguish between vanishing 9 and non-vanishing singular values σ i to determine the nullspace of a matrix. However, there is a sharp drop (several orders of magnitude) visible in the plot of the singular values of Φ, see Fig. 10, if linear dependencies exist. Figure 10 shows the singular values σ i of the unmodified reduction matrix Φ of the circular beam-like model, Fig. 9(c), in descending order; the especially distinctive singular value drop by almost six orders of magnitude occurs for the depicted example between the 72nd and 73rd singular value. The vanishing singular values after the drop indicate that Eq. (60) is fulfilled for nonzero coefficients a pi , i.e., the components of the right singular vectors v p of Φ, see Sect. 3. If Φ f is already well-conditioned (Cosine similarity or Gram-Schmidt preconditioning), linear dependencies between the rigid body motion and flexible part of the reduction matrix exist. These linear dependencies between Φ f and [Φ t Φ r ] manifest themselves as a sharp drop between two subsequent singular values, which is why no threshold for vanishing singular values needs to be determined, i.e., right singular vectors corresponding to singular values after the drop (here, 73 to 84) compose the nullspace of Φ, see Eq. (66). The corresponding right singular vectors indicate the columns of Φ involved in linear dependencies.
If linear dependencies arise, always a full set of three -for 3D problems -generalized flexible component modes m φ fkl with l = 1, 2, 3, see Eq. (20), must be removed, since the number of generalized component modes must be a multiple of three, so that the introduction of A bd is feasible, see Sect. 2.2. This is because the size of A bd is also a multiple of three for 3D problems, because the 3 × 3 rotation matrix is placed on its diagonal. The matrix multiplications involved in the term in the equation of motion, Eq. (58), make only sense if the inner dummy indices between, e.g., Φ and A T bd match, which is only true if the number of generalized flexible coordinates/modes is a multiple of three. It should be noted that this requirement is no restriction or problem, but will be satisfied a priori, because linear dependencies between, e.g., modes of the set { m φ fkl } with, say, l = 1 can only arise between another mode set with l = 1, since modes from a set with, say, l = 1 are always orthogonal to any other set where l = 2, 3, and if { m φ fkl } with, say, l = 1 is involved in a linear dependency, so is the corresponding mode set with l = 2, 3. This may be seen from the hypothetical one node and one mode example depicted in Eq. (28), which is restated with additional information below for a better understanding: In other words, the only possible dependencies for this example are but no "cross-dependencies" between modes with unequal l are possible, because of the null entries, see Eq. (91).

Further results for the beam-like models
For all analyzed beam-like models (Fig. 9), cond(Φ) is one order of magnitude larger than the condition number of the well-conditioned flexible reduction matrix. As already mentioned, it is not possible to scale the rotational part Φ r of the reduction matrix, however, it is important to ensure that the generalized component modes have similar norms. This is why the generalized translational and flexible component modes are scaled to the mean value of the norms of the columns of Φ r , which lowers the overall condition number of the already preconditioned total reduction matrix down to the order of 10 2 , and no dependencies between the rigid body motion and flexible part emerged for these examples.

Analysis of a crankshaft -a relevant engineering example
The previously analyzed beam-like examples (Sects. 4, 5 and 6.2) appear less often in practical applications, which is why this section is devoted to the analysis of a (steel) crankshaft of a reciprocal two-cylinder combustion engine. The crankshaft was discretized with 31 999 ten-noded tetrahedrals, yielding 50 289 nodes in total. The first 12 eigenmodes of vibration, covering the frequency range up to 10 000 Hz, are included in the flexible reduction basis and are visualized in Fig. 11. The Cosine similarity matrix of Φ f obtained from the crankshafts' eigenmodes of vibration is depicted in Fig. 12. The only dark-red entries are present in the upper corner of the heatmap with Cosine similarity values of approximately 0.97, which is, as already  Fig. 11 of a crankshaft. The first nine columns and rows correspond to the eigenmode depicted in Fig. 11(a), the following nine to the eigenmode of Fig. 11(b), and so forth (Color figure online)

Fig. 13
Singular values σ i of the unmodified reduction matrix Φ of the analyzed crankshaft, Fig. 11, in descending order covered in Sect. 5.1, considered too low to exclude generalized modes involved in these correlations. This is also confirmed by the "shortened" Gram-Schmidt process, since no zero vectors are generated. Hence, there are no linear dependencies within the flexible reduction matrix, which is relatively well-conditioned with cond(Φ f ) in the order of 10 2 . It may be orthonormalized by the standard Gram-Schmidt process to obtain a condition number of one.
In addition, there is no significant drop between subsequent singular values σ i of Φ, as shown in Fig. 13. Numerical experiments showed that the drop after σ 9 by approximately one order of magnitude is not caused by linear combinations between the flexible and rigid body motion part of Φ, i.e., no improvement in cond(Φ) is reached if the right singular vectors v i with i > 9 are considered to be in null(Φ) and the "pretended" linear combinations are removed via column exclusion. The condition number of the total reduction matrix Φ of the crankshaft is in the order of 10 4 , in contrast to the condition number of approximately 1.40 of the eigenmodes of vibration matrix generated from the modes depicted in Fig. 11. However, scaling of the generalized translational and flexible component modes to the mean value of the norms of the generalized rotational component modes, lowers the overall condition number of the reduction matrix Φ of the crankshaft down to the order of 10 2 and if the Gram-Schmidt process is applied to Φ f prior scaling down to approximately 41. The results are also summarized in Table 2. Table 2 Comparison of the condition numbers of the reduction matrices generated from the eigenmode selection of the crankshaft, Fig. 11. The table shows the condition of the flexible, cond(Φ f ), and full, cond(Φ), reduction matrices, the number of excluded modes n x , the condition of the scaled flexible reduction matrix cond s (Φ f ), of the full reduction matrix with scaled translational and flexible part (mean value of the norms of the rotational part), of the Gram-Schmidt orthonormalized flexible reduction matrix cond gs (Φ f ) and of the full reduction matrix with scaled and orthonormalized flexible and scaled translational part 8.75 × 10 2 6.22 × 10 4 0 6 .64 × 10 2 4.72 × 10 2 1.00 4.09 × 10 1 8 Step-by-step strategies to handle ill-conditioned reduction matrices The algorithm shown in the present section state strategies how to handle linear dependencies within the reduction matrix Φ and how to improve the overall condition number to convert ill-conditioned into well-conditioned problems.

Algorithm 1:
Pseudocode to obtain well-conditioned system matrices. For Cosine similarity preconditioning omit red instructions and for Gram-Schmidt preconditioning the blue ones (color version online) The pseudocode 10 contains the two different strategies discussed in detail in the previous sections, where one employs the Cosine similarity preconditioning (omit red instructions) and the other the "shortened" version of the Gram-Schmidt process (omit blue instructions). As already mentioned, the latter one should be preferred for practical applications, whereas the former allows more freedom during the mode selection process and provides more insight into the nature of the dependencies.

Conclusions
The present paper makes a contribution to the understanding and applicability of a promising 3D flexible multibody dynamics formulation known as generalized component mode synthesis. The novel nodal-based derivations of the governing equations of motion and the reduction basis in the already spatially discretized domain, in contrast to the standard continuum-mechanics-based derivation available in the literature, are presented in a traceable and intuitive way.
The main part of this contribution identified and resolved the weakness of the formulation, i.e., the possibility of ill-conditioned reduction matrices, which may introduce large errors or even lead to unsolvable problems -condition numbers up to nearly 10 18 -if not handled appropriately. The high condition numbers are caused by linear dependencies and Euclidean norms that deviate by orders of magnitude from each other of the generalized component modes -the reduction matrix's columns -generated from original linearly independent FE eigenmodes of vibration, undeformed FE nodal coordinates and the set of Cartesian base vectors.
It was shown by numerical experiments of simple extruded bodies with different crosssections that such linear dependencies may also arise for fully asymmetric geometries and FE meshes, which is why it is impossible to identify these dependencies a priori. It was shown how the singular value decomposition, the matrix condition number, the Cosine similarity and a "shortened" version of the classical Gram-Schmidt process may be used to identify and eliminate these linear dependencies, yielding not only sufficiently well-conditioned systems, but may also less DOFs and, therefore, a further gain in efficiency. The crankshaft example suggests that complex real-world engineering systems are potentially less likely to result in ill-conditioned reduction matrices, however, a neat preprocessing is suggested in any case to avoid unexpected errors.
Finally, step-by-step preprocessing procedures were derived to convert ill-into wellconditioned problems, which shall improve the formulation's applicability. Fig. 14 Illustration of orthonormal q * r versus non-orthonormal q r rotational DOFs. The generalized rotational coordinates from (b) were extracted from an elementary rotation matrix, i.e., q r = [cos(π/8) − 1 sin(π/8) − sin(π/8) cos(π/8) − 1] T , which is why the triangle is simply rotated by π/8 rad counter-clockwise. The generalized rotational coordinates from (c) were arbitrary chosen, so that Q r + I is not orthogonal, i.e., q r = 0.1[1 12 4 2] T , and the orthogonal contribution from (c) shown in (d) is obtained by Gram-Schmidt orthonormalization. Node 1 is the center of rotation, since x (1)  form (x, y). Hence, according to Eq. (3) and Eqs. (9) to (13),  (12), for different sets of rotational DOFs. The figure clearly indicates that the generalized rotational component modes contained in Φ r contain in general rigid body rotation, as well as shear and stretch deformations, and represent exclusively rigid body rotation only if Q r + I is orthogonal, Eq. (36). Therefore, part of the displacement c r = Φ r q r represents flexible deformation, i.e., Φ r (q r − q * r ), and contributes to the strain energy, Eq. (44). Figures 14(c)-(e) also show that a proper orthonormalization procedure separates the flexible part (e) and the rigid body motion part (d) contained in Φ r q r . In summary, Φ r q r → rotation and deformation.