A numerical algorithm for nonlinear dynamic problems based on BEM

Engineering Analysis with Boundary Elements 23 (1999) 503–513

A numerical algorithm for nonlinear dynamic problems based on BEM H.J. Holl a,*, A.K. Belyaev b, H. Irschik a a

Johannes-Kepler-University, Institute of Mechanics and Machine Design, Division of Technical Mechanics, Altenbergerstrasse 69, A-4040 Linz, Austria b Technical University of St. Petersburg, St. Petersburg, Russia

Abstract A semi-analytical time-integration procedure for the integration of discretized dynamic mechanical systems is presented. This method utilizes the advantages of the boundary element method (BEM), well known from quasi-static field problems. Motivated by these spatial formulations, the present dynamic method is based on influence functions in time, and gives exact solutions in the linear time-invariant case. Similar to domain-type BEM’s for nonlinear field problems, the method is extended for different nonlinear dynamic systems having nonclassical damping and time-varying mass. The numerical stability and accuracy of the semi-analytical method is discussed in two steps for the nonclassical damping and for the nonlinear restoring forces, e.g. of the Duffing type. The damped Duffing oscillator and a linear oscillator with time-varying mass are used as representative model problems. For a nonlinear rotordynamic system, a comparison is given to other conventionally used time integration procedures, which shows the efficiency of the present method. 䉷 1999 Elsevier Science Ltd. All rights reserved. Keywords: Time-integration; Nonlinear dynamic systems; Time-varying mass; Rotordynamic system

1. Introduction During the last decades, various powerful computational methods have been developed for nonlinear structural dynamics. A main stage of the computational analysis is the space-wise discretization of the governing field equations leading to a system of ordinary differential equations in time. Especially, a semi-discretization by the Finite Element Method (FEM) leads to a dynamic system of second-order equations of the form  ⫹ CX _ ⫹ KX ˆ F; MX

…1†

compare Zienkiewicz, Taylor [1], where X denotes the vector of generalised coordinates, M stands for the mass matrix, K and C are matrices of stiffness and damping, respectively, and F gathers the imposed forcing terms. In many practical situations nonlinearities do occur, the matrices in Eq. (1) then depending on the generalised coordinates   _ C ˆ c ⫹ C…X; _ KX ˆ kX ⫹ Q…X; X†; _ M ˆ m ⫹ M…X; X†; X†; …2†

where m, c and k denote the symmetric and constant * Corresponding author. Tel.: ⫹43 732 2468 9766; fax: ⫹43 732 2468 9763. E-mail address: [email protected] (H.J. Holl)

portions of the matrices in Eq. (1). The damping matrix c is of the Rayleigh type. It has been shown in Zienkiewicz, Taylor [1] that the numerical integration of Eq. (1) by a finite element discretization in time provides a strategy unifying many existing single- and multi-step algorithms and providing a variety of new ones. When discussing stability problems of these schemes, it is useful to revert to modally uncoupled equations, see [1] and Hughes [2]: _ X†;  X i ⫹ 2zi vi X_ i ⫹ v2i Xi ˆ pi ˆ fi ⫹ qi …X; X;

…3†

where the index (i) refers to the number of the mode, which is supressed in the following. The modal analysis transforming Eq. (1) into Eq. (3) is performed with respect to the constant symmetric matrices m, c and k of Eq. (2), while the remaining unsymmetric and nonlinear portions are gathered in the modal forces qi. The modal force arising from F is denoted by fi, v is the natural frequency, and z denotes slight modal damping. With respect to Eq. (3), it is noted that our group has been successful to solve various structural elasto-viscoplastic problems by utilizing an analogy between the inelastic strains and a distribution of eigenstrains in a linear-elastic background-structure and by taking the continuous modes of that linear structure for a modal reduction to equations of the type of Eq. (3). The nonlinear restoring forces qi were calculated from the constitutive equations in a time-stepping procedure, and a semi-analytic time-integration of Eq. (3) was performed, see

0955-7997/99/$ - see front matter 䉷 1999 Elsevier Science Ltd. All rights reserved. PII: S0955-799 7(98)00105-2

