# A Direct Eigenanalysis of Multibody System in Equilibrium

Authors: Cheng Yang; Dazhi Cao; Zhihua Zhao; Zhengru Zhang; Gexue Ren
Journal: Journal of Applied Mathematics
 (2012)

---

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


---

## Body

## 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. Section2 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
(2.1)ddt(∂L∂ẋ)-∂L∂x=(∂Φ∂x)Tλ+f,Φ(x)=0,
where x∈ℝn is the vector of general coordinates, λ∈ℝm is the vector of Lagrange multipliers, L is the Lagrange function associated with (ẋ,x), f is the vector of general nonconservative forces which are functions of (ẋ,x,t), Φ is the vector of constraint functions associated with (x), and (∂Φ/∂x)Tλ is the vector of generalized constraint forces. Denote(2.2)F(ẍ,ẋ,x,t)=(-ddt(∂L∂ẋ)+∂L∂x+f)(ẍ,ẋ,x,t),G(x)=(∂Φ(x)∂x)T.
Then (2.1) becomes
(2.3)F(ẍ,ẋ,x,t)+G(x)λ=0,Φ(x)=0.
The equilibrium state (x0,λ0,t0) satisfies(2.4)F(0,0,x0,t)+G(x0)λ0=0,Φ(x0)=0.
Assume small vibration about the equilibrium (δẍ,δẋ,x0+δx,λ0+δλ,t), then
(2.5)F(δẍ,δẋ,x0+δx,t)+G(x0+δx)(λ0+δλ)=0,Φ(x0+δx)=0.
Carry out the Taylor expansion and reserve linear terms, we have
(2.6)M̂(x0)δÿ+Ĉ(x0)δẏ+K̂(x0,λ0)δy=0,
where M̂(x0)=[M(x0)000], Ĉ(x0)[C(x0)000], K̂(x0,λ0)=[K(x0,λ0)G(x0)GT(x0)0], M(x0)=(∂F/∂ẍ)(x0), C(x0)=(∂F/∂ẋ)(x0), K(x0,λ0)=(∂(F+Gλ)/∂x)(x0,λ0), δy=[δxδλ].Write the solution asδy(t)=ertV, where r is one of the eigenvalues, V is the corresponding eigenvector, and the eigenvalue problem is
(2.7)(r2M̂+rĈ+K̂)V=0.
A frequency shift is applied to (2.7) by let r̅=r-α, α∈ℝ+, and then
(2.8)(r̅2M̅+r̅C̅+K̅)V=0,
where M̅=M̂, K̅=α2M̂+αĈ+K̂, C̅=2αM̂+Ĉ. The frequency shift makes us calculate the eigenvalues about α and makes sure K̅’s reversibility (the probability of singular is zero when α≠0) 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:
(2.9)r̅M̃Y=K̃Y,
where Y=(yr̅y), M̃=[-C̅-M̅I0], K̃=[K̅00I].


