Journal of Applied MathematicsVolume 2012 (2012), Article ID 638546, 12 pagesdoi:10.1155/2012/638546

Research Article

A Direct Eigenanalysis of Multibody System in Equilibrium

Cheng Yang,¹ Dazhi Cao,² Zhihua Zhao,² Zhengru Zhang,¹ and Gexue Ren²

¹School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China²School of Aerospace, Tsinghua University, Beijing 100084, China

Received 4 September 2011; Accepted 2 November 2011

Academic Editor: Massimiliano Ferronato

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

Abstract.

This paper presents a direct eigenanalysis procedure for multibody system in equilibrium. The first kind Lagrange’s equation of the dynamics of multibody system is a set of differential algebraic equations, and the equations can be used to solve the equilibrium of the system. The vibration of the system about the equilibrium can be described by the linearization of the governing equation with the generalized coordinates and the multipliers as the perturbed variables. But the multiplier variables and the generalize coordinates are not in the same dimension. As a result, the system matrices in the perturbed vibration equations are badly conditioned, and a direct application of the mature eigensolvers does not guarantee a correct solution to the corresponding eigenvalue problem. This paper discusses the condition number of the problem and proposes a method for preconditioning the system matrices, then the corresponding eigenvalue problem of the multibody system about equilibrium can be smoothly solved with standard eigensolver such as ARPACK. In addition, a necessary frequency shift technology is also presented in the paper. The importance of matrix conditioning and the effectiveness of the presented method for preconditioning are demonstrated with numerical examples.

1. Introduction

Modal analysis of multibody system [1, 2] has important application in structure analysis and modal reduction method [3–6]. The eigenanalysis methods for constrained multibody system can be mainly divided into two categories, namely, the eigenanalysis method by transforming the equations of motion from DAEs to ODEs by selecting independent coordinates [7–9] and the direct eigenanalysis method [10–12]. Compared with structural systems, the governing equation of constrained multibody system is differential algebraic equations, instead of the ordinary differential systems. In the reduction methods, the differential algebraic equations of multibody system are reduced into ordinary differential equations for independent general coordinates, through looking for a set of orthogonal basis of constraint Jacobian matrix to eliminate the constraints from the system governing equations. The reduction method is accurate, efficient, and stable provided that the degree of freedom of the system and the number of constraints are small. However, the reduction method is considered not suitable for large scale problems due to the spoilage of sparse matrices [13].

In the direct methods, the nonlinear differential and algebraic equations are simultaneously linearized about a given state of the system, that is, the algebraic constraint equations are kept during the eigenanalysis and the sparsity of system matrices is preserved. Depending on how the problem is actually formulated, the enlarged mass matrix may become structurally singular (eigenvalues associated to constraints become infinite); other approaches may result in a structurally singular enlarged stiffness matrix (eigenvalues associated to constraints become 0). Usually a frequency shift is applied, and structurally nonsingular enlarged matrices are resulted. But in most cases, matrices’ condition number of the frequency shifted system is very large and numerical singular due to the constraints and a corresponding precondition is needed. The root cause for numerical singularity can be understood by the fact that the dynamic equations and the Lagrange multipliers carry a dimension of force while the constraint equations and the general coordinate a dimension of length or angular. The literature [14, 15] suggested multiplying a parameter to the constraints’ corresponding rows and columns of system matrices as a preconditioning. The time step integration takes a parameter associate with time step and integration format to do the precondition. But the eigenanalysis about a static equilibrium has nothing to do with time integration, and we need to find an appropriate parameter.

For the direct eigenanalysis of constrained multibody system in static equilibrium, this paper proposes a procedure and a necessary preconditioning method of system matrices to guarantee correct eigensolution. This paper is organized as follows. Section 2 schematically describes the formulation of constrained multibody system based on Lagrange equation with multipliers and the direct eigenvalue problem about equilibrium. Section 3 presented a detailed condition number analysis of the system matrices and a preconditioning scaling technique. Section 4 introduces the selection of eigensolver and numerical examples.

2. Governing Equations of Multibody System and the Direct Eigenvalue Problem

The dynamic equations of a general multibody system are usually established with Lagrange’s equations with multipliers [16]. The system equations are where is the vector of general coordinates, is the vector of Lagrange multipliers, is the Lagrange function associated with , is the vector of general nonconservative forces which are functions of , is the vector of constraint functions associated with , and is the vector of generalized constraint forces. Denote Then (2.1) becomes The equilibrium state satisfies Assume small vibration about the equilibrium , then Carry out the Taylor expansion and reserve linear terms, we have where , , , , , , .