504

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Irschik, Ziegler [3,4], Fotiu, Irschik, Ziegler [5] for reviews, and Fotiu, Irschik, Ziegler [6] for a recent application relating modal analysis to integral equations for elastic-plastic plate vibrations. For the discretized case of nonlinear rotordynamics, a semi-analytic time-integration scheme analogous to the scheme in Refs. [3–6] has been developed by Holl [7,8], see also Holl, Irschik [9,10]. Seen from the point of the FEM, however, the time-integration of Eqs. (1) or (3) by means of finite elements in time of course is a concise strategy. One should remind, nevertheless, that a formulation of the type of Eq. (1) can also be obtained by utilising the Boundary Element Method (BEM) for the space-wise semi-discretization, Nardini and Brebbia [11], see Manolis and Beskos [12] for a review. In that connection, it is the scope of the present contribution to discuss the integration of Eq. (3) by means of a BEM in time. It is shown that a time-wise BEM corresponds to an exact integration of the modal equations in the linear case, while a domain-boundary element formulation (DBEM) is obtained for the case of nonlinear modal restoring forces, compare [13] for DBEMs of quasi-static inelastic problems. The time-domain integrals of the present formulation coincide with Duhamel’s convolution integral, well-known from the theory of structural vibrations, see e.g. Ziegler [14]. In contrast to an approximation with FEM in time, the present BEM does not approximate the generalised coordinates, but only the integrands of the domain integrals occurring in the nonlinear case are approximated. It is noted that the present formulation conceptually coincides with the above cited integration scheme of Refs. [3–6], which now is connected to the BEM. In the present article special emphasis is given to the question of numerical stability in the nonlinear case, where several approximations of the domain-integrals are compared. Results are shown for the Duffing oscillator and are compared to the outcomes of direct integration schemes. Further examples are given for systems with time-varying mass and compared to analytical solutions derived. In the present article additionally a rotordynamic system is studied in some detail. Short accounts of the present formulation have been presented in [15–17].

2. Boundary integral formulation in time Following the fundamental concepts of BEM presented by Brebbia, Dominguez [18], we start with a weighted residual formulation for Eq. (3) within the time interval: ZDt 0

 ⫹ 2zvX…t† _ ⫹ v2 X…t† ⫺ p…t†ŠW…t†dt ˆ 0; ‰X…t†

…4†

0 ⱕ t ⱕ Dt; W(t) denotes a sufficiently continuous weight function.

Integrating by parts twice gives ZDt 0

 ⫺ 2zvW…t† _ ⫹ v2 W…t††X…t† ⫺ p…t†W…t†Šdt ‰…W…t†

Dt Dt _ Dt _ ⫹ X…t†W…t†兩 0 ⫹ 2zvX…t†W…t†兩0 ⫺ X…t†W…t†兩0

ˆ 0:

…5†

In order to relate Eq. (5) to BEM, we follow the singular approach indicated by Brebbia, Dominguez [18]. In the sense of this latter approach, the weight function W is identified by the fundamental solution W1 of the modal oscillator owing to an unit impulse acting at the inner site of the boundary t ˆ 0. W1 is governed by the differential equation W 1 ⫺ 2zvW_ 1 ⫹ v2 W1 ˆ d…t⫹ †;

…6†

where W1 …t† ˆ

1 zvt e sinvd t H…t⫹ †; vd

q vd ˆ v 1 ⫺ z2 :

…7†

Dirac’s delta function is denoted by d , and t⫹ ˆ lim1!0 …t ⫺ 1†. W1 satisfies homogeneous initial conditions at t ˆ 0. As a second fundamental solution, the solution W2 of the modal oscillator as a result of a unit dipole d_ is introduced: W 2 ⫺ 2zvW_ 2 ⫹ v2 W2 ˆ d_ …t⫹ †;

…8a†

where W2 …t† ˆ ezvt ‰cosvd t ⫹

vz sinvd tŠH…t⫹ †: vd

…8b†

Applying both, Eqs. (7), (8a) and (8b), to Eq. (5), considering the properties of d and d_ , ZDt 0

ZDt 0

X…t†d…t†dt ˆ X…t ˆ 0† ˆ X0 ; …9† _ ˆ 0† ˆ ⫺X_ 0 ; X…t†d_ …t†dt ˆ ⫺X…t

using some trigonometric transformation rules and solving for the state-vector at the end of the time-interval, the following transition matrix formulation is obtained: ( ) ( ) ( ) ZDt 0 X…Dt† X0 A…Dt ⫺ t† p…t†dt; …10† ⫹ ˆ A…Dt† _ 0 X_ 0 1 X…Dt† where

2

vzsinvd t 6 cosvd t ⫹ vd 6 A…t† ˆ e⫺zvt 6 6 4 v2 sinvd t ⫺ vd q and vd ˆ v 1 ⫺ z2 :

sinvd t vd

3

7 7 7 7 vzsinvd t 5 cosvd t ⫺ vd

…11†

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