## 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
(3.1)[r̅2M̅+K̅]V=0.
It is more convenient to consider the standard eigenvalue problem than the original generalized one
(3.2)[K̅-1M̅+r̅-2I]V=0.
As the inverse of matrix K̅ is used, cases in large condition number of K̅ could make the direct solution of K̅-1M̅ 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α2M+K=cD where c=∥α2M+K∥∞. Without loss of generality, let the first m (the number of constraints) rows of G be linearly independent and be denoted as G1∈ℝm×m. Divide the matrix as
(3.3)K̅=[cD11cD12G1cD21cD22G2G1TG2T0],
where D11∈ℝm×m, D22∈ℝ(n-m)×(n-m), D12∈ℝm×(n-m), D21∈ℝ(n-m)×m, G2∈ℝ(n-m)×m, n is the number of freedom of general coordinates.Then(3.4)K̅-1=c-1[G1-TG2TP-1G2G1-1-G1-TG2TP-1c(G1-T+G1-TG2TP-1R)-P-1G2G1-1P-1-cP-1Rc(G1-1+QP-1G2G1-1)-cQP-1c2(QP-1R-N)],
where N=G1-1D11G1-T, P=D22+G2NG2T-D21G1-TG2T-G2G1-1D12, Q=G1-1D12-NG2T, R=D21G1-T-G2N. See appendix for detail of Gaussian elimination.Assume(3.5)‖N‖∞≈‖P-1‖∞≈‖Q‖∞≈‖R‖∞=O(1).
Then it shows ∥K̅-1∥∞=O(c) and by taking condition number cond(·)=∥·∥∞∥·-1∥∞(3.6)cond(K̅)≈O(c2).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λ̂=k-1λ and multiply the constraint equations by k in (2.1) we have
(3.7)ddt(∂L∂ẋ)-∂L∂x=k(∂Φ∂x)Tλ̂+f,kΦ(x)=0.
The corresponding block matrices are(3.8)K̅=[cD11cD12kG1cD21cD22kG2kG1TkG2T0],K̅-1=c-1[G1-TG2TP-1G2G1-1-G1-TG2TP-1k-1c(G1-T+G1-TG2TP-1R)-P-1G2G1-1P-1-k-1cP-1Rk-1c(G1-1+QP-1G2G1-1)-k-1cQP-1k-2c2(QP-1R-N)].
Then the condition number analysis shows
(3.9)cond(K̅)≈O(max⁡⁡(k,c)×max⁡⁡(k-1c,k-2c2)).
Choosing k≈c, the estimation of condition number decreases from about O(c2) to about O(1).The elements of matrixK stand for the elastic coefficients of the system and the elements of G and GTwork similarly as thestiffness of constraints. Taking K̅ as an enlarged stiffness matrix, the scaling is equivalent to changing the stiffness of constraints from O(1) 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(x0,λ0) from (2.4).Step 2.
CalculateM̂, Ĉ, K̂ at (x0,λ0).Step 3.
PreconditionK̂ with c which is the maximum absolute element of K̂.Step 4.
Frequency shift with parameterα(1~10) 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 c.


