Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 52, 3, pp. 687-697, Warsaw 2014 WAVELET APPROXIMATION OF ADOMIAN’S DECOMPOSITIONAPPLIED TO THE NONLINEAR PROBLEM OF A DOUBLE-BEAM RESPONSE SUBJECT TO A SERIES OF MOVING LOADS Piotr Koziol Koszalin University of Technology, Koszalin, Poland e-mail: piotr.koziol@wbiis.tu.koszalin.pl The dynamic response of a double-beam resting on a nonlinear viscoelastic foundation and subjected to a finite series of moving loads is analysed. The beams are connected by a viscoelastic layer and the load moving along the upper beam represents motion of a train on the rail track. Themathematical model is described by a coupled system of fourth order partial differential equations with homogeneous boundary conditions. The nonlinearity is included in the foundation stiffness of medium supporting a lower beam. The coiflet based approximation combined with Adomian’s decomposition is adopted for the displacements derivation. The developed approach allows one to overcome difficulties related to direct calculation of Fourier integrals as well as the small parameter method. The conditions for correctness of the approximate solution are defined. The influence of some factors on the system sensitivity is discussed, with special focus on the distance between the separated loads. Numerical examples are presented for a certain system of physical parameters. Keywords: infinite double-beam, nonlinear problem, Adomian’s decomposition, coiflet approximation, moving load 1. Introduction The subject of moving load problems, especially the analysis of beam vibrations arising from dynamic excitations, is a very important direction of studies associated with modern railway engineering. New conditions, never met in the past, such as high speed rails or intensive freight transport, lead tonecessity of bettermodelling andprediction of possible scenario for operational transport systems (Bogacz and Frischmuth, 2009; Bogacz and Krzyzynski, 1991; Krylov, 2001; Thompson, 2009). In order to achieve this aim, new computational approaches are needed. The new techniques should allow effective parametrical analysis and show phenomena that can appear in certain situations. One can find a number of published results showing closed form solutions for special cases of linearmodels representing the rail track on a rigid foundation. The presented analyses usually involve classical approaches such as Fourier and Laplace transforms or the Green function technique (Fryba, 1999; Hussein and Hunt, 2006; Oniszczuk, 2003). It was shown that the nonlinearmodel reflects behaviour of the beam deflection better than the linearmodelwhen comparedwithmeasurements (Dahlberg, 2002). In the case of a nonlinear foundation, classical procedures such as the finite elementmethod or the perturbation approach are often employed, and they usually give results that are insufficiently exact, making unable an efficient parametrical study (Kargarnovin et al., 2005; Kuo and Lee, 1994; Sapountzakis and Kampitsis, 2011; Wu and Thompson, 2004). The coiflet based approximation combined with Adomian’s decomposition (Adomian, 1989; Koziol, 2010; Wang et al., 2003), adapted in this paper to thedoublebeamproblemsolution, appearedmoreefficient and leading to reliable results enabling parametrical analysis of nonlinear systems of the beam-foundation type (Hryniewicz and Kozioł, 2013; Koziol and Hryniewicz, 2012). The method has been validated in a number 688 P. Koziol of studies, e.g. by FEM and wavelet FEM in the case of the Euler-Bernoulli beam resting on a linear viscoelastic foundation (Musuva et al., 2014) and by the Lindestedt–Poincare regular perturbationmethod in the case of the Euler-Bernoulli beam resting on a nonlinear foundation (Kargarnovin, 2005; Koziol, 2013a) both subjected to moving loads. This paper introduces nonlinear assumptions to the linear model of a double-beam (Onisz- czuk, 2003; Janget al., 2008) assumingthat the rail andslabsaremodelledby theEuler-Bernoulli beam and the bottom layer has nonlinear characteristics. Preliminary results for the proposed model solved by coiflet estimation were presented at the Railways 2012 conference (Hryniewicz and Koziol, 2012), and this paper is an essential extension of the past work. It is assumed that the load is represented by a finite series of harmonically varying loads distributed on separated intervals, this being a reliable representation of train movement. The study carried out shows that the developedmethod is efficient enough for the parametrical analysis. The undertaken ef- fort is focused on the investigation of the influence of the distance between certain loadsmoving along the upper beam on the nonlinear response of the structure, and also on difficulties in the method application arising from nonlinear assumptions. It is shown that the conditions needed for a solution exact enough in the linear case (Koziol, 2010; Wang et al., 2003) are insufficient for the nonlinear system. 2. Double-beam on a nonlinear foundation The discussed problem of the dynamic response of the double-beam system supported by a viscoelastic nonlinear foundation subject to a series of loads moving along the upper beam is very important in railway engineering. There is a need for structural dynamics prediction for constructions associated with train transportation, and these can be modelled analytically or numerically. The analysed model is related to the rail track consisting of rails, railpads, the floating slab and slab bearings. The theoretical model is composed of an upper beam to account for the rails (with the mass m1 and the bending stiffness EI1), a lower beam representing the slab (with the mass m2 and the bending stiffness EI2) and two viscoelastic layers; the first one linear (with the stiffness k1 and the viscous damping factor c1) between the rail and the slab, and the second one nonlinear and supporting the structure (with the linear stiffness k2, the nonlinear part of stiffness kN and the viscous damping factor c2). The two couplednonlinearpartial differential equations ofmotiondescribevertical vibrations of the double-beam system EI1 ∂4U ∂x4 +m1 ∂2U ∂t2 + c1 (∂U ∂t − ∂W ∂t ) +k1(U−W)−kNW 3 =P(x,t) EI2 ∂4W ∂x4 +m2 ∂2W ∂t2 + c2 ∂W ∂t +k2W − c1 (∂U ∂t − ∂W ∂t ) −k1(U−W)+kNW 3 =0 (2.1) where U and W are the transverse displacements of the upper and lower beam, respectively, at position x and time t. The load moving along the upper beam is modelled by a finite series of loads harmonically varying in time and distributed on separated intervals in the space domain P(x,t)= L−1∑ l=0 P0 2a cos2 (π[x−Vt− (2a+s)l] 2a ) H ( a2− [x−Vt− (2a+s)l]2 ) eiΩt (2.2) where H(·) is the Heaviside step function, 2a is the span of the moving load, Ω=2πfΩ is the frequency of the load, V is the velocity of the moving excitation and s is the distance between the separated loads. In order to obtain a steady-state response, amoving coordinate system can be introduced (x,t)→ (x=x−Vt,t) (2.3) Wavelet approximation of Adomian’s decomposition applied to the nonlinear problem... 689 The response in the following classical form is analysed in this paper U(x,t)=u(x)eiΩt W(x,t)=w(x)eiΩt (2.4) Applying representation (2.4) and new variables (2.3), one can rewrite equations (2.1) in form represented by differential operators Q4u+R1w= kNw 3e2iΩt+P(x) Q1u+R4w=−kNw 3e2iΩt (2.5) where Q1 =α1 d dx +α0 Q4 =α4 d4 dx4 +α2 d2 dx2 +α1 d dx +α0 R1 =β1 d dx +β0 R4 =β4 d4 dx4 +β2 d2 dx2 +β1 d dx +β0 (2.6) and P(x)= L−1∑ l=0 P0 2a cos2 (π[x− (2a+s)l] 2a ) H ( a2− [x− (2a+s)l]2 ) (2.7) with the new coefficients α4 =EI1 α2 =m1V 2 α1 =−V (c1+2im1Ω) α0 = k1−m1Ω 2+ic1Ω α1 = c1V α0 =−k2− ic1Ω β1 =−c1V β0 =−k1+iΩ(c2− c1) β4 =EI2 β2 =m2V 2 β1 =−V (c1+c2+2im2Ω) β0 = k1+k2−m2Ω 2+ic1Ω (2.8) The deflection, slope and curvaturemust tend to zero far from the excitation and, therefore, the boundary conditions can be assumed as follows lim x→±∞ u(x)= lim x→±∞ w(x) = 0 lim x→±∞ du dx = lim x→±∞ dw dx =0 lim x→±∞ d2u dx2 = lim x→±∞ d2w dx2 =0 (2.9) 3. Adomian’s procedure The method developed by Adomian assumes that the solution to any nonlinear problem de- scribed by differential equations can be represented by an infinite series of functions (Adomian, 1994; Koziol, 2010; Koziol andHryniewicz, 2012). The first of these functions is a solution to the linear problem associated with the problem considered, and the other terms describe nonlinear factors influencing the response. The functions to be evaluated are represented byAdomian po- lynomials.One canfindmanyprocedures for the evaluation of these polynomials in the literature (Adomian, 1989; Adomian 1994; Hosseini and Nasabzadeh, 2006; Pourdarvish, 2006; Wazwaz, 1999;Wazwaz andEl-Sayed, 2001). Some of themoffer a bit better convergence in the sense that a smaller number of polynomials is needed for obtaining results with accuracy sufficient for the solution exact enough (Wazwaz, 1999;Wazwaz andEl-Sayed, 2001). The performed simulations showthat this feature is kept also for theproblemconsidered, but computational difficulties grow due to complexity of the formulas. Therefore, a good balance between computational cost and effectiveness must be foundwhen choosing the numerical procedure. In this paper, the following form of the polynomials is taken (Das, 2009; Pourdarvish, 2006) Aj(x)= 1 j! ( dj dλj ∞∑ k=0 λkwk(x) )∣∣∣∣ λ=0 for j=0,1,2, . . . (3.1) 690 P. Koziol and for practical calculations, the first four polynomials are used A0 =w 3 0 A1 =3w 2 0w1 A2 =3(w0w 2 1 +w 2 0w2) A3 =w 3 1 +6w0w1w2+3w 2 0w3 A4 =3(w 2 1w2+w0w 2 2 +2w0w1w3+w 2 0w4) (3.2) For the effective solution of the problem considered (Eqs. (2.5)), one assumes that the beams vibrations are represented by the following series u(x)= ∞∑ j=0 uj(x) w(x) = ∞∑ j=0 wj(x) (3.3) and the nonlinear cubic term w3(x) is characterized by the Adomian polynomials w3(x)= ∞∑ j=0 Aj(x) (3.4) The problemof the series convergence (Eqs. (3.3)) was solved anddiscussed in the literature. In order to control the accuracy of the Adomian approximation, one can use the convergence condition defined by the following formula (Hosseini and Nasabzadeh, 2006) 0< ‖uj+1‖ ‖uj‖ < 1 0< ‖wj+1‖ ‖wj‖ < 1 (3.5) with the norm ‖uj‖=max x |Re[uj(x)]| and ‖wj‖=max x |Re[wj(x)]|. A reliable representation of the response can be described by the n-th order approximation Sn(u,x)= n∑ j=0 uj(x) Sw(w,x) = n∑ j=0 wj(x) (3.6) depending on the accuracy assumed. Combining equations (2.5), (3.3) and (3.4) leads to an effective algorithm for the evaluation of the terms uj(x) and wj(x) (j=1,2,3, . . .) Q4u0+R1w0 =P(x) Q1u0+R4w0 =0 Q4uj +R1wj = kNAj−1(x)e 2iΩt Q1uj +R4wj =−kNAj−1(x)e 2iΩt (3.7) Thepractical formofAdomianpolynomials (3.2) needed for consecutive evaluation of appro- ximate series terms (Eq. (3.3)) is computed in this paperbyusingawavelet based approximation allowing derivation of the Fourier transformwithout numerical approach and, in the same time, keeping the assumed accuracy. Essential details of the adoptedmethod are presented in the next section. 4. The Fourier transform and the coiflet approximation The above system of equations (3.7) can be solved by applying the Fourier transforms f̂(ω)= ∞∫ −∞ f(x)e−iωx dx f(x)= 1 2π ∞∫ −∞ f̂(ω)eiωx dω (4.1) Thus, one can obtain the systemof formulas that are needed for theAdomian series computation in the transform domain û0(ω)= R̂4P̂(ω) Q̂4R̂4− Q̂1R̂1 ûj(ω)= kNÂj−1(ω)(R̂4+ R̂1) Q̂4R̂4− Q̂1R̂1 e2iΩt ŵ0(ω)= −Q̂1P̂(ω) Q̂4R̂4− Q̂1R̂1 ŵj(ω)= −kNÂj−1(ω)(Q̂4+ Q̂1) Q̂4R̂4− Q̂1R̂1 e2iΩt (4.2) Wavelet approximation of Adomian’s decomposition applied to the nonlinear problem... 691 where j=1,2,3, . . . and P̂(ω)=P0 iπ2[1−exp(2iaω)] 4aω(π2−a2ω2) L−1∑ l=0 exp[−i(a+2la+ ls)ω] (4.3) Formulas (4.2) and (4.3)must be retransformed to obtain the dynamic response of the system in the physical domain. Because of complexity of the integrands appearing in the Fourier integrals, classical ornumericalmethodsof the inverseFourier transformcalculationbecome ineffective and might give wrong results. An alternative method of derivation based on the coiflet expansion of functions is adopted in this paper (Koziol, 2010;Wang et al., 2003). It uses thefilter family called coiflets, possessing the property of vanishing moments for both wavelet and scaling functions, which allows one to estimate multiresolution coefficients relatively easy (Monzon et al., 1999; Wang et al., 2003) Ψ(x)= KC∑ k=0 (−1)kpKC−kΦ(2x−k) Φ(x)= KC∑ k=0 pkΦ(2x−k) (4.4) ΨC and ΦC are thewavelet function and the scaling function, respectively and KC is the number of filter coefficients belonging to the applied family of coiflets pk. Using the properties of coiflets and relations between the wavelets, the Fourier analysis leads to the formulas allowing one to approximate analytically the Fourier transforms of every function belonging to the L2 space f̂(ω)= lim n→∞ f̂n(ω)= lim n→∞ 2−n−1 K̃∏ k=1 ( KC∑ k=0 pke ikω/2n+k ) kmax(n)∑ k=kmin(n) f((k+M)2−n)e−iωk2 −n f(x)= lim n→∞ fn(x)= lim n→∞ 1 2n+2π K̃∏ k=1 ( KC∑ k=0 pke −ikx/2n+k ) kmax(n)∑ k=kmin(n) f̂((k+M)2−n)eixk2 −n (4.5) where M = ∑KC k=0kpk and the range of summation kmin(n) = ωmin2 n −KC − 1, kmax(n) = =ωmax2 n−1 must be determined in such a way that at least the interval [ωmin,ωmax] covers the set of the variable ω having strong influence on the behaviour of the original function that guaranties proper application in the linear case. In the next section, it will be shown that this criterion, given in previous publications for linear problems (Koziol, 2010; Wang et al., 2003), becomes insufficient for the problem considered. In this case, mainly due to the nonlinear characteristics of the system and the applied excitation consisting of a number of loads, the analysis of the power spectrum does not give all information needed for the determination of coiflet approximation, as opposed to thementioned criterion applied tomuch simpler problems, i.e. a single load or linear models. The proper condition appropriate for the model analysed is the stabilisation of the solutionwith an increasing order of the approximation combinedwith the analysis of the power spectrumproposed before and a deep knowledge about possible behaviour of the structure. Combining Adomian’s decomposition with the coiflet approximation (Eqs. (4.5)) gives a semi-analytical procedure for the displacements evaluation for j=0,1,2, . . . uj(x)= 1 2n+1π K̃∏ k=1 KC∑ j=0 pj exp(ijx/2 n+k) kmax(n)∑ k=kmin(n) ûj (k+M 2n ) exp(ikx/2n) wj(x)= 1 2n+1π K̃∏ k=1 KC∑ j=0 pj exp(ijx/2 n+k) kmax(n)∑ k=kmin(n) ŵj (k+M 2n ) exp(ikx/2n) Âj(ω)= 1 2n K̃∏ k=1 KC∑ j=0 pj exp(−ijω/2 n+k) kmax(n)∑ k=kmin(n) Aj (k+M 2n ) exp(−ikω/2n) (4.6) 692 P. Koziol An arbitrarily chosen point in the space can be taken for numerical simulations. Further analysis is carried out for x=0 and the n-th order Adomian’s approximation (Eqs. (3.6)) u(t)≈Sn(u,t)= n∑ j=0 uj(−Vt) w(t)≈Sn(w,t) = n∑ j=0 wj(−Vt) (4.7) 5. Simulations and discussion The following parameters are considered for numerical examples (Kargarnovin et al., 2005; Abu- Hilal, 2006; Hryniewicz and Koziol, 2013): KC = 17, K̃ = 10, n = 5, P0 = 5 · 10 4N/m, EI1 =10 7Nm2,m1 =100kg/m, k1 =4 ·10 7N/m2, c1 =6.3 ·10 3Ns/m2,EI2 =1.43 ·10 9Nm2, m2 = 3.5 ·10 3kg/m, k2 = 5 ·10 7N/m2, c2 = 4.18 · 10 4Ns/m2, kN = 4 ·10 13N/m4, fΩ = 5Hz, V =20m/s. One shouldnote that theused systemofparameters is adjusted in order tohighlight themain features of the approach developed, andmore realistic assumptions could be considered for real structures behaviour investigation. Nevertheless, the implemented values secure the convergence of the approximate solutions and show that the adopted method indeed enables parametrical analysis of the considered nonlinear dynamic system. Answers to questions: how to describe the deflection of the beams after the load passed the observation point and how to characterize the behaviour of extreme values of the vibrations amplitude are of importance in railway engineering. For this study, the complexmodulus of the solution can be examined, called in the literature “themaximal response” (Kim andCho, 2006; Koziol andHryniewicz, 2012; Sun, 2002). Themaximal response describes changes of the system dynamic sensitivity in time. This kind of approach is helpful in investigation of the nonlinear system because for such cases, the real and imaginary part of the solution cannot be so easily interpreted as it is possible for linear modelling. Numerical simulations show that the 5-th order Adomian approximation can be taken for the analysis as a reliable representation of the solution for the considered model and system of parameters (Eq. (3.6)). The order of the coiflet estimation can be fixed as n=5, and above this index the solution stabilizes without showing significant changes when n increases. Figures 1 and 2 present how the beamsdeflection changes in time, depending on the distance between the separated loads being parts of the moving excitations (Eq. (2.2)). One can observe that the nonlinear solution differs from the linear onemuch stronger when the distance s beco- mes smaller. The nonlinear response accumulates in such cases, especially for the upper beam representing the rail. The response of the upperbeam is stronger than vibrations of the lower be- am in each case. Themore noticeable nonlinear effects can also be observed for the upper beam, whichmeans that the influence of the nonlinear factor included in the layer supporting the lower beampropagates through the system to the surfacewhere it causes the strongest variations.One can see that for the lower beam, the difference between the linear and the nonlinear response is hardly visible (Figs. 1d-f). It was shown previously that the increasing number of loads makes the nonlinear effects stronger in the systems similar to the model considered (Koziol, 2013b). There aremany factors with an important influence on the developed approximate solution. A good balance between the efficiency and cost effectiveness of the algorithm must be found in order to achieve a good procedure allowing parametrical analysis. Themost important features of this approachwere discussed in past publications, e.g. (Koziol, 2010; Koziol 2013b). However, in the case of a nonlinear problem, newquestions concerning thewavelet approximation appear. Figures 3 and 4 show that the criterion regarding the range of summation in the coiflet approximation formula (Eqs. (4.6)) formulated previously is insufficient for the double-beam system. It was shown that for the linear modelling it is enough to determine the range of Wavelet approximation of Adomian’s decomposition applied to the nonlinear problem... 693 Fig. 1. Vibrations of beams in the case of linear (solid) and nonlinear (dashed) system (L=2; 3 loads): (a), (d) s=2.5m; (b), (e) s=3m; (c), (f) s=3.5m Fig. 2. The maximal response in the case of linear (solid) and nonlinear (dashed) model (L=2; 3 loads): (a), (d) s=2.5m; (b), (e) s=3m; (c), (f) s=3.5m summation kmin(n) = ωmin2 n −KC −1, kmax(n) = ωmax2 n−1 on the basis of the behaviour of the transformed solution, i.e. by taking into account the fact that the interval [ωmin,ωmax] includes all points ω important for the original solution.These can be recognized by the analysis of the supportof thepower spectrumfor each termappearing consecutively in the approximating sequence (Eqs. (4.6)). Figure 3 shows the solution for thedouble-beamsystemobtainedbyusing the criterion based on the power spectrum analysis in the case of a relatively big distance of 15m between 2 loads. One can see that the linear part of the solution clearly shows two separated excitations, whereas thenonlinearpart doesnot reflect our expectations.One couldpresumethat the second load also affects the nonlinear solution, especially when its distance from the first excitation is big enough topossiblyneglect the combined accumulation of the loads.Figure 4 presents the solution for the same system of parameters but with an increased range of summation: −818