It is thus seen that the evaluation of the weighted residual statement of Eq. (5) in the sense of BEM leads to the exact integration of the modal Eq. (3) by means of Duhamel’s integral, as A is the exact transition matrix of the linear oscillator, compare Hughes [2], and see Refs. [3–6], where a semi-analytical time-stepping procedure analogous to Eqs. (10) and (11) has been used. It is certainly not surprising that the use of BEM in time leads to exact formulations, as the boundaries of the time-interval reduce to a two-point problem, with algebraic relations instead of boundary integrals. In the case of nonlinear problems, the portion q of p in Eq. (3) is not known in advance. An iterative procedure has to be established in the sense of the DBEM by approximating the time-evolution of q in Eq. (10), see the following section. The solution for the load transfer matrix of the modal equations of motion owing to the excitation forces can be computed using (

X…Dt†

ˆ A…Dt†

_ X…Dt†

b11

b22

X0 X_ 0

)

( ) ZDt 0 A…Dt ⫺ t† fkm …t†dt: ⫹ 0 1 …12†

parameters of the interpolation function in the form:

f km

8 9 fkm > > > > < = ˆ Dfkm : > > > > : ; DDfkm

…14†

The above Eqs. (12) and (13) are valid for the k-th modal coordinate. In the following index k, which refers to the number of the generalized coordinate, is omitted as it is obvious in these equations. The evaluation of the convolution integral of Eq. (12) for the above stated interpolation function (13) gives the corresponding load transfer matrix involving dynamic Green’s functions. After evaluation of Eq. (12) for all parts of the solution, we get the following formulas for the elements of the load transfer matrix: " B…Dt† ˆ

b11

b12

b13

b21

b22

b23

# :

…16a†

1 ⫺zvDt e sin…vd Dt†; vd

…16b†

" # vDt ⫺ 2z e⫺zvDt 1 ⫺ 2z 2 ˆ ⫹ 3 2zcos…vd Dt† ⫺ p sin…vd Dt† ; v Dt v3 Dt 1 ⫺ z2 " #) ( 1 zsin…vd Dt† ⫺zvDt ˆ 2 cos…vd Dt† ⫹ p 1⫺e ; v Dt 1 ⫺ z2

b13 ˆ 8

…15†

The used abbreviations can be expressed by the following equations:

" #) ( 1 zsin…vd Dt† ⫺zvDt ; ˆ 2 1⫺e cos…vd Dt† ⫹ p v 1 ⫺ z2

b21 ˆ

b12

(

)

505

…16c†

…16d†

2 2 1 ⫹ zvt ⫺ 4z2 4 8 ⫺zvDt …1 ⫺ 2z †vDt ⫹ 2z…3 ⫺ 4z † p   sin…vd Dt† ⫺ 4 2 e⫺zvDt …1 ⫺ zvDt ⫺ 4z2 †cos…vd Dt†; ⫺ e v4 Dt2 v4 Dt2 v Dt 1 ⫺ z2

…16e† b23 ˆ

4 4 2…1 ⫺ 2z2 † ⫺ zvDt p sin…vd Dt†: ‰…4z ⫺ vDt† ⫺ e⫺zvDt …4z ⫹ vDt†cos…vd Dt†Š ⫹ 3 2 e⫺zvDt 2 v Dt v Dt 1 ⫺ z2 3

The excitation force in the integral is approximated in the form of a parabolic interpolation function of the force within the m-th time step according to fkm …t† ˆ fkm0 ⫹ Dfkm

 2 t t ⫹ 4DDfkm : Dt Dt

…13†

We can define a corresponding forcing vector with the three

…16f†

For the above given BEM-type algorithm the unified transfer matrix method of Eq. (12) including the transfer matrix A and the load transfer matrix B is defined. These multiplications have to be performed for every time step. If we write the result for the generalized coordinate at the end of the m-th time step, we get a linear system m X Am⫺p B f~ p : …17† q~ m ˆ Am q~ 0 ⫹ pˆ1

506

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Fig. 1. (a) Absolute values of the eigenvalue of the transition matrix for the approximation of Eq. (20a); (b) Absolute values of the eigenvalue of the transition matrix for the approximation of Eq. (20b); (c) Absolute values of the eigenvalue of the transition matrix for the approximation of Eq. (20c).

This is the exact solution if the parabolic interpolation of the excitation force is exact. Again in Eq. (23) the index k of the generalized coordinate is omitted and the index m corresponds to the m-th time step. The application of modal reduction gives additional advantages for the computational effort, which does only have a neglectable influence on the accuracy. A study of the number of modes required to be included in the computational analysis is discussed in [19] and [20]. Modal reduction can be included with some small changes in the algorithm. As modal analysis is involved only when computing the dynamic parts of the solution, the implementation of modal reduction is a consequence when considering the computational effort. In this procedure the modal reduction is done in a consistent way as the results from the higher modes are very small and hence negligible.