Write the solution as , where is one of the eigenvalues, is the corresponding eigenvector, and the eigenvalue problem is A frequency shift is applied to (2.7) by let , , and then where , , . The frequency shift makes us calculate the eigenvalues about and makes sure ’s reversibility (the probability of singular is zero when ) as the system’s energy must be in the form of kinetic energy or deformation energy. The eigenvalue problem can be transformed into the state space form as the following: where , , .

3. Ill-Condition in the Matrices of the Perturbed Vibration Equations and the Preconditioning

Consider the shifted eigenvalue problem (2.8) in the undamped case It is more convenient to consider the standard eigenvalue problem than the original generalized one As the inverse of matrix is used, cases in large condition number of could make the direct solution of result in erroneous eigensolution. As shown in Section 4, the eigenvalue problems (2.8) obtained by the linearization of the Lagrange equations with multipliers for practical multibody systems usually encounter serious condition number problem. The reason mainly lies in the different dimensions of the multiplier variables and constraint equations with that of generalized coordinates and governing equations.

Inspired by the eigenanalysis method through transforming the equations of motion from DAEs to ODEs by selecting independent coordinates, we divide the matrix into blocks and calculate its inverse by Gaussian elimination. Let where . Without loss of generality, let the first (the number of constraints) rows of be linearly independent and be denoted as . Divide the matrix as where , , , , , is the number of freedom of general coordinates.

Then where , , , . See appendix for detail of Gaussian elimination.

Assume Then it shows and by taking condition number

Equation (3.5) can be expected when the elements of eigenvector of the underlying ODE system have about the same magnitude, that is, the elastic parameters of the physical system are about the same magnitude in physical description, and this is often the case. The elastic parameters mean the Young’ modulus of flexible bodies, stiffness of spring forces, and so forth.

The condition number of the original system can be changed by scaling the multiplier variables and constraint equations. Let and multiply the constraint equations by in (2.1) we have The corresponding block matrices are Then the condition number analysis shows Choosing , the estimation of condition number decreases from about to about .

The elements of matrix stand for the elastic coefficients of the system and the elements of and work similarly as the stiffness of constraints. Taking as an enlarged stiffness matrix, the scaling is equivalent to changing the stiffness of constraints from to the same magnitude of physical stiffness.

Note that the condition cannot be improved by only scaling multipliers or by only scaling constraint equations. It is suggested that any different dimension in the procedure of modeling could lead to bad condition of numerical calculation.

The complete process for direct eigenanalysis of constrained deformable multibody system is thus as follows.

Step 1. Calculate the static equilibrium from (2.4).

Step 2. Calculate , , at .

Step 3. Precondition with which is the maximum absolute element of .

Step 4. Frequency shift with parameter with (2.8).

Step 5. Enlarge the matrix into state space as (2.9).

Step 6. Solve the eigenvalues and eigenvectors form (2.9) with eigensolver such as ARPACK.

Step 7. Add to the eigenvalues and multiply the eigenvector components of Lagrange multipliers with .

4. Selection of Eigensolver and Numerical Examples

The famous iteration method Implicitly Restarted Arnoldi Method (IRAM) [17, 18] is used in the numerical solution of eigenvalue problem (2.9) in the paper. Arnoldi method is one of the Krylov subspace methods [19] which calculate the largest few eigenvalues of a matrix. The presented numerical examples have verified the efficiency of IRAM in practical calculation. The famous parallel sparse linear system solver PARDISO [20, 21] is used to provide of (2.9) in the paper. A cantilever and a four-bar linkage are calculated to illustrate the effect of the condition number of the system matrix and the effectiveness of the presented preconditioning method. It is shown that the numerical results are physically meaningless without scaling and the precondition presented in the paper is very effective.

4.1. A Cantilever

In this example, beam element and 3-dimensional solid element are separately used to model the cantilever in the multibody frame, as shown in Figure 1. The parameters of the cantilever are shown in Table 1. The models are computed with the presented direct eigenanalysis method in two cases, that is, with and without axial tension, and the results are compared with analytical results based on the mathematical physics, see Tables 2 and 3.

Table 1: Parameters for the cantilever.

Item
Density
Young’s modulus
Poisson’s ratio
Beta damping
Size
Area
Inertia moment

Value
7860

0
0.0

Unit
kg/m³
N/m

m²
m⁴

Table 2: Eigenvalues comparison for cantilever ( Newton).

Frequency order
Unscaled model 1
Unscaled model 2
Euler-Bernoulli beams Analytical
Scaled model 1
Scaled model 2

value
error
value
error

1st-X
i
i
0.2067i
0.2065i

