- Email: [email protected]
Nonlinear vibrations in a covered roll system with viscoelastic contact
Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Nonlinear vibrations in a covered roll system with viscoelastic contact L. Yuan *, V.-M. Järvenpää Department of Mechanics and Design, Tampere University of Technology, Korkeakoulunkatu 6, P.O. Box 589, 33101 Tampere, Finland
a r t i c l e
i n f o
Article history: Received 1 October 2008 Received in revised form 25 October 2008 Accepted 25 October 2008 Available online 14 November 2008 PACS: 05.45.a Keywords: Viscoelastic cover Time delay Contact vibrations Nonlinear dynamics Perturbation analysis
a b s t r a c t A theoretical analysis on the contact vibration problem of a nip rolling system composed with two paper machine rolls is presented. It is modelled at the contact region by a twodegree-of-freedom (DOF) mass-spring-damper system. The rolling contact force in a case of viscoelastic polymer covered roll is formulated by integration of the contact stress. Regenerative chatter effects to the stability of the vibrations are considered. The perturbation technique is applied through the multi-scale method to calculate the nonlinear normal form of the governing equations to determine the stability behavior of the system. The numerical results in time domain are consistent with the bifurcation diagrams. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction The contact vibration problem of a roll nip system is presented. Roll nips are used for finishing the surface of processing paper and extensively exist in large scale paper machines of modern paper manufacturing industry. One nip unit consists of two rolls rotating with contact. The traditional hard-nip machine calenders, in which a paper web is calendered to a uniform thickness, cause the variations in paper surface properties [1]. In supercalenders, the nip is typically formed by one hard metal roll and one soft roll coated with polymer layer, so-called soft-nip. In this paper, the soft-nip will be studied. The soft cover may distribute the contact pressure in the nip and produce a paper with more constant density rather than constant caliper of the paper. However, the polymer material has such strong viscoelastic characters that can considerably affect the contact stresses. Consequently a regenerative vibration condition evidently originates from the polymer cover due to its finite viscoelastic recovery time. How to deeply understand the vibration mechanics of the nip contact and substantially mitigate the harmful vibrations (including noises accordingly) are still challenging tasks in the field of the paper manufacturing engineering. In the related literatures, many of researches are concentrated on the contact stress between elastic or viscoelatic bodies in different methods [2–11]. These are mostly for static analysis. On the other hand, some interesting works on dynamic modeling of rolling contact intrigue us. For example, Nayak [12] in earlier 1972 presented a rolling contact dynamic model with two elastic bodies solving by harmonic balance method. Much later in 1997, using the same model as Nayak, Narayana and Sekar [13] solved its equation by Fourier–Galerkin–Newton method and obtained the bifurcation structure of contact vibrations. Gray and Johnson [14] investigated the response of the rolling elastic bodies through the theory and techniques of random vibration. They claimed that contact load and damping are significant parameters. Turner [15] studied nonlinear
* Corresponding author. E-mail address: lihong.yuan@tut.fi (L. Yuan). 1007-5704/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2008.10.028
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
3171
vibrations with Cantilever–Hertzian contact boundary conditions in a beam application. In recent researches, the dynamic contact force has been analyzed and might be nonlinear from the contact stiffness [16,17]. But it is still needed to discovery more for the aspect that the contact force can be applied into a dynamic model. Such kind of problems belongs to contact vibrations. Authors have studied on the time delay dynamic system [18] whereas some tests and observations on a nip installation have been done and analyzed [19]. Furthermore, to model and predict the vibration motion of the nip system, formulating of the contact force based on the contact theory is needed to explore the soft-nip contact mechanics. The intent of this work is to develop the model of the contact vibrations of the nip system. Firstly, we will present the contact mechanics involved in the deformation of the corrugations on the polymer covered roll surface. The normal contact force is integrated from the pressure distribution between the contact areas. The nonlinear relation between the normal force and the penetration is examined. The variation of the contact force causes the vibrations and appears as a nonlinear term in the mathematical equations of the model. Secondly, we will establish the equations of motion for the dynamic contact vibrations essentially nonlinear, which is of also theoretical interest. Moreover, the time delay term owing to the viscoelastic polymer cover is considered and makes the entire system under a regenerative vibration condition. These two reasons might lead to severe difficulties in the application of approximate techniques for the analysis of these nonlinear time delay differential equations including the solutions and stability. As refer to the stability of time delay system is more complex and interesting especially when it involves nonlinear elements. The analysis of the Hopf bifurcation in different models can be carried out by using analytical methods like multi-scale, harmonic and Floquet theory [20–25]. From the perturbation theory [26], the multi-scale method in this work will be executed for analyzing the mathematical model to obtain the normal form and the amplitude of the limit cycle of the motion. 2. Mathematic modeling of nip contact For this model, there are several assumptions. The time delay is constant, there is no loss of contact in the nip, the small torsional vibration is neglected, the bearing effects are not considered and the temperature generated from the rolling contact is not included. The metal rolls are big (the contact length is 4.4 m and roll radius 2.75 m) and the rolling speed is high (500 m/min). To establish the dynamic vibration model of the nip system, the normal contact force between the rigid roll and the polymer covered roll should be formulated at first. The distribution function p(z) of the normal contact pressure between the two rolls in Fig. 1 is unsymmetrical due to the viscoelastic material behavior on layered contact [27]
pðzÞ ¼
Ea2 1 z2 z þ bnð1 þ nÞð1 eð1þz=aÞ=n Þ : 1 2 bn 1 þ a Rs 2 a
ð1Þ
The static normal force over the contact area can be computed by
P¼
Z
b
pðzÞdz:
ð2Þ
a
We get
" # 3 2 Ea3 1 b b 1 1 2 2 2 bþa an ; P¼ þ a þ abn b 2 bn þ bbn þ abn ð1 þ nÞe n 3 2 Rs 2 6a 2a
ð3Þ
where the equivalent radius R = R1R2/(R1 + R2), Deborah number n = vwebsrelaxation/a, and a R. The dimension coefficients R1 and R2 are radiuses of the two rolls, the sum of the contact widths (a + b) is in the z-direction of the loaded area, s is the thickness of the thin polymer cover layer of the upper roll and E is equivalent elastic modulus [27]. vweb is the speed of the paper roll 1, metal
s polymer cover
R1 p(z)
δ R2
x a
b
roll 2, metal z
Fig. 1. The pressure distribution in the contact region.
3172
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
web through the nip, which corresponds to the running speed of the rolls; srelaxation is the viscoelastic relaxation time of the polymer material and b is the viscoelastic material parameter. The compressive strain ec in an element of the contact at horizontal z-direction is given by
ec ¼ ðd z2 =2RÞ=s; a 6 z 6 b;
ð4Þ
where d is the maximum depth of the penetration of the covered roll. By assuming b a in a symmetric contact, which is the typical case in practical applications
P¼
E 2 2 2 þ bn2 1 þ en þ bn3 1 þ en a3 Rs 3
ð5Þ
or one can rewrite like
P ¼ hF U a3 ;
ð6Þ
where the parameter
FU ¼
E 2 2 2 þ bn2 1 þ en þ bn3 1 þ en R 3
ð7Þ
And h is the reciprocal of the thickness of the polymer cover s, h = 1/s. Fig. 2 illustrates the equivalent spring-mass system for the nip contact rolls. Let x1s = x1(t s). Symbol s is the delay time relevant to the rotational speed of the polymer covered roll. Now we need to seek the relation between the contact width a and the dynamic displacements. In the formula Eq. (4), the strain is zero at z=a; thus one can get
a2 ¼ 2R d:
ð8Þ
The maximum depth of the penetration of the roll d in the dynamical model should be related to the displacements x1, x2 and expressed like
d ¼ d0 þ x1 x2 þ DðtÞ;
ð9Þ
where d0 is the nominal penetration in the contact region and the external excitation is D(t) = Dsin(Xt) from the paper web, which is assumed as zero here.
c1
k1 x1 roll 1 m1
δ0
Δ (t)
Pdynamic
c nip
x2
m2 roll 2 c2
k2
Fig. 2. The equivalent spring-mass system.
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
3173
Then the width of the contact a can be solved as the function of x1 and x2 as
a2 ðtÞ ¼ 2Rðd0 þ x1 x2 Þ;
ð10Þ
and the formula Eq. (10) can be transferred to
aðtÞ ¼
ffi pffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x1 x2 2Rd0 1 þ : d0 d0
ð11Þ
Now the normal force P in formula Eq. (6) becomes the function of x1 and x1s due to the formula Eq. (10), thereby the governing equations of the motion of the rolling system are expressed as
€x1 þ 2f1 x_ 1 þ p21 x1 ¼ Nðx1 ; x1s ; x2 ; x2s ; tÞ=m1 ; €x2 þ 2f2 x_ 2 þ p22 x2 ¼ Nðx1 ; x1s ; x2 ; x2s ; tÞ=m2 ;
ð12Þ ð13Þ
where
2fi ¼ ci =mi ;
p2i ¼ ki =mi ;
and i ¼ 1; 2:
ð14Þ
Considering the time delay effect in formula Eq. (6), the entire dynamic contact force N including the nip damping in the contact is
_ _ Nðx1 ; x1s ; x2 ; x2s ; tÞ ¼ Pdynamic þ h cnip dðtÞ ¼ h F U ½a3 ðtÞ ce a3 ðt sÞ þ h cnip dðtÞ:
ð15Þ
The relaxation parameter ce is related to the time delay and the polymer material characteristics. 3. Perturbation analysis To do the perturbation analysis, the nonlinear cubic term a3(t) obtained from Eqs. (10) and (11) can be expanded as a Taylor series
ffi pffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x1 x2 2Rd0 1 þ d0 d0 " 2 3 # pffiffiffiffiffiffiffiffiffiffiffi 1 x1 x2 1 x1 x2 1 x1 x2 ; þ ffi 2Rðd0 þ x1 x2 Þ 2Rd0 1 þ 2 d0 d0 8 d0 d0 16 d0 d0
a3 ðtÞ ¼ 2Rðd0 þ x1 x2 Þ
ð16Þ
where higher orders are omitted and xd10 xd20 6 1 must be satisfied. Now the equations of the model Eqs. (12) and (13) evidently consist of nonlinear terms due to the formula Eq. (16). In this system, the masses of the two rolls are close to equal. For simplicity, let m1 = m2. It is assumed that the term FU is constant. Then we can give the notes like
pffiffiffiffiffiffiffiffiffiffiffi F ¼ F 1 ¼ F 2 ¼ ðF U =m1 Þ 2R 2Rd0 ;
cn ¼ cn1 ¼ cn2 ¼ cnip =m1 :
ð17Þ
The perturbation analysis can be performed into the delay type equations using the method of multiple scales [28]. In this method, a fast time-scale T0 = t and a slow time-scale T1 = e2t are introduced. We seek a second-order expansion for the Eqs. (12 and 13) in the form
xi ¼ exi0 ðT 0 ; T 1 Þ þ e2 xi1 ðT 0 ; T 1 Þ þ e3 xi2 ðT 0 ; T 1 Þ þ Oðe4 Þ;
i ¼ 1; 2:
ð18Þ
The time derivatives are
@xi0 @xi1 @xi2 @xi0 þ e2 þ e3 þ 2e3 þ ; @T 0 @T 0 @T 0 @T 1 @ 2 xi0 @ 2 xi1 @ 2 xi2 @ 2 xi0 €xi ¼ e þ e2 þ e3 þ 2e3 þ 2 2 2 @T 0 @T 1 @T 0 @T 0 @T 0 x_ i ¼ e
ð19Þ ð20Þ
Then we perturb the contact parameter (relevant to the effective thickness of the polymer layer), let
h ¼ hc þ e2 h1 ;
ð21Þ
where hc is the reciprocal of thickness value s at the stability boundary. Substituting the Eqs. (18)–(21) into the governing Eqs. (12) and (13), and equating the coefficients of like power of e, we obtain results as following. Order e1: for the first governing Eq. (12)
@ 2 x10 @T 20
þ 2f1
@x10 1 1 þ p21 x10 ¼ F hc ðx10 x20 Þ þ hc ðx10 x20 Þ hc ce ðx10s x20s Þ hc ce ðx10s x20s Þ 2 2 @T 0 @x10 @x20 hc c n ; @T 0 @T 0
ð22aÞ
3174
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
for the second Eq. (13)
@ 2 x20 @T 20
þ 2f2
@x20 1 1 þ p22 x20 ¼ F hc ðx10 x20 Þ þ hc ðx10 x20 Þ hc ce ðx10s x20s Þ hc ce ðx10s x20s Þ 2 2 @T 0 @x10 @x20 þ hc c n : @T 0 @T 0
ð22bÞ
Order e2:
@ 2 x11 @T 20
@ 2 x21 @T 20
þ 2f1
þ 2f2
@x11 1 1 1 þ p21 x11 ¼ F hc ðx10 x20 Þ2 þ hc x10 x20 þ ðx10 x20 Þ2 ðx11 x21 Þ þ 2 8d0 2d0 @T 0 1 1 1 ðx10s x20s Þ2 hc ce x11s x21s þ ðx10s x20s Þ2 ðx11s x21s Þ þ hc ce 2 8d0 2d0 @x11 @x21 þh1 d0 ð1 ce Þ hc cn ; @T 0 @T 0
@x21 1 1 1 þ p22 x21 ¼ F hc ðx10 x20 Þ2 þ hc x10 x20 þ ðx10 x20 Þ2 ðx11 x21 Þ þ 2 8d0 2d0 @T 0 1 1 1 hc ce ðx10s x20s Þ2 hc ce x11s x21s þ ðx10s x20s Þ2 ðx11s x21s Þ þ 2 8d0 2d0 @x11 @x21 h1 d0 ð1 ce Þ þ hc cn : @T 0 @T 0
ð23aÞ
ð23bÞ
Order e3:
@ 2 x12 @T 20
þ2
@ 2 x10 @x12 @x10 1 1 þ 2f1 þ ðx10 x11 þ x20 x21 x10 x21 x11 x20 Þ ðx12 x22 Þ þ p21 x12 ¼ F hc 2 4d0 @T 0 @T 1 @T 0 @T 1 ! 1 3 2 2 3 þ ðx 3x x þ 3x x x Þ 10 10 20 20 10 20 16d20 1 þ hc x12 x22 ðx10 x11 þ x20 x21 x10 x21 x11 x20 Þ d0 ! 1 3 2 2 3 2 ðx10 3x10 x20 þ 3x20 x10 x20 Þ 8d0 1 1 ðx10s x11s þ x20s x21s x10s x21s x11s x20s Þ hc ce ðx12s x22s Þ 2 4d0 ! 1 3 2 2 3 ðx 3x x þ 3x x x Þ þ 10s 10s 20s 20s 10s 20s 16d20 hc ce ðx12s x22s þ
1 ðx10s x11s þ x20s x21s x10s x21s x11s x20s Þ d0
) 2 3 2 2 3 ðx 3x x þ 3x x x ÞÞ þ h x Þð1 c Þ ðx 1 10 20 e 10s 10s 20s 20s 10s 20s 3 8d20 @x10 @x12 @x20 @x22 @x10 @x20 þ þ c n hc þ cn1 hc h1 cn1 : ð24Þ @T 1 @T 0 @T 1 @T 0 @T 0 @T 0
1
The solutions of the Eqs. (22a) and (22b) can be expressed as
x10 ¼ A1 ðT 1 Þeixc T 0 þ A1 ðT 1 Þeixc T 0 ;
x20 ¼ A2 ðT 1 Þeixc T 0 þ A2 ðT 1 Þeixc T 0 ;
x10s ¼ A1 ðT 1 Þeixc ðT 0 sÞ þ A1 ðT 1 Þeixc ðT 0 sÞ ;
x20s ¼ A2 ðT 1 Þeixc ðT 0 sÞ þ A2 ðT 1 Þeixc ðT 0 sÞ ;
ð25Þ
where xc is the chatter frequency at the Hopf bifurcation [20]. Also let
x11 ¼ P1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 ; x21 ¼ R1 e2ixc T 0 þ R2 e2ixc T 0 þ R3 :
ð26Þ
3175
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
By substituting the Eqs. (25) and (26) into the Eq. (23a), we obtain
4x2c ðP1 e2ixc T 0 þ P2 e2ixc T 0 Þ þ 4f1 ixc ðP1 e2ixc T 0 P2 e2ixc T 0 Þ þ p21 ðP1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 Þ
1 1 ¼ F hc ðA1 eixc T 0 þ A1 eixc T 0 A2 eixc T 0 A2 eixc T 0 Þ2 ðP 1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 R1 e2ixc T 0 R2 e2ixc T 0 R3 Þ 2 8d0 1 2ixc T 0 2ixc T 0 2ixc T 0 2ixc T 0 þ P2 e þ P 3 R1 e R2 e R3 Þ ðA1 eixc T 0 þ A1 eixc T 0 A2 eixc T 0 A2 eixc T 0 Þ2 Þ hc ððP1 e 2d0 1 1 þ hc ce ðP1 e2ixc ðT 0 sÞ þ P2 e2ixc ðT 0 sÞ þ P3 R1 e2ixc ðT 0 sÞ R2 e2ixc ðT 0 sÞ R3 Þ ðA1 eixc ðT 0 sÞ þ A1 eixc ðT 0 sÞ 2 8d0 A2 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ Þ2 þ hc ce ððP 1 e2ixc ðT 0 sÞ þ P2 e2ixc ðT 0 sÞ þ P3 R1 e2ixc ðT 0 sÞ R2 e2ixc ðT 0 sÞ R3 Þ 1 þ ðA1 eixc ðT 0 sÞ þ A1 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ Þ2 Þ h1 d0 þ h1 ce d0 hc cn 2ixc ðP1 e2ixc T 0 P2 e2ixc T 0 Þ 2d0 þ hc cn 2ixc ðR1 e2ixc T 0 R2 e2ixc T 0 Þ:
ð27Þ
For Eq. (27), equating each of the coefficients of exp(2ixcT0), exp(2ixcT0) and exp(0) on both sides, we obtain
P1 ¼ /1 ðA1 A2 Þ2 ; P 2 ¼ /2 ðA1 A2 Þ2 ;
P3 ¼ /3 ðA1 A2 ÞðA1 A2 Þ
ð28Þ
and
R1 ¼ h1 ðA1 A2 Þ2 ;
R2 ¼ h2 ðA1 A2 Þ2 ;
R3 ¼ h3 ðA1 A2 ÞðA1 A2 Þ;
ð29Þ
where
D2 C 3 þ D3 C 1 ; D1 C 1 D2 C 2 D9 C 8 þ D7 C 9 h3 ¼ ; D7 C 7 D8 C 8
/1 ¼
/2 ¼
D5 C6 þ D6 C4 ; D4 C4 D5 C5
/3 ¼
D8 C 9 þ D9 C 7 ; D7 C 7 D8 C 8
h1 ¼
D3 C 2 þ D1 C 3 ; D1 C 1 D2 C 2
h2 ¼
D6 C5 þ D4 C6 ; D4 C4 D5 C5 ð30Þ
and the groups of the parameters
D1 ¼ 4x2c þ 4f1 ixc þ p21 þ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; D2 ¼ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; D3 ¼ 38 Fhc ð1 ce e2ixc s Þ; D4 ¼ 4x2c 4f1 ixc þ p21 þ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; D5 ¼ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; D6 D7 D8 D9
¼ 38 Fhc ð1 ce e2ixc s Þ;
C1 C2 C3 C4 C5 C6 C7 C8 C9
¼ 4x2c þ 4f2 ixc þ p22 þ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ;
ð31Þ
¼ p21 þ 32 Fhc ð1 ce Þ; ¼ 32 Fhc ð1 ce Þ; ¼ 4d30 Fhc ð1 ce Þ;
¼ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; ¼ 38 Fhc ð1 ce e2ixc s Þ; ¼ 4x2c 4f2 ixc þ p22 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; ¼ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; ¼ 38 Fhc ð1 ce e2ixc s Þ; ¼ p22 þ 32 Fhc ð1 ce Þ;
ð32Þ
¼ 32 Fhc ð1 ce Þ; ¼ 4d30 Fhc ð1 ce Þ:
Introducing the polar transformations and the derivative term
A1 ¼
1 id ae ; 2
A01 ¼
1 0 id 1 a e þ i ad0 eid : 2 2
ð33Þ
Moreover, it becomes simpler if we assume
A2 ¼ h b A1 ; where the coefficient hb can be estimated by Eq. (13) in a linear case [28].
ð34Þ
3176
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178 0.012
τ = 0.2160 0.01
a [m]
0.008
τ = 0.2033
0.006
τ = 0.2116 0.004
0.002
0 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
s [mm] Fig. 3. The bifurcation diagram when the s = 0.2160 s, 0.2116 s and 0.2033 s.
Now we can substitute the solutions (26), (28) and (29) into Eq. (24), eliminating the terms that lead to secular terms yields
1 1 ð2ix þ 2f1 ÞA01 ¼ F hc ð1 ce eixc s Þð/3 h3 þ /1 h1 Þð1 hb Þ3 a3 eid 4d0 8 þhc
1 16d20
ð1 ce e
ixc s
1 3 1 Þ 3ð1 hb Þ a3 eid h1 ð1 ce eixc s Þ ð1 hb Þ aeid 8 2 2
)
3
1 1 1 cn hc ð1 hb Þ ð a0 þ i ad0 Þeid cn hc ixc ð1 hb Þ aeid : 2 2 2
ð35Þ
Substituting Eq. (33) into Eq. (35) and separating the real and imaginary parts, we obtain the normal form
a0 ¼ g 1 h1 a þ g 2 a3 ; d0 ¼ g 3 h1 þ g 4 a2 ;
ð36Þ ð37Þ
where
g 3 F 34 ð1 hb Þ ; f1 þ 12 cn hc ð1 hb Þ g þ a3r ð1 ce cos xc sÞ a3i ce sin xc s g2 ¼ 4 ; f1 þ 12 cn hc ð1 hb Þ xc F 34 ð1 hb Þ ðf1 þ 12 cn hc ð1 hb ÞÞ 12 cn h1 xc ð1 hb Þ
g3 ¼ 2 ; x xc a3i ce sin xc s þ f1 þ 12 cn hc ð1 hb Þ 12 cn hc ð1 hb Þ c
g 4 ¼ xc a3i ce sin xc s xc a3r ð1 ce cos xc sÞ þ f1 þ 12 cn hc ð1 hb Þ 2 ða3r ce sin xc s þ a3i ð1 ce cos xc sÞÞg= xc xc a3i ce sin xc s;
þ f1 þ 12 cn hc ð1 hb Þ 12 cn hc ð1 hb Þ ; g1 ¼
ð38Þ
0.012
τ = 0.2033
0.01
a [m]
0.008
τ = 0.1920 0.006
0.004
τ = 0.1994
0.002
0 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
s [mm] Fig. 4. The bifurcation diagram when the s = 0.2033 s, 0.1994 s and 0.1920 s.
3177
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
and the complex parameter a3 = a3r + ia3i that is
3
a3 ¼ hc Fð1 hb Þ3
2
128d0
! 1 ð/3 h3 þ /1 h1 Þ : 32d0
ð39Þ
The approximate solutions of the system are expressed by Eqs. (36) and (37). The qualitative behavior near the Hopf bifurcation can be determined by the sign of g2. For the following result g2> 0, it is subcritical [29]. The amplitude of the steady-state limit cycles can be obtained from the Eq. (36) by setting the derivative term a0 = 0. It yields the expression of determining the constant amplitude of the limit cycles
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g1 ðhch Þ: g2
a¼
ð40Þ
4. Numerical results and discussion Using the MATLABÒ SIMULINK environment the simulation model of the nip system can be set up for the numerical simulations. Several cases have been tested with different rolling speeds. The bifurcation diagram presenting the amplitude of the limit cycle vs the bifurcation parameter from the Eq. (40) can be depicted in Figs. 3. As mentioned before, the layer thickness is s = 1/h. Here the value sc = 1 mm is with respect to hc on the Hopf bifurcation. In Fig. 3, the values of the time delays are s = 0.2160 s (rolling speed 480 m/min) and 0.2116 s (490 m/min), 0.2033 s (510 m/min), respectively. When the rolling speed is 490 m/min, the amplitude is lower than the two others owing to the characteristics of the stability lobe of time delay
-3
-1
x 10
20
40
60
80
0 x1 [m]
x1 [m]
0
0
-3
-3
x 10
-0 .5
-0 .5
-1
100
20
40
20
40
60
80
1 0 -1
100
4
0
20
40
t [s]
60
80
100
3
0
0 0. 005 a [m]
0. 01
0. 015
x 10
4
1 0 -1 -0.01
0.02
-0.005
0
0. 005 a [m]
0. 01
x1 [m]
x1 [m]
0
-1 0
20
40
60
0. 015
80
-1 0
20
40
4 Force [N]
Force [N]
60
80
100
60
80
100
0.015
0.02
t [s]
x 10
1 0 60
80
x 10
2 0 -2 0
100
20
40 t [s]
t [s] 5
5
x 10
Force [N]
4
2 1 0 -1 -0.01
0. 01
-0.5
100
2
3
0
5
40
-0.005
0
0.005 a [m]
0.01
100
0
x 10
5
20
80
2
-2 -0 .01
0.02
t [s]
-1 0
60
x 10
-3
x 10
-0.5
3
40
a [m]
-3
0
Force [N]
0
20
t [s ]
2
Force [N]
1
100
5
Force [N]
x 10
80
0
t [s]
2
60
x 10
5
-0.005
40
2
-2
5
-1 -0.01
20 5
x 10
Force [N]
0
3
0
t [s ]
2
Force [N]
Force [N]
3
1
0
100
5
x 10
2
-1
80
t [s]
5
3
60
x 10
-0 .5 -1
0
t [s]
Force [N]
x1 [m]
0
0.015
0.02
x 10
2 0 -2 -0.01
-0.005
0
0.005 a [m]
0.01
Fig. 5. Time history responses at the rolling speeds 480, 490, 510 (above), 520 and 540 m/min (below).
0.02
3178
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
system [21]. Similarly in Fig. 4, the values of the time delays are s = 0.2033 s (510 m/min) and 0.1994 s (520 m/min), 0.1920 s (540 m/min), respectively. And the amplitude in 520 m/min case is lower than the two others as well. Thus in the rolling speed range of from 480 to 520 m/min, there are jumps from the stable to the unstable alternatively. Let’s check the stability behaviors from the time domain responses next. The numerical results of Eqs. (12) and (13) are illustrated in Fig. 5. Here, we only plot the displacement response x1 of the polymer covered roll since x2 has the same value but an opposite phase of x1. Totally, there are still five cases in different speeds (480, 490, 510, 520, 540 m/min). Each picture includes three subplots as the displacement vs time, the normal force vs. time and the normal force vs the contact width. The nominal contact force is about 60 kN. The displacement amplitudes at the rolling speeds 490 and 520 m/min are much smaller than the others. The stability characteristics are consistent with the bifurcation plots in Figs. 3 and 4. 5. Conclusions The dynamic behaviors of the contact vibrations of the rolling nip system have been investigated through the two-DOF model considering the nonlinear contact stiffness and the time delay effects. The nonlinear normal form of the equations of motion is obtained by the multi-scale method. Thereby the bifurcation diagrams are plotted. The time history responses are consistent with the bifurcation diagrams. Both at the same rolling speeds indicate the same characteristics of the time delay system jumping from the stable to unstable alternatively. Acknowledgment The research is funded by the Academy of Finland. References [1] Smook GA. Handbook for pulp and paper technologists. Angus Wilde Publications Inc.; 1992. [2] Gonzalez Jose A, Abascal R. Efficient stress evaluation of stationary viscoelastic rolling contact problems using the boundary element method: application to viscoelatic coatings. Eng Anal Bound Elem 2006;30:426–34. [3] Karacay T, Akturk N. Vibrations of a grinding spindle supported by angular contact ball bearings. Proc IMech Part K: J Multi-Body Dynam 2008; 222. [4] Bapat C, Batra RC. Indentation of a viscoelastic rubber covered roll by a rigid plane surface. Mech Res Commun 1982;9(4):264–72. [5] Naghieh GR, Jin ZM, Rahnejat H. Contact characteristics of viscoelastic bonded layers. Appl Math Model 1998;22:569–81. [6] Elsharkawy Abdallah A. Visco-elastohydrodynamic lubrication of line contacts. Wear 1996;199:45–53. [7] Digiovanni AA, Chan HM, Harmer MP. The use of hertzian contact in determining coating thickness. J Mater Sci Lett 1997;16:363–7. [8] Van Dooren R. An additional contribution on the vibrations of two elastic bodies in rolling contact. J Sound Vib 2001;242(5):923–8. [9] Thompson DJ. The influence of the contact zone on the excitation of wheel/rail noise. J Sound Vib 2003;267:532–5. [10] Kong XA, Wang Q. A boundary element approach for rolling contact of viscoelastic bodies with friction. Comput Struct 1995;54(No.3):405–13. [11] Bapat C, Batra RC. Indentation of a viscoelastic rubber covered roll by a rigid plane surface. Mech Res Commun 1982;9(4):265–72. [12] Ranganath Nayak P. Contact vibrations. J Sound Vib 1972;22(3):297–322. [13] Narayanan S, Sekar P. A frequency domain based numeric-analytical method for non-linear dynamical systems. J Sound Vib 1998;21(3):409–24. [14] Gray GG, Johnson KL. The dynamic response of elastic bodies in rolling contact to random roughness of their surfaces. J Sound Vib 1972;22(3):323–42. [15] Turner JA. Non-linear vibrations of a beam with Cantilever–Hertzian contact boundary conditions. J Sound Vib 2004;275:177–91. [16] Zaki K, Noah S, Rajagopal KR, Srinivasa AR. Effect of nonlinear stiffness on the motion of a flexible pendulum. Nonlinear Dynam 2002;27:1–18. [17] Oh J-E, Joe Y-G, Shin K. Analysis of out-of-plane motion of a disc brake system using a two-degree-of-freedom model with contact stiffness. Proc IMechE Part D: J Automob Eng 2005;219. [18] Yuan L, Järvenpää VM. Perturbation analysis of nip contact delay system. In: International conference on nonlinear science and complexity, Beijing, China; 2006. [19] Järvenpää VM, Järvinen V, Salmenperä P, Yuan L. Experimental analysis of non-linear roll contact. In: IMAC-XXV conference. Society for Experimental Mechanics, Inc., FL, USA; February 15–18, 2007. [20] Moon FC. Dynamics and chaos in manufacturing processes. John Wiley & Sons, Inc.; 1998. [21] Stepan G. Retarded dynamical systems: stability and characteristic functions. John Wiley & Sons, Inc.; 1989. [22] Kalmar-Nagy T, Stepen G, Moon FC. Subcritial hopf bifurcation in the delay equation model for machine tool vibrations. Nonlinear Dynam 2001;26:121–42. [23] Balachandran B. Nonlinear dynamics of milling processes. Philos Trans R Soc Lond A 2001;359:793–819. [24] Stone E, Campbell SA. Stability and bifurcation analysis of a nonlinear DDE model for drilling. J Nonlinear Sci 2004;14:27–57. [25] Polach O. On non-linear methods of bogie stability assessment using computer simulations. Proc IMechE Part F: J Rail Rapid Transit 2006;220. [26] Nayfeh AH. Problems in perturbation. John Wiley & Sons, Inc.; 1993. [27] Johnson KL. Contact mechanics. Cambridge University Press; 2003. [28] Pratt JR, Nayfeh AH. Chatter control and stability analysis of a cantilever boring bar under regenerative cutting conditions. Philos Trans R Soc Lond A 2001;359:759–92. [29] Thomsen JJ. Vibrations and stability: advanced theory, analysis and tools. Springer; 2003.
Contents lists available at ScienceDirect
Commun Nonlinear Sci Numer Simulat journal homepage: www.elsevier.com/locate/cnsns
Nonlinear vibrations in a covered roll system with viscoelastic contact L. Yuan *, V.-M. Järvenpää Department of Mechanics and Design, Tampere University of Technology, Korkeakoulunkatu 6, P.O. Box 589, 33101 Tampere, Finland
a r t i c l e
i n f o
Article history: Received 1 October 2008 Received in revised form 25 October 2008 Accepted 25 October 2008 Available online 14 November 2008 PACS: 05.45.a Keywords: Viscoelastic cover Time delay Contact vibrations Nonlinear dynamics Perturbation analysis
a b s t r a c t A theoretical analysis on the contact vibration problem of a nip rolling system composed with two paper machine rolls is presented. It is modelled at the contact region by a twodegree-of-freedom (DOF) mass-spring-damper system. The rolling contact force in a case of viscoelastic polymer covered roll is formulated by integration of the contact stress. Regenerative chatter effects to the stability of the vibrations are considered. The perturbation technique is applied through the multi-scale method to calculate the nonlinear normal form of the governing equations to determine the stability behavior of the system. The numerical results in time domain are consistent with the bifurcation diagrams. Ó 2008 Elsevier B.V. All rights reserved.
1. Introduction The contact vibration problem of a roll nip system is presented. Roll nips are used for finishing the surface of processing paper and extensively exist in large scale paper machines of modern paper manufacturing industry. One nip unit consists of two rolls rotating with contact. The traditional hard-nip machine calenders, in which a paper web is calendered to a uniform thickness, cause the variations in paper surface properties [1]. In supercalenders, the nip is typically formed by one hard metal roll and one soft roll coated with polymer layer, so-called soft-nip. In this paper, the soft-nip will be studied. The soft cover may distribute the contact pressure in the nip and produce a paper with more constant density rather than constant caliper of the paper. However, the polymer material has such strong viscoelastic characters that can considerably affect the contact stresses. Consequently a regenerative vibration condition evidently originates from the polymer cover due to its finite viscoelastic recovery time. How to deeply understand the vibration mechanics of the nip contact and substantially mitigate the harmful vibrations (including noises accordingly) are still challenging tasks in the field of the paper manufacturing engineering. In the related literatures, many of researches are concentrated on the contact stress between elastic or viscoelatic bodies in different methods [2–11]. These are mostly for static analysis. On the other hand, some interesting works on dynamic modeling of rolling contact intrigue us. For example, Nayak [12] in earlier 1972 presented a rolling contact dynamic model with two elastic bodies solving by harmonic balance method. Much later in 1997, using the same model as Nayak, Narayana and Sekar [13] solved its equation by Fourier–Galerkin–Newton method and obtained the bifurcation structure of contact vibrations. Gray and Johnson [14] investigated the response of the rolling elastic bodies through the theory and techniques of random vibration. They claimed that contact load and damping are significant parameters. Turner [15] studied nonlinear
* Corresponding author. E-mail address: lihong.yuan@tut.fi (L. Yuan). 1007-5704/$ - see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.cnsns.2008.10.028
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
3171
vibrations with Cantilever–Hertzian contact boundary conditions in a beam application. In recent researches, the dynamic contact force has been analyzed and might be nonlinear from the contact stiffness [16,17]. But it is still needed to discovery more for the aspect that the contact force can be applied into a dynamic model. Such kind of problems belongs to contact vibrations. Authors have studied on the time delay dynamic system [18] whereas some tests and observations on a nip installation have been done and analyzed [19]. Furthermore, to model and predict the vibration motion of the nip system, formulating of the contact force based on the contact theory is needed to explore the soft-nip contact mechanics. The intent of this work is to develop the model of the contact vibrations of the nip system. Firstly, we will present the contact mechanics involved in the deformation of the corrugations on the polymer covered roll surface. The normal contact force is integrated from the pressure distribution between the contact areas. The nonlinear relation between the normal force and the penetration is examined. The variation of the contact force causes the vibrations and appears as a nonlinear term in the mathematical equations of the model. Secondly, we will establish the equations of motion for the dynamic contact vibrations essentially nonlinear, which is of also theoretical interest. Moreover, the time delay term owing to the viscoelastic polymer cover is considered and makes the entire system under a regenerative vibration condition. These two reasons might lead to severe difficulties in the application of approximate techniques for the analysis of these nonlinear time delay differential equations including the solutions and stability. As refer to the stability of time delay system is more complex and interesting especially when it involves nonlinear elements. The analysis of the Hopf bifurcation in different models can be carried out by using analytical methods like multi-scale, harmonic and Floquet theory [20–25]. From the perturbation theory [26], the multi-scale method in this work will be executed for analyzing the mathematical model to obtain the normal form and the amplitude of the limit cycle of the motion. 2. Mathematic modeling of nip contact For this model, there are several assumptions. The time delay is constant, there is no loss of contact in the nip, the small torsional vibration is neglected, the bearing effects are not considered and the temperature generated from the rolling contact is not included. The metal rolls are big (the contact length is 4.4 m and roll radius 2.75 m) and the rolling speed is high (500 m/min). To establish the dynamic vibration model of the nip system, the normal contact force between the rigid roll and the polymer covered roll should be formulated at first. The distribution function p(z) of the normal contact pressure between the two rolls in Fig. 1 is unsymmetrical due to the viscoelastic material behavior on layered contact [27]
pðzÞ ¼
Ea2 1 z2 z þ bnð1 þ nÞð1 eð1þz=aÞ=n Þ : 1 2 bn 1 þ a Rs 2 a
ð1Þ
The static normal force over the contact area can be computed by
P¼
Z
b
pðzÞdz:
ð2Þ
a
We get
" # 3 2 Ea3 1 b b 1 1 2 2 2 bþa an ; P¼ þ a þ abn b 2 bn þ bbn þ abn ð1 þ nÞe n 3 2 Rs 2 6a 2a
ð3Þ
where the equivalent radius R = R1R2/(R1 + R2), Deborah number n = vwebsrelaxation/a, and a R. The dimension coefficients R1 and R2 are radiuses of the two rolls, the sum of the contact widths (a + b) is in the z-direction of the loaded area, s is the thickness of the thin polymer cover layer of the upper roll and E is equivalent elastic modulus [27]. vweb is the speed of the paper roll 1, metal
s polymer cover
R1 p(z)
δ R2
x a
b
roll 2, metal z
Fig. 1. The pressure distribution in the contact region.
3172
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
web through the nip, which corresponds to the running speed of the rolls; srelaxation is the viscoelastic relaxation time of the polymer material and b is the viscoelastic material parameter. The compressive strain ec in an element of the contact at horizontal z-direction is given by
ec ¼ ðd z2 =2RÞ=s; a 6 z 6 b;
ð4Þ
where d is the maximum depth of the penetration of the covered roll. By assuming b a in a symmetric contact, which is the typical case in practical applications
P¼
E 2 2 2 þ bn2 1 þ en þ bn3 1 þ en a3 Rs 3
ð5Þ
or one can rewrite like
P ¼ hF U a3 ;
ð6Þ
where the parameter
FU ¼
E 2 2 2 þ bn2 1 þ en þ bn3 1 þ en R 3
ð7Þ
And h is the reciprocal of the thickness of the polymer cover s, h = 1/s. Fig. 2 illustrates the equivalent spring-mass system for the nip contact rolls. Let x1s = x1(t s). Symbol s is the delay time relevant to the rotational speed of the polymer covered roll. Now we need to seek the relation between the contact width a and the dynamic displacements. In the formula Eq. (4), the strain is zero at z=a; thus one can get
a2 ¼ 2R d:
ð8Þ
The maximum depth of the penetration of the roll d in the dynamical model should be related to the displacements x1, x2 and expressed like
d ¼ d0 þ x1 x2 þ DðtÞ;
ð9Þ
where d0 is the nominal penetration in the contact region and the external excitation is D(t) = Dsin(Xt) from the paper web, which is assumed as zero here.
c1
k1 x1 roll 1 m1
δ0
Δ (t)
Pdynamic
c nip
x2
m2 roll 2 c2
k2
Fig. 2. The equivalent spring-mass system.
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
3173
Then the width of the contact a can be solved as the function of x1 and x2 as
a2 ðtÞ ¼ 2Rðd0 þ x1 x2 Þ;
ð10Þ
and the formula Eq. (10) can be transferred to
aðtÞ ¼
ffi pffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x1 x2 2Rd0 1 þ : d0 d0
ð11Þ
Now the normal force P in formula Eq. (6) becomes the function of x1 and x1s due to the formula Eq. (10), thereby the governing equations of the motion of the rolling system are expressed as
€x1 þ 2f1 x_ 1 þ p21 x1 ¼ Nðx1 ; x1s ; x2 ; x2s ; tÞ=m1 ; €x2 þ 2f2 x_ 2 þ p22 x2 ¼ Nðx1 ; x1s ; x2 ; x2s ; tÞ=m2 ;
ð12Þ ð13Þ
where
2fi ¼ ci =mi ;
p2i ¼ ki =mi ;
and i ¼ 1; 2:
ð14Þ
Considering the time delay effect in formula Eq. (6), the entire dynamic contact force N including the nip damping in the contact is
_ _ Nðx1 ; x1s ; x2 ; x2s ; tÞ ¼ Pdynamic þ h cnip dðtÞ ¼ h F U ½a3 ðtÞ ce a3 ðt sÞ þ h cnip dðtÞ:
ð15Þ
The relaxation parameter ce is related to the time delay and the polymer material characteristics. 3. Perturbation analysis To do the perturbation analysis, the nonlinear cubic term a3(t) obtained from Eqs. (10) and (11) can be expanded as a Taylor series
ffi pffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi x1 x2 2Rd0 1 þ d0 d0 " 2 3 # pffiffiffiffiffiffiffiffiffiffiffi 1 x1 x2 1 x1 x2 1 x1 x2 ; þ ffi 2Rðd0 þ x1 x2 Þ 2Rd0 1 þ 2 d0 d0 8 d0 d0 16 d0 d0
a3 ðtÞ ¼ 2Rðd0 þ x1 x2 Þ
ð16Þ
where higher orders are omitted and xd10 xd20 6 1 must be satisfied. Now the equations of the model Eqs. (12) and (13) evidently consist of nonlinear terms due to the formula Eq. (16). In this system, the masses of the two rolls are close to equal. For simplicity, let m1 = m2. It is assumed that the term FU is constant. Then we can give the notes like
pffiffiffiffiffiffiffiffiffiffiffi F ¼ F 1 ¼ F 2 ¼ ðF U =m1 Þ 2R 2Rd0 ;
cn ¼ cn1 ¼ cn2 ¼ cnip =m1 :
ð17Þ
The perturbation analysis can be performed into the delay type equations using the method of multiple scales [28]. In this method, a fast time-scale T0 = t and a slow time-scale T1 = e2t are introduced. We seek a second-order expansion for the Eqs. (12 and 13) in the form
xi ¼ exi0 ðT 0 ; T 1 Þ þ e2 xi1 ðT 0 ; T 1 Þ þ e3 xi2 ðT 0 ; T 1 Þ þ Oðe4 Þ;
i ¼ 1; 2:
ð18Þ
The time derivatives are
@xi0 @xi1 @xi2 @xi0 þ e2 þ e3 þ 2e3 þ ; @T 0 @T 0 @T 0 @T 1 @ 2 xi0 @ 2 xi1 @ 2 xi2 @ 2 xi0 €xi ¼ e þ e2 þ e3 þ 2e3 þ 2 2 2 @T 0 @T 1 @T 0 @T 0 @T 0 x_ i ¼ e
ð19Þ ð20Þ
Then we perturb the contact parameter (relevant to the effective thickness of the polymer layer), let
h ¼ hc þ e2 h1 ;
ð21Þ
where hc is the reciprocal of thickness value s at the stability boundary. Substituting the Eqs. (18)–(21) into the governing Eqs. (12) and (13), and equating the coefficients of like power of e, we obtain results as following. Order e1: for the first governing Eq. (12)
@ 2 x10 @T 20
þ 2f1
@x10 1 1 þ p21 x10 ¼ F hc ðx10 x20 Þ þ hc ðx10 x20 Þ hc ce ðx10s x20s Þ hc ce ðx10s x20s Þ 2 2 @T 0 @x10 @x20 hc c n ; @T 0 @T 0
ð22aÞ
3174
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
for the second Eq. (13)
@ 2 x20 @T 20
þ 2f2
@x20 1 1 þ p22 x20 ¼ F hc ðx10 x20 Þ þ hc ðx10 x20 Þ hc ce ðx10s x20s Þ hc ce ðx10s x20s Þ 2 2 @T 0 @x10 @x20 þ hc c n : @T 0 @T 0
ð22bÞ
Order e2:
@ 2 x11 @T 20
@ 2 x21 @T 20
þ 2f1
þ 2f2
@x11 1 1 1 þ p21 x11 ¼ F hc ðx10 x20 Þ2 þ hc x10 x20 þ ðx10 x20 Þ2 ðx11 x21 Þ þ 2 8d0 2d0 @T 0 1 1 1 ðx10s x20s Þ2 hc ce x11s x21s þ ðx10s x20s Þ2 ðx11s x21s Þ þ hc ce 2 8d0 2d0 @x11 @x21 þh1 d0 ð1 ce Þ hc cn ; @T 0 @T 0
@x21 1 1 1 þ p22 x21 ¼ F hc ðx10 x20 Þ2 þ hc x10 x20 þ ðx10 x20 Þ2 ðx11 x21 Þ þ 2 8d0 2d0 @T 0 1 1 1 hc ce ðx10s x20s Þ2 hc ce x11s x21s þ ðx10s x20s Þ2 ðx11s x21s Þ þ 2 8d0 2d0 @x11 @x21 h1 d0 ð1 ce Þ þ hc cn : @T 0 @T 0
ð23aÞ
ð23bÞ
Order e3:
@ 2 x12 @T 20
þ2
@ 2 x10 @x12 @x10 1 1 þ 2f1 þ ðx10 x11 þ x20 x21 x10 x21 x11 x20 Þ ðx12 x22 Þ þ p21 x12 ¼ F hc 2 4d0 @T 0 @T 1 @T 0 @T 1 ! 1 3 2 2 3 þ ðx 3x x þ 3x x x Þ 10 10 20 20 10 20 16d20 1 þ hc x12 x22 ðx10 x11 þ x20 x21 x10 x21 x11 x20 Þ d0 ! 1 3 2 2 3 2 ðx10 3x10 x20 þ 3x20 x10 x20 Þ 8d0 1 1 ðx10s x11s þ x20s x21s x10s x21s x11s x20s Þ hc ce ðx12s x22s Þ 2 4d0 ! 1 3 2 2 3 ðx 3x x þ 3x x x Þ þ 10s 10s 20s 20s 10s 20s 16d20 hc ce ðx12s x22s þ
1 ðx10s x11s þ x20s x21s x10s x21s x11s x20s Þ d0
) 2 3 2 2 3 ðx 3x x þ 3x x x ÞÞ þ h x Þð1 c Þ ðx 1 10 20 e 10s 10s 20s 20s 10s 20s 3 8d20 @x10 @x12 @x20 @x22 @x10 @x20 þ þ c n hc þ cn1 hc h1 cn1 : ð24Þ @T 1 @T 0 @T 1 @T 0 @T 0 @T 0
1
The solutions of the Eqs. (22a) and (22b) can be expressed as
x10 ¼ A1 ðT 1 Þeixc T 0 þ A1 ðT 1 Þeixc T 0 ;
x20 ¼ A2 ðT 1 Þeixc T 0 þ A2 ðT 1 Þeixc T 0 ;
x10s ¼ A1 ðT 1 Þeixc ðT 0 sÞ þ A1 ðT 1 Þeixc ðT 0 sÞ ;
x20s ¼ A2 ðT 1 Þeixc ðT 0 sÞ þ A2 ðT 1 Þeixc ðT 0 sÞ ;
ð25Þ
where xc is the chatter frequency at the Hopf bifurcation [20]. Also let
x11 ¼ P1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 ; x21 ¼ R1 e2ixc T 0 þ R2 e2ixc T 0 þ R3 :
ð26Þ
3175
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
By substituting the Eqs. (25) and (26) into the Eq. (23a), we obtain
4x2c ðP1 e2ixc T 0 þ P2 e2ixc T 0 Þ þ 4f1 ixc ðP1 e2ixc T 0 P2 e2ixc T 0 Þ þ p21 ðP1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 Þ
1 1 ¼ F hc ðA1 eixc T 0 þ A1 eixc T 0 A2 eixc T 0 A2 eixc T 0 Þ2 ðP 1 e2ixc T 0 þ P2 e2ixc T 0 þ P3 R1 e2ixc T 0 R2 e2ixc T 0 R3 Þ 2 8d0 1 2ixc T 0 2ixc T 0 2ixc T 0 2ixc T 0 þ P2 e þ P 3 R1 e R2 e R3 Þ ðA1 eixc T 0 þ A1 eixc T 0 A2 eixc T 0 A2 eixc T 0 Þ2 Þ hc ððP1 e 2d0 1 1 þ hc ce ðP1 e2ixc ðT 0 sÞ þ P2 e2ixc ðT 0 sÞ þ P3 R1 e2ixc ðT 0 sÞ R2 e2ixc ðT 0 sÞ R3 Þ ðA1 eixc ðT 0 sÞ þ A1 eixc ðT 0 sÞ 2 8d0 A2 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ Þ2 þ hc ce ððP 1 e2ixc ðT 0 sÞ þ P2 e2ixc ðT 0 sÞ þ P3 R1 e2ixc ðT 0 sÞ R2 e2ixc ðT 0 sÞ R3 Þ 1 þ ðA1 eixc ðT 0 sÞ þ A1 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ A2 eixc ðT 0 sÞ Þ2 Þ h1 d0 þ h1 ce d0 hc cn 2ixc ðP1 e2ixc T 0 P2 e2ixc T 0 Þ 2d0 þ hc cn 2ixc ðR1 e2ixc T 0 R2 e2ixc T 0 Þ:
ð27Þ
For Eq. (27), equating each of the coefficients of exp(2ixcT0), exp(2ixcT0) and exp(0) on both sides, we obtain
P1 ¼ /1 ðA1 A2 Þ2 ; P 2 ¼ /2 ðA1 A2 Þ2 ;
P3 ¼ /3 ðA1 A2 ÞðA1 A2 Þ
ð28Þ
and
R1 ¼ h1 ðA1 A2 Þ2 ;
R2 ¼ h2 ðA1 A2 Þ2 ;
R3 ¼ h3 ðA1 A2 ÞðA1 A2 Þ;
ð29Þ
where
D2 C 3 þ D3 C 1 ; D1 C 1 D2 C 2 D9 C 8 þ D7 C 9 h3 ¼ ; D7 C 7 D8 C 8
/1 ¼
/2 ¼
D5 C6 þ D6 C4 ; D4 C4 D5 C5
/3 ¼
D8 C 9 þ D9 C 7 ; D7 C 7 D8 C 8
h1 ¼
D3 C 2 þ D1 C 3 ; D1 C 1 D2 C 2
h2 ¼
D6 C5 þ D4 C6 ; D4 C4 D5 C5 ð30Þ
and the groups of the parameters
D1 ¼ 4x2c þ 4f1 ixc þ p21 þ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; D2 ¼ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; D3 ¼ 38 Fhc ð1 ce e2ixc s Þ; D4 ¼ 4x2c 4f1 ixc þ p21 þ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; D5 ¼ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; D6 D7 D8 D9
¼ 38 Fhc ð1 ce e2ixc s Þ;
C1 C2 C3 C4 C5 C6 C7 C8 C9
¼ 4x2c þ 4f2 ixc þ p22 þ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ;
ð31Þ
¼ p21 þ 32 Fhc ð1 ce Þ; ¼ 32 Fhc ð1 ce Þ; ¼ 4d30 Fhc ð1 ce Þ;
¼ 32 Fhc ð1 ce e2ixc s Þ þ 2ixc hc cn ; ¼ 38 Fhc ð1 ce e2ixc s Þ; ¼ 4x2c 4f2 ixc þ p22 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; ¼ 32 Fhc ð1 ce e2ixc s Þ 2ixc hc cn ; ¼ 38 Fhc ð1 ce e2ixc s Þ; ¼ p22 þ 32 Fhc ð1 ce Þ;
ð32Þ
¼ 32 Fhc ð1 ce Þ; ¼ 4d30 Fhc ð1 ce Þ:
Introducing the polar transformations and the derivative term
A1 ¼
1 id ae ; 2
A01 ¼
1 0 id 1 a e þ i ad0 eid : 2 2
ð33Þ
Moreover, it becomes simpler if we assume
A2 ¼ h b A1 ; where the coefficient hb can be estimated by Eq. (13) in a linear case [28].
ð34Þ
3176
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178 0.012
τ = 0.2160 0.01
a [m]
0.008
τ = 0.2033
0.006
τ = 0.2116 0.004
0.002
0 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
s [mm] Fig. 3. The bifurcation diagram when the s = 0.2160 s, 0.2116 s and 0.2033 s.
Now we can substitute the solutions (26), (28) and (29) into Eq. (24), eliminating the terms that lead to secular terms yields
1 1 ð2ix þ 2f1 ÞA01 ¼ F hc ð1 ce eixc s Þð/3 h3 þ /1 h1 Þð1 hb Þ3 a3 eid 4d0 8 þhc
1 16d20
ð1 ce e
ixc s
1 3 1 Þ 3ð1 hb Þ a3 eid h1 ð1 ce eixc s Þ ð1 hb Þ aeid 8 2 2
)
3
1 1 1 cn hc ð1 hb Þ ð a0 þ i ad0 Þeid cn hc ixc ð1 hb Þ aeid : 2 2 2
ð35Þ
Substituting Eq. (33) into Eq. (35) and separating the real and imaginary parts, we obtain the normal form
a0 ¼ g 1 h1 a þ g 2 a3 ; d0 ¼ g 3 h1 þ g 4 a2 ;
ð36Þ ð37Þ
where
g 3 F 34 ð1 hb Þ ; f1 þ 12 cn hc ð1 hb Þ g þ a3r ð1 ce cos xc sÞ a3i ce sin xc s g2 ¼ 4 ; f1 þ 12 cn hc ð1 hb Þ xc F 34 ð1 hb Þ ðf1 þ 12 cn hc ð1 hb ÞÞ 12 cn h1 xc ð1 hb Þ
g3 ¼ 2 ; x xc a3i ce sin xc s þ f1 þ 12 cn hc ð1 hb Þ 12 cn hc ð1 hb Þ c
g 4 ¼ xc a3i ce sin xc s xc a3r ð1 ce cos xc sÞ þ f1 þ 12 cn hc ð1 hb Þ 2 ða3r ce sin xc s þ a3i ð1 ce cos xc sÞÞg= xc xc a3i ce sin xc s;
þ f1 þ 12 cn hc ð1 hb Þ 12 cn hc ð1 hb Þ ; g1 ¼
ð38Þ
0.012
τ = 0.2033
0.01
a [m]
0.008
τ = 0.1920 0.006
0.004
τ = 0.1994
0.002
0 0.7
0.8
0.9
1
1.1
1.2
1.3
1.4
s [mm] Fig. 4. The bifurcation diagram when the s = 0.2033 s, 0.1994 s and 0.1920 s.
3177
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
and the complex parameter a3 = a3r + ia3i that is
3
a3 ¼ hc Fð1 hb Þ3
2
128d0
! 1 ð/3 h3 þ /1 h1 Þ : 32d0
ð39Þ
The approximate solutions of the system are expressed by Eqs. (36) and (37). The qualitative behavior near the Hopf bifurcation can be determined by the sign of g2. For the following result g2> 0, it is subcritical [29]. The amplitude of the steady-state limit cycles can be obtained from the Eq. (36) by setting the derivative term a0 = 0. It yields the expression of determining the constant amplitude of the limit cycles
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi g1 ðhch Þ: g2
a¼
ð40Þ
4. Numerical results and discussion Using the MATLABÒ SIMULINK environment the simulation model of the nip system can be set up for the numerical simulations. Several cases have been tested with different rolling speeds. The bifurcation diagram presenting the amplitude of the limit cycle vs the bifurcation parameter from the Eq. (40) can be depicted in Figs. 3. As mentioned before, the layer thickness is s = 1/h. Here the value sc = 1 mm is with respect to hc on the Hopf bifurcation. In Fig. 3, the values of the time delays are s = 0.2160 s (rolling speed 480 m/min) and 0.2116 s (490 m/min), 0.2033 s (510 m/min), respectively. When the rolling speed is 490 m/min, the amplitude is lower than the two others owing to the characteristics of the stability lobe of time delay
-3
-1
x 10
20
40
60
80
0 x1 [m]
x1 [m]
0
0
-3
-3
x 10
-0 .5
-0 .5
-1
100
20
40
20
40
60
80
1 0 -1
100
4
0
20
40
t [s]
60
80
100
3
0
0 0. 005 a [m]
0. 01
0. 015
x 10
4
1 0 -1 -0.01
0.02
-0.005
0
0. 005 a [m]
0. 01
x1 [m]
x1 [m]
0
-1 0
20
40
60
0. 015
80
-1 0
20
40
4 Force [N]
Force [N]
60
80
100
60
80
100
0.015
0.02
t [s]
x 10
1 0 60
80
x 10
2 0 -2 0
100
20
40 t [s]
t [s] 5
5
x 10
Force [N]
4
2 1 0 -1 -0.01
0. 01
-0.5
100
2
3
0
5
40
-0.005
0
0.005 a [m]
0.01
100
0
x 10
5
20
80
2
-2 -0 .01
0.02
t [s]
-1 0
60
x 10
-3
x 10
-0.5
3
40
a [m]
-3
0
Force [N]
0
20
t [s ]
2
Force [N]
1
100
5
Force [N]
x 10
80
0
t [s]
2
60
x 10
5
-0.005
40
2
-2
5
-1 -0.01
20 5
x 10
Force [N]
0
3
0
t [s ]
2
Force [N]
Force [N]
3
1
0
100
5
x 10
2
-1
80
t [s]
5
3
60
x 10
-0 .5 -1
0
t [s]
Force [N]
x1 [m]
0
0.015
0.02
x 10
2 0 -2 -0.01
-0.005
0
0.005 a [m]
0.01
Fig. 5. Time history responses at the rolling speeds 480, 490, 510 (above), 520 and 540 m/min (below).
0.02
3178
L. Yuan, V.-M. Järvenpää / Commun Nonlinear Sci Numer Simulat 14 (2009) 3170–3178
system [21]. Similarly in Fig. 4, the values of the time delays are s = 0.2033 s (510 m/min) and 0.1994 s (520 m/min), 0.1920 s (540 m/min), respectively. And the amplitude in 520 m/min case is lower than the two others as well. Thus in the rolling speed range of from 480 to 520 m/min, there are jumps from the stable to the unstable alternatively. Let’s check the stability behaviors from the time domain responses next. The numerical results of Eqs. (12) and (13) are illustrated in Fig. 5. Here, we only plot the displacement response x1 of the polymer covered roll since x2 has the same value but an opposite phase of x1. Totally, there are still five cases in different speeds (480, 490, 510, 520, 540 m/min). Each picture includes three subplots as the displacement vs time, the normal force vs. time and the normal force vs the contact width. The nominal contact force is about 60 kN. The displacement amplitudes at the rolling speeds 490 and 520 m/min are much smaller than the others. The stability characteristics are consistent with the bifurcation plots in Figs. 3 and 4. 5. Conclusions The dynamic behaviors of the contact vibrations of the rolling nip system have been investigated through the two-DOF model considering the nonlinear contact stiffness and the time delay effects. The nonlinear normal form of the equations of motion is obtained by the multi-scale method. Thereby the bifurcation diagrams are plotted. The time history responses are consistent with the bifurcation diagrams. Both at the same rolling speeds indicate the same characteristics of the time delay system jumping from the stable to unstable alternatively. Acknowledgment The research is funded by the Academy of Finland. References [1] Smook GA. Handbook for pulp and paper technologists. Angus Wilde Publications Inc.; 1992. [2] Gonzalez Jose A, Abascal R. Efficient stress evaluation of stationary viscoelastic rolling contact problems using the boundary element method: application to viscoelatic coatings. Eng Anal Bound Elem 2006;30:426–34. [3] Karacay T, Akturk N. Vibrations of a grinding spindle supported by angular contact ball bearings. Proc IMech Part K: J Multi-Body Dynam 2008; 222. [4] Bapat C, Batra RC. Indentation of a viscoelastic rubber covered roll by a rigid plane surface. Mech Res Commun 1982;9(4):264–72. [5] Naghieh GR, Jin ZM, Rahnejat H. Contact characteristics of viscoelastic bonded layers. Appl Math Model 1998;22:569–81. [6] Elsharkawy Abdallah A. Visco-elastohydrodynamic lubrication of line contacts. Wear 1996;199:45–53. [7] Digiovanni AA, Chan HM, Harmer MP. The use of hertzian contact in determining coating thickness. J Mater Sci Lett 1997;16:363–7. [8] Van Dooren R. An additional contribution on the vibrations of two elastic bodies in rolling contact. J Sound Vib 2001;242(5):923–8. [9] Thompson DJ. The influence of the contact zone on the excitation of wheel/rail noise. J Sound Vib 2003;267:532–5. [10] Kong XA, Wang Q. A boundary element approach for rolling contact of viscoelastic bodies with friction. Comput Struct 1995;54(No.3):405–13. [11] Bapat C, Batra RC. Indentation of a viscoelastic rubber covered roll by a rigid plane surface. Mech Res Commun 1982;9(4):265–72. [12] Ranganath Nayak P. Contact vibrations. J Sound Vib 1972;22(3):297–322. [13] Narayanan S, Sekar P. A frequency domain based numeric-analytical method for non-linear dynamical systems. J Sound Vib 1998;21(3):409–24. [14] Gray GG, Johnson KL. The dynamic response of elastic bodies in rolling contact to random roughness of their surfaces. J Sound Vib 1972;22(3):323–42. [15] Turner JA. Non-linear vibrations of a beam with Cantilever–Hertzian contact boundary conditions. J Sound Vib 2004;275:177–91. [16] Zaki K, Noah S, Rajagopal KR, Srinivasa AR. Effect of nonlinear stiffness on the motion of a flexible pendulum. Nonlinear Dynam 2002;27:1–18. [17] Oh J-E, Joe Y-G, Shin K. Analysis of out-of-plane motion of a disc brake system using a two-degree-of-freedom model with contact stiffness. Proc IMechE Part D: J Automob Eng 2005;219. [18] Yuan L, Järvenpää VM. Perturbation analysis of nip contact delay system. In: International conference on nonlinear science and complexity, Beijing, China; 2006. [19] Järvenpää VM, Järvinen V, Salmenperä P, Yuan L. Experimental analysis of non-linear roll contact. In: IMAC-XXV conference. Society for Experimental Mechanics, Inc., FL, USA; February 15–18, 2007. [20] Moon FC. Dynamics and chaos in manufacturing processes. John Wiley & Sons, Inc.; 1998. [21] Stepan G. Retarded dynamical systems: stability and characteristic functions. John Wiley & Sons, Inc.; 1989. [22] Kalmar-Nagy T, Stepen G, Moon FC. Subcritial hopf bifurcation in the delay equation model for machine tool vibrations. Nonlinear Dynam 2001;26:121–42. [23] Balachandran B. Nonlinear dynamics of milling processes. Philos Trans R Soc Lond A 2001;359:793–819. [24] Stone E, Campbell SA. Stability and bifurcation analysis of a nonlinear DDE model for drilling. J Nonlinear Sci 2004;14:27–57. [25] Polach O. On non-linear methods of bogie stability assessment using computer simulations. Proc IMechE Part F: J Rail Rapid Transit 2006;220. [26] Nayfeh AH. Problems in perturbation. John Wiley & Sons, Inc.; 1993. [27] Johnson KL. Contact mechanics. Cambridge University Press; 2003. [28] Pratt JR, Nayfeh AH. Chatter control and stability analysis of a cantilever boring bar under regenerative cutting conditions. Philos Trans R Soc Lond A 2001;359:759–92. [29] Thomsen JJ. Vibrations and stability: advanced theory, analysis and tools. Springer; 2003.