3. Analysis of systems with nonlinear restoring forces The transition-matrix formulation of Eq. (12) provides exact solutions for linear vibrations. In order to obtain a measure for the quality of various approximations of the modal force qi in the nonlinear case, however, free damped

linear vibrations are studied at first, where the damping force formally is treated as a nonlinear restoring force: _ _ ˆ ⫺2zX: X ⫹ v2 X ˆ q…X†

…18†

Using the undamped matrix A…t†zˆ0 , the solution of Eq. (18) reads, compare Eq. (10): ( ) ( ) ZDt X0 X…Dt† A…Dt ⫺ t†zˆ0 ⫺ 2z ˆ A…Dt†zˆ0 _ 0 X_ 0 X…Dt† ( ) 0 _  X…t†dt: …19† 1 Three types of approximations of the integrand in Eq. (19) are studied: Xa …t† ˆ X0 ⫹

t _ 1 _ …X…Dt† ⫹ X_ 0 †; X_ a ˆ …X…Dt† ⫹ X_ 0 †: …20a† 2 2

Xb …t† ˆ X0 ⫹

t 1 …X…Dt† ⫺ X0 †; X_ b …t† ˆ …X…Dt† ⫺ X0 †; Dt Dt …20b†

2t X_ c …t† ˆ X…Dt† ⫹ 2 …X…Dt† ⫺ X0 ⫺ X_ 0 Dt†: Dt

…20c†

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Fig. 2. Absolute values of the eigenvalues for the system of Eq. (25).

Inserting Eq. (20) into Eq. (19), integrating and rearranging terms gives the form ( ) ( ) X…Dt† X0 ~ …21† ˆ A…Dt† _ X_ 0 X…Dt† appropriate for a stability discussion. For example the transition matrix of case (b) reads ~ b …Dt† ˆ A "

1 vDt ⫹ 2z…1 ⫺ cosvDt†

2z ⫹ …vDt ⫹ 2z†cosvDt

DtsinvDt

⫺v DtsinvDt

⫺2z ⫹ …vDt ⫹ 2z†cosvDt

2

#

…22† The stability analysis is performed for a linear initial value problem, which is extended to the consideration of the above described effects. The analysis steps for the linear case are taken from [2]. The exact transfer matrix A is given in Eq. (11). The transformation into the Jourdan’s form A ˆ PJP ⫺1 yields a diagonal matrix J, see [21]. As can be seen from Eq. (17) the power m of the matrix A is the same as the power m of the diagonal elements in the matrix J, i.e. Am ˆ PJm P⫺1 . The eigenvalues of the matrix A are the diagonal elements of the matrix J. It can be stated, that the absolute value of the maximum eigenvalues has to be less than one. This must hold for every numerical algorithm if it is said to be unconditionally stable. The eigenvalues of the exact transfer matrix A are l1;2 ˆ e…⫺zv^ivd †t and the absolute value pof  the eigenvalues is computed to r ˆ max兩l1;2 兩 ˆ A2 ˆ e⫺zvt , which is called spectral radius. The presented semi-analytic procedure can be analysed with respect to its numerical behaviour. In [7] the equations for the numerical dissipation, the numerical dispersion and the accuracy are given for the approximate transfer matrix. When considering a linear modal oscillator it is confirmed, that the present BEM-type method is exact. Most of the direct time integration algorithms show convergence of second-order. The numerical dissipation is the sum of the modal damping of the oscillator and an additional damping owing to the assumptions of the algorithm. The numerical dispersion is characterized by the resulting

507

relative period error. The accuracy is computed from the local error for one time step. The advantage of the present method for the linear system results owing to the application of the exact transfer matrix, so the numerical dissipation and dispersion and the local error are zero for the linear symmetric system. The numerical stability of Eq. (21), when used in a timestepping procedure, depends on the eigenvalues l of the ~ which for unconditional stability have transition matrix A, to satisfy the condition 兩l兩 ⱕ 1. The absolute values 兩l兩 for the three cases of Eq. (20) are shown in Figs 1(a)–(c) as a function of Dt=T, where T denotes the period of the undamped vibration. Also shown are the exact solutions. It is seen that only case (b) gives an unconditionally stable algorithm. The errors accumulated within the time-interval while compared to the exact solution become: ! 1 ⫺ 2z2 _ zv4 Dt2 ⫹ O…Dt3 †; ka …Dt† ˆ zX0 ⫺ …23a† X0 v 3 ! 1 ⫺ 4z2 _ zv4 Dt2 ⫹ O…Dt3 †; kb …Dt† ˆ zX0 ⫺ X0 2v 3