## 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 K̃-1 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 Figure1. 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.
ItemDensityYoung’s modulusPoisson’s ratioBeta dampingSizeAreaInertia momentValue78602.06e1100.010*10*200102(10)4/12Unitkg/m3N/m——m*m*mm2m4Table 2
Eigenvalues comparison for cantilever (F=0 Newton).
Frequency  orderUnscaled model 1Unscaled model 2Euler-Bernoulli beams AnalyticalScaled model 1Scaled model 2valueerrorvalueerror1st-X0.15+0.43i-0.01+0.26i0.2067i0.2065i-0.10%0.2208i6.82%1st-Y-0.47+0.49i-0.49+0.33i0.2067i0.2065i-0.10%0.2208i6.82%2nd-X-0.76+0.03i-0.38+0.95i1.2957i1.2906i-0.40%1.3701i5.74%2nd-Y-0.87+0.01i-0.73+1.35i1.2957i1.2906i-0.40%1.3701i5.74%3rd-X-1.45+1.17i-0.41+1.49i3.6280i3.6049i-0.64%3.7803i4.20%3rd-Y-1.44+0.38i-2.68+0.11i3. 6280i3.6049i-0.64%3.7803i4.20%Condition number3.9409e+222.3198e+199.3813e+063.4183e+07Note: error = (Algorithm−Analytical)/Analytical.Table 3
Eigenvalues comparison for cantilever (F=1e10 Newton).
Frequency  orderUnscaled model 1Unscaled model 2Euler-Bernoulli beam AnalyticalScaled model 1Scaled model 2valueerrorvalueerror1st-X0.12+0.42i-0.03+0.28i0.2796i0.2792i-0.14%0.2906i3.93%1st-Y-0.87+0.19i-0.16+0.60i0.2796i0.2792i-0.14%0.2906i3.93%2nd-X-1.14+0.05i0.07+1.35i1.3923i1.3863i-0.43%1.4611i4.94%2nd-Y1.27+0.22i1.45+0.88i1.3923i1.3863i-0.43%1.4611i4.94%3rd-X-0.18+1.38i-1.06+1.41i3.7127i3.6876i-0.68%3.8602i3.97%3rd-Y-1.26+0.71i-0.23+2.73i3.7127i3.6876i-0.68%3.8602i3.97%Condition number4.3103e+217.9065e+185.1001e+066.4188e+07Note: 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:(4.1)EId4ydx4-Fd2ydx2-ρω2y=0,y|0=0dydx|0=0d2ydx2|l=0d3ydx3|l=adydx|l,
where x is the material coordinate, y is the lateral vibration displacement mode, E is the Young’s modulus, I is the sectional moment of inertia, F is the axial tension, ρ is the linear density, and ω is the circular frequency of vibration. l is the length of the cantilever and a=F/EI. With the conventional mathematical procedure, characteristic equation for frequency ω is
(4.2)[r12cosh⁡(r1l)+r22cos⁡(r2l)]{r13sinh⁡(r1l)-r23sin⁡(r2l)+a[-r1sinh⁡(r1l)-r2sin⁡(r2l)]}=[r12sinh⁡(r1l)+r1r2sin⁡(r2l)]{r13cosh⁡(r1l)+r1r22cos⁡(r2l)+a[-r1cosh⁡(r1l)-r1cos⁡(r2l)]},
where b=ρω2/EI, r1=(a2+4b+a)/2, and r2=(a2+4b-a)/2.It is shown in Tables2 and 3 that the results without the scaling precondition are totally erroneous  in both cases. In the two models, through computation c=1.39167e+06 for F=0 and c=1.39252e+06 for F=1e+10. Precondition parameter k is chosen the same as c and the condition numbers in the two computed cases are decreased about O(c2) 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 Figure2, 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 c=9.11047e+13 and precondition parameter k is chosen, the same value as c. The condition numbers are reduced about O(c) times after the preconditioning, and correct results are obtained with the direct eigenanalysis method.Table 4
Parameters for the linkage.
ItemDensityYoung’s modulusPoisson’s ratioBeta dampingAreaMoment of inertiaGravityValue78602.06e1100102(10)4/1210Unitkg/m3N/m——m2m4m/s2Table 5
Eigenvalues with or without scaling.
Frequency orderABAQUSUnscaled modelScaled modelvalueerrorvalueerror1st0.0805i-1.09+0.02i—0.0808i0.37%2nd0.4215i-1.06+0.34i—0.4236i0.50%3rd1.1922i-1.20+0.35i—1.2134i1.78%4th1.6457i-1.30+0.08i—1.6520i0.38%5th2.2631i-1.30+0.03i—2.2633i0.01%6th2.6006i-1.34+0.23i—2.5967i-0.15%Condition number3.3416e+191.1560e+10Note: error = (Algorithm−ABAQUS)/ABAQUS.Figure 2
Four-bar linkage discrete model.Modes of the four-bar linkage.
(a)(b)(c)(d)(e)(f)