0.2208i

1st-Y
i
i
0.2067i
0.2065i

0.2208i

2nd-X
i
i
1.2957i
1.2906i

1.3701i

2nd-Y
i
i
1.2957i
1.2906i

1.3701i

3rd-X
i
i
3.6280i
3.6049i

3.7803i

3rd-Y
i
i
3. 6280i
3.6049i

3.7803i

Condition number

Note: error = (Algorithm−Analytical)/Analytical.

Table 3: Eigenvalues comparison for cantilever ( Newton).

Frequency order
Unscaled model 1
Unscaled model 2
Euler-Bernoulli beam Analytical
Scaled model 1
Scaled model 2

value
error
value
error

1st-X
i
i
0.2796i
0.2792i

0.2906i

1st-Y
i
i
0.2796i
0.2792i

0.2906i

2nd-X
i
i
1.3923i
1.3863i

1.4611i

2nd-Y
i
i
1.3923i
1.3863i

1.4611i

3rd-X
i
i
3.7127i
3.6876i

3.8602i

3rd-Y
i
i
3.7127i
3.6876i

3.8602i

Condition number

Note: error = (Algorithm−Analytical)/Analytical.

Figure 1: Two models of cantilever.

The governing equation and boundary conditions for the single-mode-free vibration of the cantilever under axial tension is as the following: where is the material coordinate, is the lateral vibration displacement mode, is the Young’s modulus, is the sectional moment of inertia, is the axial tension, is the linear density, and is the circular frequency of vibration. is the length of the cantilever and . With the conventional mathematical procedure, characteristic equation for frequency is where , , and .

It is shown in Tables 2 and 3 that the results without the scaling precondition are totally erroneous in both cases. In the two models, through computation for and for . Precondition parameter is chosen the same as and the condition numbers in the two computed cases are decreased about times after the precondition. The results show that the presented precondition procedure could guarantee correct results for the direct eigenanalysis, and the unscaled models give out meaningless results by ARPACK as the condition numbers are too large for the nonphysical contribution of constraints. Note that model 2 is used to illustrate the precondition’s effectiveness for models including other flexible elements such as 3D solid elements, and the error of scaled model 2 to Euler beam analytical results is partly because the solid model does not agree with the plan cross section assumption of Euler-Bernoulli beam.

4.2. Four-Bar Linkage

The four-bar linkage is shown in Figure 2, and the parameters of the bar are shown in Table 4. The model is modeled with Cartesian coordinates, and the constraints in the model are all revolute joints. The beam element model is used in this example, as shown in Figure 3. The equilibrium under gravity is first calculated with dynamic relaxation method, [22] and the eigenvalues are calculated about the deformed configuration. The equilibrium configuration of the four-bar system is also shown in Figure 3; the computed results with and without the preconditioning scaling are shown in Table 5 compared with results of FEA software ABAQUS with the same model. Again the unscaled results in the case are all meaningless as they have nonphysical real parts. In this problem and precondition parameter is chosen, the same value as . The condition numbers are reduced about times after the preconditioning, and correct results are obtained with the direct eigenanalysis method.

Table 4: Parameters for the linkage.

Item
Density
Young’s modulus
Poisson’s ratio
Beta damping
Area
Moment of inertia
Gravity

Value
7860

0
0

10

Unit
kg/m³
N/m

m²
m⁴
m/s²

Table 5: Eigenvalues with or without scaling.

Frequency order
ABAQUS
Unscaled model
Scaled model

value
error
value
error

1st
0.0805i
i
—
0.0808i

2nd
0.4215i
i
—
0.4236i

3rd
1.1922i
i
—
1.2134i

4th
1.6457i
i
—
1.6520i

5th
2.2631i
i
—
2.2633i

6th
2.6006i
i
—
2.5967i

Condition number

Note: error = (Algorithm−ABAQUS)/ABAQUS.

Figure 2: Four-bar linkage discrete model.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Modes of the four-bar linkage.

5. Summary and Conclusions

In this paper, the bad condition of the eigenvalue problem formulated from the direct linearization of the differential-algebraic equations of multibody system is reviewed and then discussed. The root cause may be attributed to the different dimensions of constraint multipliers and equations with physical system. Through analysis, a simple scale of the constraint multiplier variables and equations according to the stiffness of the elastic elements in the system can reduce the condition number of the system matrix to be inverted. After the preconditioning scale, eigensolvers such as ARPACK can be directly applied to the preconditioned system and yield correct result. The soundness of proposed approach and applicability are illustrated by comparing the results of a cantilever and a four-bar linkage with those obtained from either analysis or ABAQUS.