…23b†

  5 3 4z zv Dt kc …Dt† ˆ …1 ⫺ 4z2 †X0 ⫹ …1 ⫺ 2z2 †X_ 0 v 36 ⫹ O…Dt4 †:

…23c†

It is seen that case (c) is third order accurate, while for the cases (b) and (c) the accuracy is only of second-order. (Note that these results are purely artificial, as in practice one would treat linear oscillations by means of the exact transition matrix.) Similar estimates can be obtained for the errors in the vibration period and for the phenomenon of numerical damping. E.g. in case (b) one gets T~ b ⫺ T …3 ⫺ 4z2 †z2 v2 Dt2 ⫹ O…Dt4 †; ˆ⫺ T 12…1 ⫺ z2 †

…24a†

z~ b ⫺ z …1 ⫺ 2z2 †v2 Dt2 ⫹ O…Dt4 †; ˆ⫺ 12 z

…24b†

see Holl [7] for details. Here for the above algorithm some additional results are presented for the case of the equation of motion X ⫹ 2zvX_ ⫹ v2 X ˆ ⫺av2 X:

…25†

The exact transfer matrix formulation for this case can be taken from Eqs. (11) and (15) to (17) with some modifications with respect to v and z . Eq. (21) is calculated using the convolution integral of Eq. (10) and inserting Eq. (20b) for the generalized coordinate. The approximate transfer matrix ~ results after some algebraic manipulations. Fig. 2 shows A the absolute values of the eigenvalues of the resulting algorithm with the parameters a ˆ 0:3 and z ˆ 0:1. An analogous analysis was performed for the Newmark method [22].

508

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Fig. 3. (a)–(c) Phase-plane of the non-dimensional error for the Duffing oscillator for three methods.

Based on the characteristic equations, see [23,24], an application to Eq. (25) was computed. In Fig. 2 the result for the absolute value of the eigenvalues for the unconditionally stable Newmark method can be seen for g ˆ 1=2 and b ˆ 1=4. The result for the accuracy is computed from the local error measure, see [2,7], using a series expansion. The result for the above presented algorithm is ! X_ 0 v4 Dt2 k ˆ …1 ⫺ a†X0 ⫹ 2z a…1 ⫺ zvDt† ⫹ O…Dt4 †; v 12 …26† where X0 and X_ 0 are the initial values at the beginning of the

time step. When performing the analogous analysis for the Newmark method the local error for z ˆ 0 results to   3   X_ 0 1 v Dt g 1 ⫺ ⫺b ⫺g ⫹ k ˆ …1 ⫺ a† v 2 2 2 6  …1 ⫺ a†2 X0 v4 Dt2 ⫹ O…Dt3 †:

…27†

We see, that g ˆ 1=2 gives an algorithm, which has a convergence of second-order. This result for the convergence does not change for the linear case a ˆ 0, whereas from Eq. (26) it is obvious, that for a ˆ 0 the local error disappears, i.e. k ˆ 0 Some detailed study was done in [8] for analysing a gyroscopic system and in [23] for nonclassical damping, where

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

509

Fig. 4. Resulting motion and error of the computation for a linear pendulum with time-dependent mass.

different approximations for the generalized coordinates are analysed. An extension to systems with time-varying mass is presented in [16] and [24]. 4. Examples 4.1. Free vibrations of an undamped Duffing oscillator The above BEM-type time-integration scheme is applied to free vibrations of a Duffing oscillator, q ˆ ⫺bX 3 with initial conditions bX02 ˆ 4v2 and X_ 0 ˆ 0. Scheme (b) of Eq. (20b), which has been shown to be unconditionally stable in the examples of Section 3, is used in a time-stepping procedure with vDt ˆ 2p=100. The increments of the state were calculated solving a cubic equation. Fig. 3 shows the error in the displacement 1…t† ˆ …wcalc ⫺ wexact †=w0 versus the error in velocity s…t† ˆ …wcalc ⫺ wexact †=w0 v for the present method, the Newmark method and for the Cubic Hermitian Method by Argyris [25]. It is seen that the present method is more accurate, while the CPU-time of the present method was smaller by a factor of 1/2. Comparisons with other direct time-integration schemes turned out to give an even better benefit of our method. An iterative and incremental computation is necessary within each time step for the nonlinear analysis. For the special cases of a nonlinear restoring force which is a polynom of up to the third order the solution can be found without any iteration, and we are able to solve the resulting nonlinear equation in one step. 4.2. Examples of systems with time-dependent mass

m…t† ˆ m0 e⫺avt ; …29† p where v ˆ g=l is the eigenfrequency of the linearized system with constant mass. The chosen system has the advantage, that the exact solution is known for this special function for the variable mass, so that an analysis of the computed error is possible easily. Insertion of Eq. (29) in the equation of motion (28) leads to the new differential equation for the small amplitudes