## 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 Figure1. 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.
ItemDensityYoung’s modulusPoisson’s ratioBeta dampingSizeAreaInertia momentValue78602.06e1100.010*10*200102(10)4/12Unitkg/m3N/m——m*m*mm2m4Table 2
Eigenvalues comparison for cantilever (F=0 Newton).
Frequency  orderUnscaled model 1Unscaled model 2Euler-Bernoulli beams AnalyticalScaled model 1Scaled model 2valueerrorvalueerror1st-X0.15+0.43i-0.01+0.26i0.2067i0.2065i-0.10%0.2208i6.82%1st-Y-0.47+0.49i-0.49+0.33i0.2067i0.2065i-0.10%0.2208i6.82%2nd-X-0.76+0.03i-0.38+0.95i1.2957i1.2906i-0.40%1.3701i5.74%2nd-Y-0.87+0.01i-0.73+1.35i1.2957i1.2906i-0.40%1.3701i5.74%3rd-X-1.45+1.17i-0.41+1.49i3.6280i3.6049i-0.64%3.7803i4.20%3rd-Y-1.44+0.38i-2.68+0.11i3. 6280i3.6049i-0.64%3.7803i4.20%Condition number3.9409e+222.3198e+199.3813e+063.4183e+07Note: error = (Algorithm−Analytical)/Analytical.Table 3
Eigenvalues comparison for cantilever (F=1e10 Newton).
Frequency  orderUnscaled model 1Unscaled model 2Euler-Bernoulli beam AnalyticalScaled model 1Scaled model 2valueerrorvalueerror1st-X0.12+0.42i-0.03+0.28i0.2796i0.2792i-0.14%0.2906i3.93%1st-Y-0.87+0.19i-0.16+0.60i0.2796i0.2792i-0.14%0.2906i3.93%2nd-X-1.14+0.05i0.07+1.35i1.3923i1.3863i-0.43%1.4611i4.94%2nd-Y1.27+0.22i1.45+0.88i1.3923i1.3863i-0.43%1.4611i4.94%3rd-X-0.18+1.38i-1.06+1.41i3.7127i3.6876i-0.68%3.8602i3.97%3rd-Y-1.26+0.71i-0.23+2.73i3.7127i3.6876i-0.68%3.8602i3.97%Condition number4.3103e+217.9065e+185.1001e+066.4188e+07Note: 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:(4.1)EId4ydx4-Fd2ydx2-ρω2y=0,y|0=0dydx|0=0d2ydx2|l=0d3ydx3|l=adydx|l,
where x is the material coordinate, y is the lateral vibration displacement mode, E is the Young’s modulus, I is the sectional moment of inertia, F is the axial tension, ρ is the linear density, and ω is the circular frequency of vibration. l is the length of the cantilever and a=F/EI. With the conventional mathematical procedure, characteristic equation for frequency ω is
(4.2)[r12cosh⁡(r1l)+r22cos⁡(r2l)]{r13sinh⁡(r1l)-r23sin⁡(r2l)+a[-r1sinh⁡(r1l)-r2sin⁡(r2l)]}=[r12sinh⁡(r1l)+r1r2sin⁡(r2l)]{r13cosh⁡(r1l)+r1r22cos⁡(r2l)+a[-r1cosh⁡(r1l)-r1cos⁡(r2l)]},
where b=ρω2/EI, r1=(a2+4b+a)/2, and r2=(a2+4b-a)/2.It is shown in Tables2 and 3 that the results without the scaling precondition are totally erroneous  in both cases. In the two models, through computation c=1.39167e+06 for F=0 and c=1.39252e+06 for F=1e+10. Precondition parameter k is chosen the same as c and the condition numbers in the two computed cases are decreased about O(c2) 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 Figure2, 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 c=9.11047e+13 and precondition parameter k is chosen, the same value as c. The condition numbers are reduced about O(c) times after the preconditioning, and correct results are obtained with the direct eigenanalysis method.Table 4
Parameters for the linkage.
ItemDensityYoung’s modulusPoisson’s ratioBeta dampingAreaMoment of inertiaGravityValue78602.06e1100102(10)4/1210Unitkg/m3N/m——m2m4m/s2Table 5
Eigenvalues with or without scaling.
Frequency orderABAQUSUnscaled modelScaled modelvalueerrorvalueerror1st0.0805i-1.09+0.02i—0.0808i0.37%2nd0.4215i-1.06+0.34i—0.4236i0.50%3rd1.1922i-1.20+0.35i—1.2134i1.78%4th1.6457i-1.30+0.08i—1.6520i0.38%5th2.2631i-1.30+0.03i—2.2633i0.01%6th2.6006i-1.34+0.23i—2.5967i-0.15%Condition number3.3416e+191.1560e+10Note: error = (Algorithm−ABAQUS)/ABAQUS.Figure 2
Four-bar linkage discrete model.Modes of the four-bar linkage.
(a)(b)(c)(d)(e)(f)


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



---

*Source: 638546-2012-02-16.xml*