Appendix

A. Gaussian elimination to get (3.4) from (3.3)

Acknowledgment

Cheng Yang and Zhengru Zhang were supported by National NSF of China under Grant 11071124.

References

A. A. Shabana, “Flexible multibody dynamics: review of past and recent developments,” Multibody System Dynamics, vol. 1, no. 2, pp. 189–222, 1997.

A. A. Shabana, Dynamics of Multibody Systems, Cambridge University Press, Cambridge, UK, 2005.

W. H. A. Schilders, H. A. Van Der Vorst, and J. Rommes, Model Order Reduction: Theory, Research Aspects and Applications, vol. 13 of Mathematics in Industry, Springer-Verlag, Berlin, Germany, 2008.

M. Ripepi and P. Masarati, “Reduced order models using generalized eigenanalysis,” Proceedings of the Institution of Mechanical Engineers, Part K, vol. 225, no. 1, pp. 52–65, 2011.

K. Holm-Joslashrgensen and S. R.K. Nielsen, “A component mode synthesis algorithm for multibody dynamics of wind turbines,” Journal of Sound and Vibration, vol. 326, no. 3–5, pp. 753–767, 2009.

P. Koutsovasilis and M. Beitelschmidt, “Comparison of model reduction techniques for large mechanical systems: a study on an elastic rod,” Multibody System Dynamics, vol. 20, no. 2, pp. 111–128, 2008.

Y. Zhang, B. Wen, and S. Chen, “Eigenvalue problem of constrained flexible multibody systems,” Mechanics Research Communications, vol. 24, no. 1, pp. 11–16, 1997.

M. Da Lio and R. Lot, “Analisi modale di sistemi multibody descritti in coordinate naturali,” pp. 6-C9, 1999.

V. Cossalter, R. Lot, and F. Maggio, “The modal analysis of a motorcycle in straight running and on a curve,” Meccanica, vol. 39, no. 1, pp. 1–16, 2004.

P. Masarati, “Direct eigenanalysis of constrained system dynamics,” Proceedings of the Institution of Mechanical Engineers, Part K, vol. 223, no. 4, pp. 335–342, 2010.

V. Boschi, R. De Salvo, P. Masarati, G. Quaranta, and V. Sannibale, “Seismic attenuation system synthesis by reduced order models from multibody analysis,” in Proceedings of the Multibody Dynamics, Milano, Italy, June 2007.

J. L. Escalona and R. Chamorro, “Stability analysis of vehicles on circular motions using multibody dynamics,” Nonlinear Dynamics, vol. 53, no. 3, pp. 237–250, 2008.

M. Géradin and A. Cardona, Flexible Multibody Dynamics: A Finite Element Approach, Wiley & Sons, New York, NY, USA, 2001.

O. A. Bauchau, A. Epple, and C. L. Bottasso, “Scaling of constraints and augmented lagrangian formulations in multibody dynamics simulations,” Journal of Computational and Nonlinear Dynamics, vol. 4, no. 2, pp. 1–9, 2009.

C. L. Bottasso, D. Dopico, and L. Trainelli, “On the optimal scaling of index three DAEs in multibody dynamics,” Multibody System Dynamics, vol. 19, no. 1-2, pp. 3–20, 2008.

C. Lanczos, The Variational Principles of Mechanics, Dover Publications, New York, NY, USA, 1970.

R. B. Lehoucq and D. C. Sorensen, “Deflation techniques for an implicitly restarted Arnoldi iteration,” SIAM Journal on Matrix Analysis and Applications, vol. 17, no. 4, pp. 789–821, 1996.

H. Saberi Najafi and E. Khaleghi, “A new restarting method in the Arnoldi algorithm for computing the eigenvalues of a nonsymmetric matrix,” Applied Mathematics and Computation, vol. 156, no. 1, pp. 59–71, 2004.

Y. Saad, Numerical Methods for Large Eigenvalue Problems, Manchester University Press, Manchester, UK, 1992.

O. Schenk and K. Gartner, “Solving unsymmetric sparse systems of linear equations with PARDISO,” in Computational Science-ICCS 2002, Part II, P. Sloot, C. J. K. Tan, J. J. Dongarra, and A. G. Hoekstra, Eds., vol. 2330, pp. 355–363, Springer, Berlin, Germany, 2002.

O. Schenk and K. Gartner, “On fast factorization pivoting methods for sparse symmetric indefinite systems,” Electronic Transactions on Numerical Analysis, vol. 23, pp. 158–179, 2006.

W. J. Lewis, Tension Structures: Form and Behaviour, Thomas Telford Services, 2003.