w ⫺ avw_ ⫹ v2 w ˆ f^ ˆ

f : m

…30†

The solution of the linearized equation of motion (30) for the given inhomogeneous initial conditions of w_ 0 =w0 v ˆ 2 is known. For the nonlinear case of Eq. (28) the solution could be computed by means of elliptic integrals. Additionally the force f with the amplitude fstep =v2 w0 m0 ˆ 1 is applied as a step function. The parameter of the time-depenp dent mass is a ˆ 1=10. With v ˆ v 1 ⫺ …a=2†2 the dimensionless solution is given by 2 3 a w_ ⫹ 0 6 7 w…t† 2 w0 v sin…v t†7 ˆ eavt=2 6 cos…v t† ⫹ q 4 5 ÿ
w0 1 ⫺ a=22 39 a > sin…  = v t† fstep e 7 ⫺avt=2 6 2 1⫺e ⫹ 2 4cos…v t† ⫹ p 5 : v w0 m0 > : ; 1 ⫺ …a=2†2 > avt

8 > <

2

…31†

As a further example a mathematical pendulum of the length l with time-dependent mass m(t) is considered. The equation of motion of this system is

2…mw_ † g ⫹ m sinw ˆ f : 2t l

It is assumed, that the mass of the pendulum decays in accordance with

…28†

The above semi-analytic algorithm is used to compute the solution numerically. In a first approach we use Eq. (30). As can be seen from Fig. 4(a) the calculations are accurate and show only some differences to the exact solution as a result of the chosen large time steps of Dt=T ˆ vDt=2p ˆ 1=10:

510

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Fig. 5. Dimensionless mass ratio x …t† ˆ m…t†=m0 .

Fig. 7. Solution for x …t† of Fig. 5 in the phase plane.

Five periods are considered. In a second approach Eq. (28) is reformulated in a suitable form for the application of the presented BEM-type algorithm: m0 w ⫹ m0 v2 w ˆ

f _ w_ ⫺ …m…t† ⫺ m0 †w ⫺ …m…t† ⫺ m…t† m…t† ⫺ m 0 † v w: 2

…32†

The approximation of Eq. (20b) is chosen and the resulting time-integration method is computed using the same time step again. Fig. 4(a) shows the computed results. Depending on the chosen time step, the calculated solution shows some error when compared to the exact solution. In Fig. 4(b) the calculation error 1 ˆ …wcalc ⫺ wexact †=w0 is drawn for this time step. A detailed study showed, that the smaller the time step the smaller is the error in the result. The computed results for the Newmark method are drawn in Fig. 4(a) and (b) additionally, which shows higher errors, when using the same time step. 4.3. Exponentially varying mass: Comparison to exact results In this section, an exact solution of an undamped linear oscillator, z ˆ 0, b ˆ 0, with a particular time-varying mass m…t† ˆ m0 x …t† is compared with the numerical solution. Using the abbreviation

x …t† ˆ 1 ⫹ x…t†;

…33†

a non-dimensional time scale t is introduced by Zt r c tˆ dt: x m 0 0  …t†

…34†

Note that for time-invariant systems, x (t) ˆ 0, the definition of Eq. (34) results in the natural time t ˆ v t of the background oscillator. Considering free vibrations, f ˆ 0, and using Eq. (34), the equation of motion, Eq. (3), transforms to X 00 …t† ⫹

x 0 …t† 0 X …t† ⫹ X…t† ˆ 0; 2x …t†

…35†

where prime denotes differentiation with respect to the nondimensional time t . In order to find a solution of Eq. (35), we introduce a multiplicative representation of the form X…t† ˆ Z…t†U…t†: Eq. (35) then takes on the following form: ! Z 0 …t† x 0 …t† 00 ⫹ U 0 …t† U …t† ⫹ 2 Z…t† 2x …t† ! Z 00 …t† 1 Z 0 …t† x 0 …t† ⫹ U…t† ˆ 0: ⫹ 1⫹ Z…t† 2 Z…t† x …t†

…36†

…37†

Z is chosen in such a way that the term which is linear in U 0 …t† vanishes, i.e. 2

Z 0 …t† x 0 …t† ⫹ ˆ 0: Z…t† 2x …t†

…38†

A solution of Eq. (38) follows as: 1 Z…t† ˆ p4  : x …t†

…39†

Using Eqs. (38) and (39), Eq. (37) simplifies to the form: "  # d2 1 00 …40† U…t† ˆ 0: U …t† ⫹ 1 ⫺ Z…t† 2 dt Z…t† Fig. 6. Displacement and velocity for x …t† of Fig. 5.

