Jtam.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 46, 3, pp. 483-509, Warsaw 2008 DYNAMICS OF MULTI-PENDULUM SYSTEMS WITH FRACTIONAL ORDER CREEP ELEMENTS Katica R. (Stevanović) Hedrih Faculty od Mechanical Engineering University of Nǐs, Mathematical Institute SANU, Nǐs, Serbia e-mail: katica@masfak.ni.ac.yu; khedrih@eunet.yu A survey as a short review of author’s research results in area of dynamics of hybrid systems and analytical dynamics of discretematerial particle sys- tem containing creep elements described by fractional order derivatives, is presented.Freevibrationsof amulti-pendulumsystem intercoupledby stan- dard light elements and different properties are considered. The correspon- ding system of an ordinary fractional order as well as integro-differential equations, described dynamics of the multi-pendulum system, are derived and analytically solved. For the case of one pendulum and two pendulum systems containing standard light creep elementswith the stress-strain con- stitutive relation expressed by a fractional order derivative, ordinary diffe- rential equations are analytically solved. From the analytical solutions, for the case of the homogeneous two-pendulum system, it is visible that free vibrations under arbitrary initial conditions contain three modes, one pure periodic and two aperiodic expressed by time series expansions. The obta- ined analytical solution modes are numerically analysed and characteristic vibrationmodes for different kinetic parameters are graphically presented. Key words: coupled subsystems, coupled dynamics,multi-pendulum system 1. Introduction Fast development of material science (Rzhanitsin, 1949; Savin andRuschisky, 1976) and experimental mechanics as well as ofmethods of numerical analysis led to the creation of differentmodels of real material bodies (Enelund, 1996) andmethods for studying dynamics and processes which happen in them du- ring the traveling of a disturbance through deformable bodies. The interest in the studyof coupled systems (Hedrih, 1999, 2003a, 2004b,c, 2005a,b, 2006a,b,e,f, 2007a,c) as new qualitative systems has grown exponen- tially over the last few years because of the theoretical challenges involved 484 K.R. Hedrih in the investigation of such systems with hereditary elements (Goroško and Hedrih, 2001, 2007a,b; Goroshko and Puchko, 1977) with stress-strain con- stitutive relations expressed by integral forms as well as with creep elements with the stress-strain constitutive relation described by fractional order de- rivatives(Hedrih, 2002a,b, 2003b, 2004c, 2005c, 2006c; Hedrih and Filipovski, 2002). In the References, monographs by Rzhanitsin (1949) and Savin and Ru- schisky (1976) as well as by Goroško and Hedrih (2001, 2007a,b) and paper byGoroshko andPuchko (1977), different approaches to creation of real body models are given. In basic, these approaches contain physical discretization of a body or mathematical discretisation from partial to ordinary differen- tial equations. One such an approach is represented by a model of a discrete system of material particles (Hedrih, 2001, 2003a, 2004b, 2006a), which are connected by certain ties. The number of the particles then increased to create a continuum, their motion and deformable wave propagation is described by partial differential equations. And then, due to impossibility of solving them analytically, an approximation method is used for that purpose. The book by Enelund (1996) contains the same applications with elements of fractional calculus in structural dynamics. In the monograph by Gorošhko and Hedrih (2001), analytical dynamics of discrete hereditary systems and corresponding solutions was first published as an integral theory such kind of systems. In Goroško and Hedrih (2001, 2007a,b) andGoroshko and Puchko (1977) as well as in the citedmonograph, a standard light hereditary element is defined and used as a connecting or coupling rheological element between the material particles of the system. A series of References byHedrih (2002b, 2006c, 2007b,e), Hedrih and Fili- povski (2002) and [10] present numerous results of research on the properties of vibrations of rods, plates, belts made of differentmaterials. Also, in Hedrih (1999, 2001, 2003a, 2004b, 2005a,b, 2006a,b,e,f, 2007a,c) a series of coupled subsystems and hybrid systems with different material properties or different properties of the standard light elements as discrete as well as distributed coupling elements between deformable bodies, or discrete and continuum sub- systems (seeHedrih, 2005a, 2006b) are investigated. In a series of the author’s work (Hedrih, 2003b, 2005c, 2006c, 2007b,e), results of research on vibrations of deformable bodies with creep properties described by the stress strain con- stitutive relation expressed by a fractional order derivative, are presented. All these engineering problems are, also, mathematical problems and are described by partial differential equations with integral or fractional order derivative termswhich can be discretised into a problemof solving of a system Dynamics of multi-pendulum systems... 485 of ordinary differential or integro-differential or fractional order differential equations. In the last decade, an interest to the applied fractional calculus for de- scription of material properties rose. Papers by Enelund (1996), Gorenflo and Mainardi (2000) are recommended as the primary mathematical literature containing the basics of the fractional calculus. Hedrih (2006a) studied modes of homogeneous chain signals for different kinds of homogenous connections betweenmaterialmass particles in the chain and different chain boundary conditions. A finite number of coupled fractional order differential equations of creep vibrations of connectedmulti-mass partic- les into ahomogeneous chain systemhasbeenderived.Themass particleswere connected by standard creeping light elements and the constitutive relations of the stress-strain state were expressed through terms of the fractional order derivatives. The analytical solution to the system of coupled fractional order differential equations of the corresponding dynamical free creep processes was obtained by using Laplace’s transform method and trigonometrical method (see Rašković, 1965). By using inverse Laplace’s transform, time series func- tions as particular mode components of the solution were found. By using those component visualizations, analysis of the dynamical creep component processes inmass particle displacements were done. Also, analysis and a com- parison between signals in the corresponding homogeneous chains with ideal elastic or visco-elastic standard light elements between the material particles were pointed out. 2. Light standard elements The basic elements of a discrete material system with interconnections be- tween material particles as well as the mathematical multi-pendulum system considered in this paper, are: • Material particle with mass mk having one degree of freedom, defined by the following independent and generalised angular coordinate ϕk, for k=1,2, . . . ,n. • Light standard coupling element of negligible mass in the form of an axially stressed rod without bending, having the ability to resist defor- mation under static and dynamic conditions. The constitutive relation between the restitution force P and elongation x can bewrittendown in the form fpsr(P,Ṗ,x,ẋ,x α t ,D,D α t ,J,n,c, c̃,µ,α,cα,T,U,. . .) = 0, where 486 K.R. Hedrih D,Dαt and J are fractional order differential and integral operators (Go- roško and Hedrih, 2001; Hedrih, 2006a), which find their justification in experimental verifications of material behaviour (Goroško and Hedrih, 2001; Rzhanitsin, 1949; Savin andRuschisky, 1976), while n, c, c̃, µ, cα, α, . . .are material constants, which are also experimentally determined. For every single standard coupling light element of negligible mass, we shall define a specific stress-strain constitutive relation-law of dynamics. This means thatwewill definethe stress-strain constitutive relation asdeterminants of forces and/or changes of forces with distances between two constrained- coupled material particles and with changes of the distances in time, with accuracy up to constants which dependon the accuracy of their determination through an experiment. The accuracy of those constants, forces and elongations would depend not only on the nature of an object, but also on the knowledge of very complex stress-strain relations to be dealt with (seeGoroško andHedrih, 2001; Hedrih, 2006a). In thispaper,weshall use three such light standardconstraint-coupling elements, and they will be: • Light standard ideally elastic coupling element forwhich the stress-strain relation for the restitution force as a function of element axial elongation is given by a linear relation of the form P =−cy (2.1) as well as given by a nonlinear relation of the form P =−cy− ĉy3 (2.2) where c is the rigidity coefficient or elasticity coefficient for the linear, and ĉ= εχc for the nonlinear functional stress-strain constitutive rela- tion between the force and rheological coordinate of axial deformation of the standard elastic element. In a natural state, non-stressed by a force and undeformed, the force and deformation of such an element are equal to zero. • Light standard hereditary constraint element for which the stress-strain constitutive relation for the restitution force as a function of element elongation (rheological coordinate) is given: — in a differential form DP =Cy or nhṖ(t)+P(t)=nhchy(t)+ c̃hy(t) (2.3) Dynamics of multi-pendulum systems... 487 where the following differential operators are introduced D=nh d dt +1 C=nhch d dt + c̃h (2.4) and nh is the relaxation time and ch, c̃h are rigidity coefficients – mo- mentary and prolonged. — in an integral form P(t)= ch [ y(t)− t∫ 0 R(t− τ)y(τ) dτ ] (2.5) where R(t− τ)= ch− c̃h nhch exp [ − 1 nn (t− τ) ] (2.6) is the relaxation kernel (or resolvent). — in an integral form y(t)= 1 ch [ P(t)+ t∫ 0 K(t− τ)P(τ) dτ ] (2.7) where K(t− τ)= ch− c̃h nhch exp [ − c̃h nhch (t− τ) ] (2.8) is kernel of rheology (or retardation). • Light standard creep coupling element for which the stress-strain con- stitutive relation for the restitution force as a function of element elon- gation (rheological coordinate) is given by a fractional order derivative term (see Hedrih, 2006a) in the form P(t)=−{c0x(t)+ cαD α t [x(t)]} (2.9) where Dαt [·] is the fractional order differential operator of the α-th de- rivative with respect to time t in the following form D α t [x(t)] = dαx(t) dtα =x(α)(t)= 1 Γ (1−α) d dt t∫ 0 x(τ) (t− τ)α dτ (2.10) where c0, cα are rigidity coefficients – momentary and prolonged ones, and α is a rational number between 0 and 1, 0<α< 1, depending on the material properties found experimentally. 488 K.R. Hedrih 3. Multi-pendulum system with intercoupling by standard light elements In this paper, we shall define a discrete continuum mathematical pendulum chain as a system of material particles intercoupled by the light standard co- upling elements (elastic, hereditary or creep), presented in Figures 1 and 2, and which are, in the natural state, on the defined interdistances (when the coupling elements are unstressed and without prehistory as well as without memory before the initial moment of the system motion). The discrete continuummulti-pendulum chain system is an ideally elastic chain if its material particles are interconnected by light standard ideally ela- stic coupling elements with the stress-strain constitutive relation expressed by (2.1) for the linear or (2.2) for the nonlinear case. The discrete continuum multi-pendulum chain system is a standard hereditary chain if its material particles are interconnected by light standard hereditary elements with the stress-strain constitutive relation expressed by one of sets (2.3)-(2.8). The di- screte continuum multi-pendulum chain system is a standard creep chain if itsmaterial particles are interconnected by light standard creep elements with the stress-strain constitutive relation expressed by (2.9) and (2.10). Weshall defineadiscrete homogeneousmathematicalmulti-pendulumchain system as a system of discrete material particles, same masses, which can rotate along corresponding circular arches with the same radius ℓ and centers on the one horizontal line. All the system is in the vertical plane and in the gravitational field and thematerial particles are intercoupled by the same type of sets of the parallel standard light elements and on the same distance of the corresponding fixed points of pendula. The number of degrees of freedom of each of thesemulti-pendulum chains is equal to the number ofmaterial particles in it sincewe accept the previously defined character of the system. Further, for a special rheolinear case, we introduce a hypotheses about ho- mogeneity of the discrete continual chain, about small deformations of light standard coupling elements and small displacements of the material particles. Also, we introduce a hypothesis that the homogenous discrete continuum chain, is in natural, non-stressed state, before the initial moment of motion, i.e. that the light standard coupling elements do not have prehistory nor the memory of the stress-strain state. With these hypotheses, we shall direct our research to the dynamics of chain-like homogenous multi-pendulum systems. Dynamics of multi-pendulum systems... 489 4. Thermo-rheological coupled multi-pendulum system In Fig.1, a thermo-rheological system containing a finite number of coupled pendula is presented.We take into consideration the finite number of coupled mathematical pendula presented in Fig.1 as a system with material partic- les mk with the length ℓk and with finite numbers (n) of degrees of freedom definedby thegeneralized coordinates ϕk,k=1,2, . . . ,m, and the standard li- ght thermo-visco-elastic elements thermo-modified by temperature Tk(t), and the coefficient of thermo-dilatation cTk coupling the pendula at the distan- ce ℓTk, parallelly coupled but temperature isolated, and with the standard light nonlinear springs with coefficients of the linear and nonlinear rigidity respectively denoted by ck, and ĉ= εχkck, where ε is a small parameter. Fig. 1. Systemwith four pendula interconnected by standard light thermo-modified hereditary element and nonlinear springs Now, we take into account that these standard light thermo-visco- elastic elements, with natural lengths ℓ0k thermo-modified by the tempera- ture Tk(t) are in a dynamical state, and that we do not neglect the thermo- modification of the element strain.We can write that the thermo-dilatation is ∆ℓ0k = αTkTk(t)(ℓ0k +xk) and that the constitutive relation of the thermo- elastic stress-strain state is expressed in the following form Pher(k)(t)=−cT(k)[∆ℓ0(k)+ ℓT(k)(ϕk+1−ϕk)]= (4.1) =−cT(k)ℓT(k)(ϕk+1−ϕk)[1+αT(k)T(k)(t)]− cT(k)αT(k)ℓ0(k)T(k)(t) and the forces of the nonlinear spring between pendula are in the following form Pnl(k)(t)=−c(k)[ℓc(k)(ϕk+1−ϕk)+εχkℓ 3 c(k)(ϕk+1−ϕk) 3] (4.2) and the forces of damping are Pd(k)(t)=−bkℓb(k)(ϕ̇k+1− ϕ̇k) (4.3) 490 K.R. Hedrih For the light standard creep coupling element between two pendula, the stress- strain relation for the restitution force as a function of the element elongation is given by fractional order derivatives in the form Pk.k+1(t)=−{c0k[ℓc(k)(ϕk+1(t)−ϕk(t))]+ cαkD α t [ℓc(k)(ϕk+1(t)−ϕk(t))]}= (4.4) =−{c0kℓc(k)[ϕk+1(t)−ϕk(t)]+ cαkℓc(k)〈D α t [ϕk+1(t)]−D α t [ϕk(t)]〉} where Dαt [·] is the fractional order differential operator of the α-th derivati- ve with respect to time t in form (2.10), c0k, cαk are rigidity coefficients – momentary and prolonged ones, and αk is a rational number 0<αk < 1. For the light standard hereditary constraint element between two pendula, the stress-strain relation for the restitution force as a function of the element elongation is given by integral term in the form Pher(k.k+1)(t)= chk [ ℓc(k)[ϕk+1(t)−ϕk(t)]− t∫ 0 Rk(t−τ)ℓc(k)[ϕk+1(τ)−ϕk(τ)] dτ ] (4.5) where Rk(t− τ)= chk− c̃hk nhkc exp [ − 1 nhk (t− τ) ] (4.6) is the relaxation kernel (or resolvent), and where chk, c̃hk aremomentary and prolonged rigidity coefficients and nhk is the relaxation time of an element. 5. Governing equations of the multi-pendulum system – general case Now, we take into account that the pendula are intercoupled by parallel co- upled sets of the standard light elements of different properties, as it was presented in the previous part of this paper. Suppose that there are n pendu- la, and in the equilibrium state the pendulumsystem is in the vertical position as presented in Fig.1, i.e. when all generalized coordinates are equal to zero, ϕk,eq =0. The generalized forces corresponding to the generalized coordinates ϕk between two pendula are Qher(k)(t)=Pher(k,k+1)(t)ℓT(k) = =−cT(k)ℓ 2 T(k)(ϕk+1−ϕk)[1+αT(k)T(k)(t)]− cT(k)αT(k)ℓT(k)ℓ0(k)T(k)(t) Dynamics of multi-pendulum systems... 491 Qnl(k1)(t)=Pnl(k,k+1)(t)ℓc(k) =−c(k)ℓc(k)[ℓc(k)(ϕk+1−ϕk)+ +εχkℓ 3 c(k)(ϕk+1−ϕk) 3] Qd(k)(t)=Pd(k,k+1)(t)ℓb(k) =−bkℓ 2 b(k)(ϕ̇k+1− ϕ̇k) (5.1) Qcr(k)(t)=Pcr(k.k+1)(t)ℓ0(k) = =−ℓ20(k){c0k[ϕk+1(t)−ϕk(t)]+cαk〈D α t [ϕk+1(t)]−D α t [ϕk(t)]〉} Qher(k.)(t)=Pher(k.k+1)(t)ℓh(k) = chkℓ 2 h(k) [ [ϕk+1(t)−ϕk(t)]+ − t∫ 0 Rk(t− τ)[ϕk+1(τ)−ϕk(τ)] dτ ] Qg,k = ∂Epϕn ∂ϕk =mkgℓk sinϕk ≈mkgℓk (ϕk 1! − ϕ3k 3! + ϕ5k 5! − ϕ7k 7! + ϕ9k 9! + . . . ) The system of governing differential equations of the thermo-rheological coupled multi-pendulum system presented in Fig.1, is in the following form mkℓ 2 kϕ̈k =−cTk−1ℓ 2 Tk−1(ϕk −ϕk−1)[1+αTk−1Tk−1(t)]+ −cTk−1αTk−1ℓ0l−1ℓTk−1Tk−1(t)− bk−1ℓ 2 bk−1(ϕ̇k − ϕ̇k−1)+ −mkgℓk (ϕk 1! − ϕ3k 3! + ϕ5k 5! − ϕ7k 7! + ϕ9k 9! + . . . ) + −ck−1ℓck−1[ℓk−1(ϕk −ϕk−1)+εχk−1ℓ 3 ck−1(ϕk −ϕk−1) 3]+ +cTkℓ 2 Tk(ϕk+1−ϕk)[1+αTkTk(t)]+cTkαTkℓ0kℓTkTk(t)+ +bkℓ 2 bk(ϕ̇k+1− ϕ̇k)+ckℓck[ℓck(ϕk+1−ϕk)+εχkℓ 3 ck(ϕk+1−ϕk) 3]+ +ℓc(k){c0kℓc(k)[ϕk+1(t)−ϕk(t)]+ cαkℓc(k)〈D α t [ϕk+1(t)]−D α t [ϕk(t)]〉}+ −ℓc(k−1){c0k−1ℓc(k−1)[ϕk(t)−ϕk−1(t)]+ (5.2) +cαk−1ℓc(k−1)〈D α t [ϕk(t)]−D α t [ϕk−1(t)]〉}+chkℓc(k) [ ℓc(k)[ϕk+1(t)−ϕk(t)]+ − t∫ 0 Rk(t− τ)ℓc(k)[ϕk+1(τ)−ϕk(τ)] dτ ] + −chk−1ℓc(k−1) [ ℓc(k−1)[ϕk(t)−ϕk−1(t)]+ − t∫ 0 Rk=1(t− τ)ℓc(k−1)[ϕk(τ)−ϕk−1(τ)] dτ ] where k=1,2, . . . ,n, ϕ0 =0, and ϕn+1 =0. 492 K.R. Hedrih Fig. 2. Systemwith ”chain” pendula interconnected by standard light thermo-modified hereditary elements and creep elements A secial case is a homogeneousmulti-pendulum system shown in Fig.2 for all equal lengths ℓ and the same coupling setswithparallelly coupled standard light elements between the pendula. After introducing the following notations ω20 = c m ω20T = cT m ω̃20 = g ℓ γ =αTT0 2δ= b m h0 = αTℓ0T0 ℓ T̃(t)= 1 T0 T(t) χ̈=χℓ2 ω200 = c0 m ω20α = cα m ω20h = ch m the previous system of equations can be transformed into the following ϕ̈1+ ω̃ 2 0ϕ1−ω 2 0(ϕ2−ϕ1)−ω 2 00(ϕ2−ϕ1)−ω 2 0T(ϕ2−ϕ1)[1+γT̃(t)]+ −2δ(ϕ̇2− ϕ̇1)=ω 2 0Th0T̃(t)− ω̃ 2 0 ( − ϕ31 3! + ϕ51 5! − ϕ71 7! + ϕ91 9! + . . . ) + +εω20χ̃(ϕ2−ϕ1) 3+ω20α{D α t [ϕ2(t)−D α t [ϕ1(t)]}+ +ω20h [ [ϕ2(t)−ϕ1(t)]− t∫ 0 R(t− τ)[ϕ2(τ)−ϕ1(τ)] dτ ] . . . ϕ̈k + ω̃ 2 0ϕk+ω 2 0(ϕk −ϕk−1)+ω 2 00(ϕk −ϕk−1)−ω 2 00(ϕk+1−ϕk)+ −ω20(ϕk+1−ϕk)+ω 2 0T(ϕk −ϕk−1)[1+γT̃(t)]+ −ω20T(ϕk+1−ϕk)[1+γT̃(t)]+2δ(ϕ̇k − ϕ̇k−1)−2δ(ϕ̇k+1− ϕ̇k)= =−ω20Th0T̃(t)+ω 2 0Th0T̃(t)++ω̃ 2 0 (ϕ3k 3! − ϕ5k 5! + ϕ7k 7! − ϕ9k 9! + . . . ) + −εω20χ̃(ϕk −ϕk−1) 3+εω20χ̃(ϕk+1−ϕk) 3−ω20α{D α t [ϕk(t)]−D α t [ϕk−1(t)]}+ Dynamics of multi-pendulum systems... 493 +ω20α{D α t [ϕk+1(t)]−D α t [ϕk(t)]}+ (5.3) +ω20hℓ [ [ϕk+1(t)−ϕk(t)]− t∫ 0 R(t− τ)[ϕk+1(τ)−ϕk(τ)] dτ ] + −ω20h [ [ϕk(t)−ϕk−1(t)]− t∫ 0 R(t− τ)[ϕk(τ)−ϕk−1(τ)] dτ ] . . . ϕ̈n+ ω̃ 2 0ϕn+ω 2 0(ϕn−ϕn−1)+ω 2 0T(ϕn−ϕn−1)[1+γT̃(t)]+ +2δ(ϕ̇n− ϕ̇n−1)=−ω 2 0Th0T̃(t)+ ω̃ 2 0 (ϕ3n 3! − ϕ5n 5! + ϕ7n 7! − ϕ9n 9! + . . . ) + −εω20χ̃(ϕn−ϕn−1) 3−ω20α{D α t [ϕn(t)]−D α t [ϕn−1(t)]}+ −ω20h [ [ϕn(t)−ϕn−1(t)]− t∫ 0 R(t− τ)[ϕn(τ)−ϕm−1(τ)] dτ ] The basic linear ordinary differential equations of the previous system for the homogeneous case are in the following form ϕ̈1+(ω̃ 2 0 +ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T)ϕ1− (ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T)ϕ2 =0 . . . ϕ̈k− (ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T)ϕk−1+(ω̃ 2 0 +2ω 2 0 +2ω 2 00+2ω 2 0h+2ω 2 0T)ϕk+ −(ω20 +ω 2 00+ω 2 0h+ω 2 0T)ϕk+1 =0 (5.4) . . . ϕ̈n− (ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T)ϕ+(ω̃ 2 0 +ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T)ϕn =0 By introducing the following notations ω20cT =(ω 2 0 +ω 2 00+ω 2 0h+ω 2 0T) u= ω̃2− ω̃20 ω̃20cT (5.5) formally for obtaining eigen amplitude vectors of previous system (5.4), it is possible to write matrix equations (see Hedrih, 2004b, 2006a) (C̃−uÃ)A=    1−u −1 0 0 0 −1 2−u −1 0 0 0 −1 2−u −1 0 0 0 −1 2−u −1 0 0 0 −1 1−u       A1 A2 A3 A4 A5    =0 (5.6) 494 K.R. Hedrih and by the use of the trigonometrical method (see Rašković, 1965; Hedrih, 2006a), for a free homogeneous coupled pendulum system, one obtains the eigen numbers in the form us = ω2s − ω̃ 2 0 ω̃20cT =4sin2 sπ 2n s=0,1,2, . . . ,n−1 (5.7) with the eigen frequencies ωs = √ ω̃20 +4(ω 2 0c+ω 2 0T)sin 2 sπ 2n s=0,1,2, . . . ,n−1 (5.8) For a three-pendulum system the eigen-frequencies are ω2s = ω̃ 2 0cTus+ ω̃ 2 0 =    ω̃20 ω̃20cT + ω̃ 2 0 3ω̃20cT + ω̃ 2 0 =    ω̃20 ω20T +ω 2 0 + ω̃ 2 0 +ω 2 00+ω 2 0h 3ω20T +3ω 2 0c+3ω 2 00+3ω 2 0h+ ω̃ 2 0 (5.9) and the solution for a free three-pendulum system is in the form ϕ1(t)=C1cos(ω1t+α1)+C2cos(ω2t+α2)+C3cos(ω3t+α3) ϕ2(t)=C1cos(ω1t+α1)−2C3cos(ω3t+α3) (5.10) ϕ3(t)=C1cos(ω1t+α1)−C2cos(ω2t+α2)+C3cos(ω3t+α3) where C1,C2,C3, α1, α2 and α3 are constants. By the use of different methods of constants variation in the previous solution to the linear basic system with equations (5.4) corresponding to the obtained system of governing equations (5.3) with additional restrictions to the system parameters, it is possible to obtain some partial approximation of the solution. Also, it is possible to study different dynamical properties of the system as well as some phenomena in the rheolinear, thermo-visco-elastic, nonlinear, hereditaryor creepproperties of themulti-pendulumsystemdefined in our paper. In the beginning, we consider solutions of the system with one (Fig.3a), and two-pendulum system (Fig.3b) with only one simpler set of the parallelly coupled standard light elements. Dynamics of multi-pendulum systems... 495 Fig. 3. Systemwith one pendulum coupled (a) and with two pendula intercoupled (b) by standard light creep element 6. Analytical solution to governing equations of the pendulum system for special cases 6.1. One pendulum oscillator In the case of one pendulum (Fig.3a), coupled of a fixed point by a set of parallel standard light elements: nonlinear spring, hereditary, creep and thermo-modified by temperature T(t), a differential equation of motion is given in the following form ϕ̈1+(ω̃ 2 0 +ω 2 0 +ω 2 00+ω 2 0h)ϕ1+ω 2 0T [1+γT̃(t)]ϕ1+2δϕ̇1 = (6.1) =ω20Th0T̃(t)+ ω̃ 2 0 (ϕ31 3! ) −εω20χ̃ϕ 3 1−ω 2 0αD α t [ϕ1(t)]+ω 2 0h t∫ 0 R(t− τ)ϕ1(t) dτ A solution to a separate case of the previous governing equation, can be ob- tained through the following simpler tasks: —Hereditary one pendulum oscillator (Fig.3a) governed by ϕ̈1+(ω̃ 2 0 +ω 2 0h)ϕ11 =ω 2 0h t∫ 0 R(t− τ)ϕ1(t) dτ (6.2) where R(t− τ)= ch− c̃h nc exp [ − 1 n (t− τ) ] which is analogous to the problem solved approximately, and the solution is presented in Goroško and Hedrih (2001) and also in Goroško and Hedrih (2007a,b). 496 K.R. Hedrih —Creep one pendulum oscillator (Fig.3a) governed by an ordinary fractional- order differential equation in the form ϕ̈1+(ω̃ 2 0 +ω 2 00)ϕ1 =−ω 2 0αD α t [ϕ1(t)] (6.3) In the case when α ∈ (0,1), we solve previous ordinary fractional-order differential equation (6.3) through Laplace’s transformations. After transfor- ming previous ordinary fractional-order differential equation (6.3) with the fractional-order derivative and having in mind that we introduced the nota- tions L{ϕ1(t)} for Laplace’s transformations as well as L {dαϕ1(t) dtα } = pαL{ϕ1(t)}− dα−1ϕ1(t) dtα−1 ∣∣∣ t=0 = pαL{ϕ1(t)} (6.4) and also having inmind that we accepted the hypothesis that the initial con- ditions of the fractional order derivatives of the system are given through the use of: dα−1ϕ1(t)/dt α−1|t=0 =0 and that L {d2ϕ1(t) dt2 } = p2L{ϕ1(t)}− [pϕ01+ ϕ̇01] (6.5) where ϕ01 and ϕ̇01 are the initial conditions of the system we can write the following solution to the equation with unknown Laplace’s transform L{ϕ1(t)}= pϕ01+ ϕ̇01 p2+ω20αp α+ ω̃20 +ω 2 00 (6.6) To obtain the inverse to the Laplace transform, we can use the result by Gorenflo and Mainardi (2000) as well as by Hedrih (2006a). For that reason and for the case when ω̃20+ω 2 00 6=0, we rewrite the previous expression in the following form L{ϕ1(t)}=(pϕ01+ ϕ̇01) 1 p2 [ 1+ ω20α p2 ( pα+ ω̃20 +ω 2 00 ω20α )] −1 = (6.7) = ( ϕ01+ ϕ̇01 p )1 p [ 1+ ω20α p2 ( pα+ ω̃20 +ω 2 00 ω20α )] −1 Then Laplace transform solution (5.7) can be expanded into series by the following way L{ϕ1(t)}= ( ϕ01+ ϕ̇01 p )1 p ∞∑ k=0 (−1)kω2k0α p2k ( pα+ ω̃20 +ω 2 00 ω20α )k (6.8) Dynamics of multi-pendulum systems... 497 or L{ϕ1(t)}= ( ϕ01+ ϕ̇01 p )1 p ∞∑ k=0 (−1)kω2k0α p2k k∑ j=0 ( k j ) pαjω 2(j−k) 0α (ω̃20 +ω 2 00) j (6.9) In (6.8), it is assumed that the expansion leads to convergent series. The inverse Laplace transform of the previous transform of solution (6.9) in term- by-term steps is based on theknown theorem, andyields the following solution to differential equation (6.3) of the time function in the following form of time series ϕ1(t)= L −1 L{ϕ1(t)}= =ϕ01 ∞∑ k=0 (−1)kω2k0αt 2k k∑ j=0 ( k j ) ω 2j 0αt −αj (ω̃20 +ω 2 00) jΓ(2k+1−αj) + (6.10) +ϕ̇01 ∞∑ k=0 (−1)kω2k0αt 2k+1 k∑ j=0 ( k j ) ω −2j 0α t −αj (ω̃20 +ω 2 00) jΓ(2k+2−αj) or ϕ1(t)= L −1 L{ϕ1(t)}= ∞∑ k=0 (−1)kω2k0αt 2k · (6.11) · k∑ j=0 ( k j ) ω 2j 0αt −αj (ω̃20 +ω 2 00) j [ ϕ01 Γ(2k+1−αj) + ϕ̇01t Γ(2k+2−αj) ] For two special cases of the solution (for α = 0 and α = 1), we have classical conservative or nonconservative pendulum oscillators. By using expression (6.10) obtained for the time solution ϕ1(t) with corre- sponding particular solutions, we can conclude that the solution contains two particular solutions in the following forms T1(t)= ∞∑ k=0 (−1)kω2k0αt 2k k∑ j=0 ( k j ) ω 2j 0αt −αj (ω̃20 +ω 2 00) jΓ(2k+1−αj) (6.12) T̃1(t)= ∞∑ k=0 (−1)kω2k0αt 2k+1 k∑ j=0 ( k j ) ω −2j 0α t −αj (ω̃20 +ω 2 00) jΓ(2k+2−αj) which are two vibration ”creeping” modes, T1(t,α) and T̃1(t,α), of the fractional-order dynamical properties of the one-pendulum system. By using 498 K.R. Hedrih these particular solutions, we made a numerical experiment for characteristic cases. The ratio of the pendulum creep system kinetic parameters, coefficient of creeping standard light element with the constitutive relation expressed by the fractional-order derivative and the results are presented inFigures 4 and5. It is visible that some types of modes are present as in the longitudinal vibra- tions of the rod with changeable cross-sections and built by a creep material with the stress strain constitutive relation expressed by the fractional-order derivative (see Hedrih, 2004c, 2005c). Fig. 4. Time function surfaces of T1(t,α), (6.12)1, for different kinetic and creep parameters of the-one pendulum system; (a) for ω2 0α /(ω̃2 0 +ω2 00 )=1, (b) for ω2 0α /(ω̃2 0 +ω2 00 )= 1/16, (c) for ω2 0α /(ω̃2 0 +ω2 00 )= 1/9. (d) for ω2 0α /(ω̃2 0 +ω2 00 )=9 The time functions T1(t,α) and T̃1(t,α) are surfaces found for differentpa- rameters ofkinetic andstandard light creepelements in the space (T(t,α), t,α) for the interval 0¬α¬ 1. In Fig.4, numerical simulations and graphical presentations of the parti- cular solution mode T1(t,α), (6.12)1, of fractional-differential equation (6.3) for different kinetic parameters are presented. In Fig.5, the particular solution T̃1(t,α), (6.12)2, of fractional-differential equation (6.12)2 of the system for different kinetic parameters in the interval 0¬α¬ 1 are given. Dynamics of multi-pendulum systems... 499 Fig. 5. Time function surfaces of T̃1(t,α), (6.12)2, for different kinetic and creep parameters of the-one pendulum system; (a) for ω2 0α /(ω̃2 0 +ω2 00 )=1, (b) for ω2 0α /(ω̃2 0 +ω2 00 )= 1/16, (c) for ω2 0α /(ω̃2 0 +ω2 00 )= 1/9. (d) for ω2 0α /(ω̃2 0 +ω2 00 )=9 6.2. Creep double pendulum oscillator Creep double pendulum oscillator (Fig.3b), is governed by ordinary fractional-order differential equations in the form ϕ̈1+ ω̃ 2 0ϕ1−ω 2 00(ϕ2−ϕ1)=ω 2 0α{D α t [ϕ2(t)]−D α t [ϕ1(t)]} (6.13) ϕ̈2+ ω̃ 2 0ϕ2+ω 2 00(ϕ2−ϕ1)=−ω 2 0α{D α t [ϕ2(t)]−D α t [ϕ1(t)]} An analogy with the result presented for the chain system in Hedrih (2006a) is useful to obtain the solution. After applying Laplace’s transfor- mations to previous equations (6.13) and having in mind that we introduced the notations L{ϕk(t)}, k=1,2 as well as that L {dαϕk(t) dtα } = pαL{ϕk(t)}− dα−1ϕk(t) dtα−1 ∣∣∣ t=0 = pαL{ϕk(t)} k=1,2 (6.14) 500 K.R. Hedrih and also having in mind that we accepted the hypothesis that the ini- tial conditions of the fractional-order derivatives of the system are given as (dα−1ϕk(t)/dt α−1)|t=0 =0 as well that L {d2ϕk(t) dt2 } = p2L{ϕk(t)}− [pϕ0k+ ϕ̇0k] k=1,2 (6.15) where ϕ0k and ϕ̇0k, k=1,2 are the initial conditions of the double-pendulum system, we can write the following system equations with unknown Laplace’s transforms (1+v)L{ϕ1(t)}−L{ϕ2(t)}= pϕ01+ ϕ̇01 ω20αp α+ω200 =h1 (6.16) −L{ϕ1(t)}+(2+v)L{ϕ2(t)}−L{ϕ3(t)}= pϕ02+ ϕ̇02 ω20αp α+ω200 =h2 we introduce v= p2+ ω̃20 ω20αp α+ω200 (6.17) The determinant of the previous system is ∆= (p2+ ω̃20 +2ω 2 0αp α+2ω200)(p 2+ ω̃20) (ω20αp α+ω200) 2 (6.18) Solutions to system equations (6.16) with respect to L{ϕk(t)}, k = 1,2, i.e. Laplace transforms of fractional-order differential equations (6.13) are in the following forms L{ϕ1(t)}= (pϕ01+ ϕ̇01)(p 2+ω20αp α+ ω̃20 +ω 2 00)+(pϕ02+ ϕ̇02)(ω 2 0αp α+ω200) (p2+ ω̃20 +2ω 2 0αp α+2ω200)(p 2+ ω̃20) (6.19) L{ϕ2(t)}= (pϕ02+ ϕ̇02)(p 2+ω20αp α+ ω̃20 +ω 2 00)+(pϕ01+ ϕ̇01)(ω 2 0αp α+ω200) (p2+ ω̃20 +2ω 2 0αp α+2ω200)(p 2+ ω̃20) For special cases of the double-pendulum system initial conditions, when at the initial moment the second pendulum is in the equilibrium position, the solutions are L{ϕ1(t)}= (pϕ01+ ϕ̇01)(p 2+ω20αp α+ ω̃20 +ω 2 00) (p2+ ω̃20 +2ω 2 0αp α+2ω200)(p 2+ ω̃20) (6.20) L{ϕ2(t)}= (pϕ01+ ϕ̇01)(ω 2 0αp α+ω200) (p2+ ω̃20 +2ω 2 0αp α+2ω200)(p 2+ ω̃20) Dynamics of multi-pendulum systems... 501 Taking into account that the sum and difference between the solutions to (6.16), the Laplace transforms of fractional-order differential equations (6.13) are L{ξ1(t)}= L{ϕ1(t)+ϕ2(t)}= pϕ01+ ϕ̇01 p2+ ω̃20 + pϕ02+ ϕ̇02 p2+ ω̃20 (6.21) L{ξ2(t)}= L{ϕ1(t)−ϕ2(t)}= (pϕ01+ ϕ̇01)− (pϕ02+ ϕ̇02) p2+ ω̃20 +2ω 2 0αp α+2ω200 The inverse Laplace transform to L{ξ1(t)} of the sum ϕ1(t)+ϕ2(t) of so- lution (6.21)1 yields the following sumof solutions to the systemof differential equations (6.13) ξ1(t)=ϕ1(t)+ϕ2(t)= L −1 L{ϕ1(t)+ϕ2(t)}= (6.22) = (ϕ̇01 ω̃0 + ϕ̇02 ω̃0 ) sin(ω̃0t)+(ϕ01+ϕ02)cos(ω̃0t) The inverse Laplace transform to L{ξ2(t)} of the difference ϕ1(t)−ϕ2(t) of solution (6.21)2 yields the following difference of solutions differential equ- ations (6.13) ξ2 =ϕ1(t)−ϕ2(t)= L −1 L{ϕ1(t)−ϕ2(t)}= L −1 {p(ϕ01−ϕ02)+(ϕ̇01− ϕ̇02) p2+2ω20αp α+ ω̃20 +2ω 2 00 } (6.23) We can see that the obtained Laplace transform L{ϕ1(t)−ϕ2(t)} is the same as the Laplace transform of the solution for the case of one pendulum withonecreep standard light element expressedby (6.6), andthen it is possible to use previous expression (6.7) and expansions of series (6.8) and (6.9) aswell as (6.10) for obtaining the corresponding solution to the necessary modes for the double-pendulum system by replacing the following parameters ϕ01 →ϕ01−ϕ02 ϕ̇01 → ϕ̇01− ϕ̇02 ω20α → 2ω 2 0α ω̃ 2 0 +ω 2 00 → ω̃ 2 0 +2ω 2 00 In the case when ω̃20+2ω 2 00 6=0, the Laplace transform can be expanded into following series L{ξ2}= L{ϕ1(t)−ϕ2(t)}= (6.24) = ( ϕ01−ϕ01+ ϕ̇01− ϕ̇02 p )1 p ∞∑ k=0 (−1)k2kω2k0α p2k k∑ j=0 ( k j ) pαj2(j−k)ω 2(j−k) 0α (ω̃20 +2ω 2 00) j 502 K.R. Hedrih In (6.24) it is assumed that the expansion leads to convergent series. The inverse Laplace transform of the previous Laplace transform of solution (6.24) in term-by-term steps is based on the known theorem, and yields the following solution to differential equations (6.13) of the time function in the following form of time series ξ2(t)= L −1 L{ϕ1(t)−ϕ2(t)}= (6.25) = (ϕ01−ϕ02) ∞∑ k=0 (−1)k2kω2k0αt 2k k∑ j=0 ( k j ) 2jω 2j 0αt −αj (ω̃20 +2ω 2 00) jΓ(2k+1−αj) + +(ϕ̇01− ϕ̇02) ∞∑ k=0 (−1)k2kω2k0αt 2k+1 k∑ j=0 ( k j ) 2−jω −2j 0α t −αj (ω̃20 +2ω 2 00) jΓ(2k+2−αj) For two special cases of the solution (for α = 0 and α = 1), we have classical conservative ornonconservative (Hedrih,2006d)pendulumoscillators, respectively. Byusingexpressionobtained for the timesolution ξ2(t)withcorresponding particular solutions, we can conclude that the solution contains two particular solutions in the following forms T1(t)= ∞∑ k=0 (−1)k2kω2k0αt 2k k∑ j=0 ( k j ) 2jω 2j 0αt −αj (ω̃20 +2ω 2 00) jΓ(2k+1−αj) (6.26) T̃1(t)= ∞∑ k=0 (−1)k2kω2k0αt 2k+1 k∑ j=0 ( k j ) 2−jω −2j 0α t −αj (ω̃20 +2ω 2 00) jΓ(2k+2−αj) which are two vibration ”creeping” modes, T1(t,α) and T̃1(t,α), of the frac- tional one-pendulum system oscillations. By using these particular solutions, wemade a numerical experiment for characteristic cases. The ratio of the pen- dulum creep system kinetic parameters, coefficient of creeping standard light element with the constitutive relation expressed by the fractional order deri- vative and the results are presented in Figures 4 and 5. It is visible that some types of modes is present as in the longitudinal vibrations of the rod with changeable cross-sections and built by a creep material with the stress strain constitutive relation expressed by the fractional order derivative presented in Hedrih and Filipovski (2002). Time functions T1(t,α) and T̃1(t,α) are surfaces found for different kinetic and standard light creep element parameters in the space (T(t,α), t,α) for the interval 0¬α¬ 1. Dynamics of multi-pendulum systems... 503 InFig.4, numerical simulations and graphical presentations of the particu- lar solutionmode T1(t,α) of fractional-differential equation (6.13) for different system kinetic parameters are presented. In Fig.5, the time particular solu- tionmode, T̃1(t,α), of fractional-differential equation (6.13) of the system for different system kinetic parameters in the interval 0¬α¬ 1 are shown. Solution to the normal modes of the system of fractional-order differential equations (6.13) are in the form ξ1(t)=ϕ1(t)+ϕ2(t)= L −1 L{ϕ1(t)+ϕ2(t)}= = (ϕ̇01 ω̃0 + ϕ̇02 ω̃0 ) sin(ω̃0t)+(ϕ01+ϕ02)cos(ω̃0t) (6.27) ξ2(t)=ϕ1(t)−ϕ2(t)= = ∞∑ k=0 (−1)k2kω2k0αt 2k k∑ j=0 ( k j ) 2jω 2j 0αt −αj (ω̃20 +2ω 2 00) j [ ϕ01−ϕ02 Γ(2k+1−αj) + (ϕ̇01− ϕ̇02)t Γ(2k+2−αj) ] Then the solutions to equations (6.13) are ϕ1(t)= 1 2 [(ϕ̇01 ω̃0 + ϕ̇02 ω̃0 ) sin(ω̃0t)+(ϕ01+ϕ02)cos(ω̃0t) ] + + 1 2 ∞∑ k=0 (−1)k2kω2k0αt 2k k∑ j=0 ( k j ) 2jω 2j 0αt −αj (ω̃20 +2ω 2 00) j [ ϕ01−ϕ02 Γ(2k+1−αj) + (ϕ̇01− ϕ̇02)t Γ(2k+2−αj) ] (6.28) ϕ2(t)= 1 2 [(ϕ̇01 ω̃0 + ϕ̇02 ω̃0 ) sin(ω̃0t)+(ϕ01+ϕ02)cos(ω̃0t) ] + − 1 2 ∞∑ k=0 (−1)k2kω2k0αt 2k k∑ j=0 ( k j ) 2jω 2j 0αt −αj (ω̃20 +2ω 2 00) j [ ϕ01−ϕ02 Γ(2k+1−αj) + (ϕ̇01− ϕ̇02)t Γ(2k+2−αj) ] 6.3. Thermo-visco-elastic double pendulum oscillator excited by random temperature By using the results by Hedrih (2007e) presented at DSTA 2007, we can add some discussion with respect to the mode of the thermo-visco-elastic double-pendulum system with comparison to the one-pendulum system both excited by random temperature applied to the thermo-elastic standard light element. 504 K.R. Hedrih For the double-pendulum system, the governing rheolinear equations expressed in terms of normal coordinates of the corresponding linear system can be obtained in the following form ξ̈1+ ω̃ 2 0ξ1 =0 (6.29) ξ̈2+ ω̃ 2 0ξ2+2ω 2 0ξ2+2ω 2 0T [1+γT̃(t)]ξ2+4δξ̇2 =−2ω 2 0Th0T̃(t) with two eigen-frequencies ω21,2 = ω̃ 2 0 +ω 2 0 +ω 2 0T ∓ (ω 2 0 +ω 2 0T) of the linear system. Fig. 6. (a) Systemwith one free pendulum. (b) Thermo-rheological system; partial oscillator 1 For the double pendulum system, the first equation of rheo-nonlinear sys- tem (6.29) in the linearised form represents a pure partial harmonic oscillator presented in Fig.6a and 6b, with the eigen frequency ω21 = ω̃ 2 0 = g/ℓ of free one-mode vibrations. This linearised case is when both pendula oscillate with the same frequency, ω̃20 = g/ℓ, as the decoupled pendula (single ma- thematical pendula). Then, the standard light thermo-visco-elastic element thermo-modified by temperature T(t) does not influence this normal coordi- nate composed by sum ξ1 = ϕ1 +ϕ2. Along this normal (main) coordinate, the oscillation in the linearised approximation is free, without temperature influence. This is right for all cases of the multi-pendulum systems presented in Fig.6b. For the double-pendulum system, the second equation of rheo-nonlinear system (4.4), on the normal coordinate ξ2 = ϕ1 −ϕ2 in the linearised form is the Mathieu-Hill equation, and represents mathematical description of the thermo-rheological oscillator presented inFig.7a or 7c,with parallelly coupled two light standard thermo-visco-elastic elements thermo-modified by the same temperature T(t) and one linear elastic spring with rigidity c0 = mg/ℓ. For the coordinate ξ2 =ϕ1−ϕ2, we can separate twomain cases. For both cases, we take into consideration the asymptotic approximation of the amplitude and Dynamics of multi-pendulum systems... 505 Fig. 7. (a) Systemwith one pendulum coupled by the standard light visco-thermoelastic element. (b) Systemwith two pendula intercoupled by the standard light thermorheological element. (c) Thermo-rheological system; partial oscillator 2 phase of the dynamic process on this coordinate ξ2 = ϕ1 −ϕ2 close around, firstly, main resonance when Ω ≈ ω2 = √ ω̃20 +2(ω 2 0 +ω 2 0T) and, secondly, around the parametric resonance when Ω ≈ 0.5ω2 = 0.5 √ ω̃20 +2(ω 2 0 +ω 2 0T). Then, we can conclude that along this coordinate under the corresponding kinetic parameters, there can appear, firstly, regimes closest to themain reso- nant state as well as onemain resonant state, and secondly, regimes closest to the parametric resonant state as well as one resonant state under the thermo- visco-elastic temperature single frequency excitation. This second mode has the same character as vibration of the one-pendulum system presented in Fig.7a. For details, see Hedrih (2007e). 7. Concluding remarks We can conclude that between multi-pendulum systems and chain dynami- cal systems there exists a mathematical analogy in descriptions as well as in vibration phenomena depending on the character of standard light coupling elements between pendula or, analogously, between material particles in the chain. Also, there is a mathematical analogy between corresponding modes in a multi-beam system or multi-plate systems with the corresponding cha- 506 K.R. Hedrih racter of a light distributed coupling layer between the beams or plates in multi-deformable body systems. Themathematical description leads to the same ordinary differential equ- ations, or ordinary integro-differential or fractional-order differential equations governing both analogous types of problems. For a homogeneous sandwich multi-plate, or a multi-beam system, it is possible to identify some analogies between with mechanical multi-material particle chains andmulti-pendulumsystemswith interconnections by standard light elements of different properties. Acknowledgment Some parts of this researchwere supported by theMinistry of Sciences andEnvi- ronmental Protection of Republic of Serbia through Mathematical Institute SANU, Belgrade, Grant ON144002 ”Theoretical and Applied Mechanics of Rigid and Solid Body. Mechanics of Materials” and Faculty of Mechanical Engineering University of Nǐs. References 1. Enelund M., 1996,Fractional Calculus and Linear Viscoelasticity in Structu- ral Dynamics, Division of Solid Mechanics, Chalmers Tekniska Hogskola, Go- teborg, Sweden, 1-27 2. Gorenflo R., Mainardi F., 2000, Fractional calculus, integral and differen- tial equations of fractional order, CISM Lecture Notes, Udine, Italy, Preprint 54 pages, 223-276 3. Goroshko O.A., Puchko N.P., 1997, Lagrangian equations for the multi- bodies hereditary systems, Facta Universitatis, Series Mechanics, Automatic Control and Robotics, 2, 7, 209-222 4. Goroško O.A., Hedrih (Stevanović) K., 2001,Analitička dinamika (me- hanika) diskretnih naslednih sistema (Analytical Dynamics (Mechanics) of Di- screte Hereditary Systems), University of Nǐs, Monograph, p.426, YU ISBN 86-7181-054-2 5. Goroško O.A., Hedrih (Stevanović) K., 2007a, Construction of the La- grange’smechanics of the hereditary systems,APMSaint Petersburg 2007 pp., Minisymposium Oppening Lecture, The International Summer School APM – Advanced Problem in Mechanics, Saint Petersburg, 133-156 6. Goroško O.A., Hedrih (Stevanović) K., 2007b, Construction of the La- grange’s mechanics of the hereditary systems, Facta Universitatis Series Me- chanics, Automatic Control and Robotics, 6, 1, 1-23 Dynamics of multi-pendulum systems... 507 7. Hedrih (Stevanović)K., 1999,Thermorheological hereditary pendulum, In: Thermal Stresses 99, Edited by J.J. Skrzypek and R.B. Hetnarski, Cracow, 199-202 8. Hedrih (Stevanović) K., 2001, Differential equations of two mass partic- les, constrainedwith a piezo-thermo-rheological hereditary element, dynamics, Proceedings of full papers, 5th International Conference on Applied Electroma- gnetics, PES 2001, Edited by D. Veličković, 77-80 9. Hedrih (Stevanović) K., 2002a, The dissipation function of a nonconserva- tive system of mass particles,Tensor, N.S., 63, 2, 176-186, ISSN 0040-3504 10. Hedrih (Stevanović)K., 2002b,Transversal creepvibrations of a beamwith fractionalderivativeconstitutive relationorder. I –Partial fractional-differential equation. II – Stochastic stability of the beam dynamic shape, under axial bounded noise excitation, Proceedings of Forth International Conference on Nonlinear Mechanics (ICNM-IV), Edited byWei Zang Chien et al., Shanghai, P.R. China, 584-595 11. Hedrih (Stevanović) K., 2003a, Discrete continuum’s models and thermo- rheological elements – basic idea and tensor equations, homogeneous linear cha- in and plane/spacematerial nets,Proceedings of full papers, 6 th International Conference on Applied Electromagnetics, PES 2003, Edited by D. Veličković, 127-130, 131-134 12. Hedrih (Stevanović) K., 2003b, The longitudinal and transversal creep vi- brations of a fractional order constitutive relation beams, Scientific Bulletin of the Politehnica, University of Timisoara, Transaction on Mechanics, 48, 62, 5-12, 13-22, ISSN 1224-6077, http://www.utt.ro/english/pbseng.shtml 13. Hedrih (Stevanović) K., 2004a, Creep vibrations of a fractional derivati- ve order constititive relation deformable bodies, Proceedings Eighth American Congress of AppliedMechanics PACAMVIII, LaHabana,Cuba,SeriesApplied Mechanics in Americas, 10, 548-551, ISBN 959-7056-20-8 14. Hedrih (Stevanović) K., 2004b, Discrete continuum method, Computatio- nal Mechanics, WCCM VI in conjunction with APCOM’04, Beijing, China, Tsinghua University Press & Springer-Verlag, 1-11, CD 15. Hedrih (Stevanović) K., 2004c, Partial fractional order differential equ- ations of transversalvibrationsof creep connecteddoubleplates systems,Work- shop Preprints/Proceedings No. 2004-1 IFAC FDA 04, ENSEIRB, Bordeaux France, 299-304 16. Hedrih (Stevanović) K., 2005a, Eigen amplitude vectors and functions extended orthogonality of small oscillations mixed systems of the coupled di- screte and continuous subsystems, Facta Universitatis, Series Mechanics, Au- tomatic Control and Robotics, 4, 17, 225-243, YU ISSN 0534-2009 17. Hedrih (Stevanović) K., 2005b, Integrity of dynamical systems, Journal Nonlinear Analysis, 63, 854-871 508 K.R. Hedrih 18. Hedrih (Stevanović) K., 2005c, Partial Fractional order differential equ- ations of transversal vibrations of creep-connected double plate systems, In: Fractional Differentiation and its Applications, Edited by A. Le Mahaute, J.A.T. Machado, J.C. Trigeassou and J. Sabatier, U-Book, 289-302 19. Hedrih (Stevanović) K., 2006a, Modes of the homogeneous cha- in dynamics, Signal Processing, 86, 2678-2702, ISSN: 0165-1684, www.sciencedirect.com/science/journal/01651684 20. Hedrih (Stevanović) K., 2006b, The frequency equation theorem of small oscillations of a hybrid system containing coupled discrete and continuous sub- systems, Facta Universitatis Series: Mechanics, Automatic Control and Robo- tics, 5, 1, 25-41, http://facta.junis.ni.ac.yu/facta/ 21. Hedrih (Stevanović) K., 2006c, The transversal creeping vibrations of a fractional derivative order constitutive relation of nonhomogeneous beam,Ma- thematical Problems in Engineering, Special issue: Nonlinear Dynamics and their Applications in Engineering Sciences, Edit.: J.M. Barhesar, 2006, 5, 61- 78, www.hindawi.com 22. Hedrih (Stevanović) K., 2006d, Transversal forced vibrations of an axial- ly moving sandwich belt system, Archive of Applied Mechanics, Springer, http://springerlink.com/content/?k=Hedrih 23. Hedrih (Stevanović) K., 2006e, Transversal vibrations of double-plate sys- tems,Acta Mechanica Sinica, 22, 487-501 24. Hedrih (Stevanović) K., 2006f, Transversal vibrations of the axial- ly moving sandwich belts, Archive of Applied Mechanics, 77, 7, 523-539, http://springerlink.com/content/?k=Hedrih 25. Hedrih (Stevanović) K., 2007a, Double plate system with discontinuity in the elastic bonding layer,Acta Mechanica Sinica, 23, 2, 221-229 26. Hedrih (Stevanović)K., 2007b,Energyanalysis in thenonlinear hybrid sys- tem containing linear and nonlinear subsystem coupled by hereditary element, Nonlinear Dynamics, 51, 1, 127-140 27. Hedrih (Stevanović) K., 2007c, Energy analysis of the double plate system, Acta Mechanica Sinica, DOI 10.1007/s10409-007-0124-z 28. Hedrih (Stevanović) K., 2007d, Hybrid systems and hybrid dynamics: the- ory and applications, Invited Plernary Lecture, 8th HSTAM International Con- gress on Mechanics, Patras, Greece, Edited by N. Bazwos, D.L. Karabalis, D. Polyzos, D.E. Beskos and J.T. Katsikadelis, I, 77-86 29. Hedrih (Stevanović)K., 2007e, Stochastic dynamics of hybrid systemswith thermo-rheological hereditary elements,Proceedings of IX International confe- rence on Dynamical Systems – Theory and Application, Edited byK.Awrejce- wicz, P. Olejnik andMrozowski, Łódź, Poland, University of Łódź, 193-202 Dynamics of multi-pendulum systems... 509 30. Hedrih (Stevanović) K., Filipovski A., 2002, Longitudinal vibration of a fractional derivative order rheological rod with variable cross section, Fac- ta Universitatis, Series Mechanics, Automatic Control and Robotics, 3, 12, 327-350, YU ISSN 0534-2009, http://facta.junis.ni.ac.yu/facta/macar/macar2002/macar2002-02.html http://facta.junis.ni.ac.yu/facta/macar/macar200501/macar200501-04.html 31. RaškovićD., 1965,Teorija oscilacija (TheoryofOscillations),NaučnaKnjiga, 503p. [in Serbian] 32. Rzhanitsin A.R., 1949, Some Questions of the Mechanics of Deforming in Time Systems, Moscow, GTTI, 248p. [in Russian] 33. SavinG.N., RuschiskyYu.Ya., 1976,Elements of HereditaryMediaMecha- nics, Kyiv, Vyscha Shkola, 250p. [in Ukrainian] 34. TorvikP.J., BagleyR.L., 1984,On the appearance of the fractional deriva- tives in thebehaviorof realmaterials,Journal ofAppliedMechanics (Trasaction ASME), 51, 294-298 Dynamika układów wielowahadłowych z efektem pełzania opisanym elementami ułamkowego rzędu Streszczenie Wpracy zaprezentowano krótki przegląd rezultatów badań autora nad dynamiką układów hybrydowych i dyskretnych, złożonych z punktów materialnych sprzęgnię- tych standardowymi elementami odpowiadającymi za pełzaniewmateriale i opisywa- nychpochodnąułamkowegorzędu.Rozważonodrgania swobodneukładówwielowaha- dłowych z elementami o różnychwłaściwościach zdefiniowanych równaniempomiędzy stanem naprężenia a odkształcenia. Wyprowadzone równania różniczkowo-całkowe ułamkowego rzędu rozwiązano analitycznie. Przedstawiono szczegółowo przypadek układu z pojedynczymwahadłem i układu dwuwahadłowego zawierającego elementy pełzania opisane równaniemkonstytutywnym stanu naprężenia i odkształcenia o rzę- dzie ułamkowym.Na podstawie otrzymanych rozwiązań analitycznych zauważono, że drgania swobodne wykazują charakter okresowy i nieokresowy, przy czym te ostatnie mają dwa różne przebiegi (w tymprzypadku rozwiązanie podanowpostaci rozwinięć w szeregi potęgowe). Wyniki badań teoretycznych i numerycznych różnego rodzaju drgań przy zmiennych parametrach kinetycznych tych układów przedstawiono gra- ficznie. Manuscript received January 18, 2008; accepted for print January 31, 2008