A closed form solution of Eq. (40) is possible for various functions x . In the following, an exponential variation of the

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

511

A comparison of the numerically computed results is shown in Fig. 8 for different time steps. The convergence of the solution with smaller time steps is clearly demonstrated. It can be seen that the convergence of the present method is of the order of two, which states that a reduction of the applied time step by a factor of 2 gives a reduction of the error of the solution by a factor of 4.

4.4. Rotordynamic system with variable mass

Fig. 8. Solutions with different time steps.

mass is studied:

x …t† ˆ Ke⫺2avt ;

…41†

where K and a are constants. Inserting Eq. (41) into Eqs. (39) and (40), and performing the further transformation r m0 K avt ˆ eavt ; …42† xˆ1⫹ c we finally end up with a Bessel type differential equation for U(x): "  # 2 x 2 3 2d U U ˆ 0: …43† ⫺ x c 4 dx2 r m0 K . Following Kamke [26], the final soluwith c ˆ av c tion is      1 x x X…x† ˆ p4  x C1 J1 ⫹ C2 Y1 ; …44† c c K where J1 and Y1 are the Bessel functions of integer order one. This solution is discussed for the parameters c ˆ a ˆ 1=20. Fig. 5 shows the time evolution of the dimensionless mass ratio x …t† of the oscillator. The resulting vibrations are shown in Fig. 6 for the dimensionless displacement and velocity and in Fig. 7 in the phase plane for the initial conditions X(t ˆ 0) ˆ X0 and _ ˆ 0† ˆ X_ 0 ˆ 0. X…t

As a further example the nonlinear rotordynamic system of Fig. 9 is analysed. The system is discretized with 20⬚ of freedom using four rotating beam elements and three rotor elements. The elements of the stiffness matrix which correspond to the degrees of freedom at the bearing are k11 ˆ k22 ˆ 10 9 N/m and k12 ˆ kÿ21 ˆ 10
7 N/m. A cubic nonlinearity of the form rN ˆ kii xi b xi =xstat 2 with i ˆ 1,2,13,14 is chosen for the bearings with the parameter b ˆ 0.1. For homogeneous initial conditions and a rotating speed of 1500 rpm a sudden excentricity is introduced with the force F ˆ 5 × 10 4 N. For this case we get the computed results of Figs. 10–13. The time variable mass was assumed to be M…t† ˆ M0 ‰1 ⫹ xsin…nt†Š, where n is the angular velocity and x denotes to the variable portion of the mass, which is indicated in Figs. 10–13. Figs. 10 and 11 show the motion of the center of the journal in the bearings in a plane perpendicular to the undeformed rotor axis. The effect of the nonlinear restoring force is shown for these two positions for b ˆ 0 and b ˆ 0.1. The dimensionless displacements are scaled using the corresponding static deformations for b ˆ 0 owing to the weight itself. Figs. 12 and 13 show the motion of the center at the left and the right rotor position. It is seen, that the influence of the variable mass considerably influences the motion of the rotor. A study of the application of modal reduction showed that the solution gave only small errors when considering at least five generalized coordinates. The higher errors occured at the position of the bearings. In Table 1 the relative computational time is shown for this problem. In correspondance to the results of [8] the Newmark method used about two times larger computation time than the present method without modal reduction.

Fig. 9. Rotordynamic model.

512

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513

Fig. 10. Motion of the center at the left bearing. Fig. 13. Motion of the center at the right rotor.

Table 1

Modal DOF

relative CPU-time in %

20 15 10 5 3

100 85 60 36 27

5. Conclusion

Fig. 11. Motion of the center at the right bearing.

The presented semi-analytic algorithm makes use of the methods of linear structural dynamics and applies them in a consistent way to the analysis of transient vibrations of nonlinear rotordynamic systems. It is demonstrated, that the resulting algorithm is efficient and accurate. A local error is given for a special problem and the absolute values of the eigenvalues are shown as a measure of stability. The computation time for solving example problems shows an advantage for the present method when compared to direct time integration methods like the Newmark method. When using the same time step, the error in the computed results, which is compared with a converged solution, showed an advantage as well.

References

Fig. 12. Motion of the center at the left rotor.

[1] Zienkiewicz OC, Taylor RL. The finite element method, Vol. 2, 4th ed. New York: McGraw Hill, 1991. [2] Hughes TJR. Analysis of transient algorithms with particular reference to stability behavior. In: Belytschko T, Hughes TJR, editors. Computational Methods for Transient Analysis, North-Holland: Elsevier, 1983, pp. 67–155. [3] Irschik H, Ziegler F. Dynamics of linear elastic structures with selfstress: a unified treatment for linear and nonlinear problems. ZAMM 1988;68(6):199–205.

H.J. Holl et al. / Engineering Analysis with Boundary Elements 23 (1999) 503–513 [4] Irschik H, Ziegler F. Dynamic processes in structural thermo-viscoplasticity. Applied Mechanics Review 1995;48(6):301–315. [5] Fotiu P, Irschik H, Ziegler F. Material science- and numerical aspects in the dynamics of damaging structures. In: Schue¨ller GI, editor. Structural Dynamics. Berlin: Springer, 1991, pp. 235–255. [6] Fotiu P, Irschik H, Ziegler F. Modal analysis of elastic-plastic plate vibrations by integral equations. In: Engineering Analysis with Boundary Elements, Vol. 14. Amsterdam: Elsevier, 1994, pp. 81–97. [7] Holl HJ. Ein effizienter Algorithmus fu¨r nichtlineare Probleme der Strukturdynamik mit Anwendung in der Rotordynamik. Dissertation, University of Linz, 1994. [8] Holl HJ. An effizient time domain formulation for nonlinear rotordynamic systems using modal reduction. Proc. of the 13th IMAC, Nashville, 1995, pp. 856–865. [9] Holl HJ, Irschik H. A substructure method for the transient analysis of nonlinear rotordynamic systems using modal analysis. In: Proc. of the 12th IMAC, Honolulu, 1994, pp. 1638–1643. [10] Holl HJ, Irschik H. Eine effiziente Substrukturmethode fu¨r transiente Probleme der nichtlinearen Rotordynamik. In: Irretier H, Nordmann R editors. SIRM III, Berlin: Springer, 1995, pp. 297–305. [11] Nardini D, Brebbia, CA. Transient dynamic analysis by the boundary element method. In: Brebbia CA, Futagami T, Tanaka M, editors. Boundary Elements. Berlin: Springer, 1983, pp. 719–730. [12] Manolis GD, Beskos DE. Boundary element methods in elastodynamics. Unwin Hyman. 1988. [13] Mukherjee S. Boundary element methods in creep and fracture. Barking: Applied Science Publisher, 1982. [14] Ziegler F. Mechanics of solids and fluids, Berlin: Springer, 1991. [15] Holl HJ, Irschik H. A boundary element method for the time-integration of nonlinear dynmic systems. In: Morino L, Wendland WL editors. Boundary Integral Methods for Nonlinear Problems. Dordrecht: Kluwer Academic Publishers, 1997, pp. 101–107.

513

[16] Holl HJ, Belyaev AK, Irschik H. Simulation of the Duffing-oscillator with time-varying mass by a BEM in time. In: Topping, BHV, editor. Advances in Boundary Element Methods Civil Comp, 1996, pp. 113– 120. [17] Holl HJ, Belyaev AK, Irschik H. A time integration algorithm for nonlinear problems in rotordynamic systems with time-varying parameters. In: Muszynska A, editor. Proc. of the 7th International Symposium on Transport Phenomena and Dynamics of Rotating Machinery, Vol. B. Bird Rock Publishing House, 1998, pp. 930–940. [18] Brebbia CA, Dominguez J. Boundary Elements, 2nd edition, New York: McGraw Hill, 1992. [19] Maddox NR. On the number of modes necessary for accurate responce and resulting forces in dynamic analysis. Journal of Applied Mechanics 1975;42:516–517. [20] Waller H, Schmidt R. Schwingungslehre fu¨r Ingenieure. B.I. Wissenschaftsverlag, 1989. [21] Bathe K-J. Finite element procedures in engineering analysis. Englewood Cliffs, NJ: Prentice Hall, 1982. [22] Newmark NM. A method of computation for structural dynamics. Journal of the Engineering Mechanics Division 85 ASCE EM 1959;3(2094):67–94. [23] Holl, H.J. An time integration algorithm for time-varying systems with nonclassical damping based on modal methods. In: Proc. of the 15th IMAC, Orlando, 1997, pp. 1558–1564. [24] Holl HJ, Irschik H. Integration of a nonlinear dynamic system with time-varying mass using a boundary element formulation in time. In: Strganac T, Tinker M, editors. Proc. of the Dynamics Specialist Conference, Reston: AIAA, Salt Lake City, 1996, pp. 297–305. [25] Argyris J, Mlejnek HP, Dynamics of structures. Texts on Computational Mechanics, Vol. 5, North-Holland: Elsevier, 1991. [26] Kamke, E., Differentialgleichungen. Band 1, 10. Auflage, Teubner Verlag, 1983.