. This workstudies thedynamics ofaonedimensional elastic barwith random elastic modulus andpre- scribed boundary conditions, say,ﬁxedatoneend, andattached toalumped mass andtwosprings (one linear andanother nonlinear) ontheother extreme. Thesystem analysis assumes thattheelastic mod- ulus hasgamma probability distrib ution anduses Monte Carlo simulations tocompute thepropag ation ofuncertainty inthiscontinuous–discrete system. After describing thedeterministic andthestochastic modeling ofthesystem, some conﬁgurations ofthemodel areanalyzed inorder tocharacterize theeffect ofthelumped mass intheoverall behaviorofthisdynamical system. Mecánica Computacional Vol XXXI, págs. 2673-2683 (artículo completo) Alberto Cardona, Paul H. Kohan, Ricardo D. Quinteros, Mario A. Storti (Eds.) Salta, Argentina, 13-16 Noviembre 2012 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 1 INTRODUCTION The dynamics of a mechanical system depends on some parameters such as physical and ge- ometrical properties, constraints, external and internal loading, initial and boundary conditions. Most of the theoretical models used to describe the behavior of a mechanical system assume nominal values for these parameters, such that the model gives one response for a given partic- ular input. In this case the system is deterministic and its behavior is described by a single set of differential equations. However, in real systems they do not have a ﬁxed value since they are subjected to uncertainties of measurement, imperfections in manufacturing processes, change of properties, etc. This variability in the set of system parameters leads to a large number of possible system responses for a given particular input. Now the system is stochastic and there is a family of differential equations sets (one for each realization of the random system) associated to it. This work aims to study the propagation of uncertainty in the dynamics of a nonlinear con- tinuos random system with a discrete element attached to it. In this sense, this work considers a one dimensional elastic bar, with random elastic modulus, ﬁxed on the left extreme and with a lumped mass and two springs (one linear and another nonlinear) on the right extreme (ﬁxed- mass-spring bar). This paper is organized as follows. In section 2 is presented the deterministic modeling of the problem, the discretization procedure and the algorithm used to solve the equation of interest. The stochastic modeling of the problem is shown in section 3, as well as the construction of a probability distribution for the elastic modulus, using the maximum entropy principle, and a brief discussion on the Monte Carlo method. In section 4, some conﬁgurations of the model are analyzed in order to characterize the effect of lumped mass in the system dynamical behavior. Finally, in section 5, the main conclusions are emphasized and some directions for future work outlined. 2 DETERMINISTIC APPROACH The continuous system of interest is the one-dimensional ﬁxed-mass-spring bar shown in Figure 1. xu(x, t) Lk kNLm Figure 1: Sketch of a bar ﬁxed at one and attached to two springs and a lumped mass on the other extreme.A. BARBOSA DA CUNHA JR, R. SAMPAIO 2674 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 2.1 Strong Formulation The displacement of this system uevolves according to the following partial differential equation A@2u @t2(x;t) +c@u @t(x;t) =@ @x EA@u @x(x;t) +f(x;t); (1) which is valid for 0<x<L and0<t<T , beingLthe bar unstretched length and Ta ﬁnite instant of time. In this equation is the mass density, Eis the elastic modulus, Ais the circular cross section area, cis the damping coefﬁcient, and f(x;t)is an external force depending on positionxand instantt. The left side of the bar is ﬁxed at a rigid wall while the right side is attached to a lumped massmand two springs ﬁxed to a rigid wall. The ﬁrst spring (of stiffness k) is linear and exerts a restoring force proportional to the stretching on the bar. The second spring (of stiffness kNL) is nonlinear and its restoring force is proportional to the cube of the stretching. The force which the lumped mass exerts on the bar is proportional to acceleration. These boundary conditions read as u(0;t) = 0 andEA@u @x(L;t) =ku(L;t)kNL u(L;t)3m@2u @t2(L;t):(2) Initially, any point xof the bar presents displacement and a velocity respectively equal to u(x;0) =u0(x) and@u @t(x;0) = _u0(x); (3) for0xL. In these equations u0and_u0are given functions of position x. Moreover, it is noteworthy that uis assumed to be as regular as needed for the initial– boundary value problem of Eqs.(1), (2), and (3) to be well posed. 2.2 Variational Formulation LetUtbe the class of (time dependent) basis functions and Wbe the class of weight func- tions. These sets are chosen as the space of functions with square integrable spatial derivative, which satisfy the essential boundary condition deﬁned by Eq.(2). The variational formulation of the problem under study says that one wants to ﬁnd u2U t that satisfy, for all w2W , the weak equation of motion given by M(u;w) +C( _u;w) +K(u;w) =F(w) +FNL(u;w); (4) whereMis the mass operator, Cis the damping operator, Kis the stiffness operator, Fis the ex- ternal force operator, and FNLis the nonlinear force operator. These operators are, respectively, deﬁned as M(u;w) =ZL 0 Au(x;t)w(x) dx+mu(L;t)w(L); (5) C( _u;w) =ZL 0c_u(x;t)w(x)dx; (6)Mecánica Computacional Vol XXXI, págs. 2673-2683 (2012) 2675 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar K(u;w) =ZL 0 EAu0(x;t)w0(x) dx+ku(L;t)w(L); (7) F(w) =ZL 0f(x;t)w(x)dx; (8) FNL(u;w) =kNL u(L;t)3w(L); (9) where _is an abbreviation for temporal derivative and0is an abbreviation for spatial derivative. The variational formulations for the initial conditions of Eq.(3), which are valid for all w2 W, are respectively given by fM(u(;0);w) =fM(u0;w); (10) and fM( _u(;0);w) =fM( _u0;w); (11) wherefMis the associated mass operator, deﬁned as fM(u;w) =ZL 0Au(x;t)w(x)dx: (12) 2.3 An Eigenvalue Problem Now consider the following generalized eigenvalue problem associated to Eq.(4), 2M( ;w) +K( ;w) = 0; (13) whereis a natural frequency and is an associated mode shape. In order to solve Eq.(13), the technique of separation of variables is employed, which leads to a Sturm-Liouville problem (Al Gwaiz, 2007), with denumerable number of solutions. There- fore, this generalized eigenvalue problem has a denumerable number of solutions, all of then such as the following eigenpair (2 n; n), wherenis then-th bar natural frequency and nis then-th bar mode shape. It is important to observe that, the eigenfunctions f ng+1 n=1span the space of functions which contains the solution of the Eq.(13) (Brezis, 2010). As can be seen in Hagedorn and DasGupta (2007), these eigenfunctions satisfy, for all m6=n, the orthogonality relations given by M( n; m) = 0; (14) and K( n; m) = 0; (15) which made then good choices for the basis function when a weighted residual procedure (Fin- layson and Scriven, 1966) is used to approximate the solution of a nonlinear variational equa- tion, such as Eq.(4).A. BARBOSA DA CUNHA JR, R. SAMPAIO 2676 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 2.4 Mode Shapes and Natural Frequencies According to Blevins (1993), a ﬁxed-mass-spring bar has its natural frequencies and the corresponding orthogonal modes shape given by n=nc L; (16) and n(x) = sin nx L ; (17) where c=p E= is the wave speed, and the nare the solutions of cot (n) +kL AE1 nm AL n= 0: (18) The ﬁrst six orthogonal modes shape of the ﬁxed-mass-spring bar with m= 1:5kg, whose the other parameters are presented in the beginning of section 4, are illustrated in Figure 2. In this ﬁgure each sub-caption indicates the approximated natural frequency associated with the corresponding mode. 00.20.40.60.8 100.20.40.60.81 x (m)φ1(x) (m) (a)10:07105rad=s 00.20.40.60.8 1−1−0.500.51 x (m)φ2(x) (m) (b)20:22105rad=s 00.20.40.60.8 1−1−0.500.51 x (m)φ3(x) (m) (c)30:37105rad=s 00.20.40.60.8 1−1−0.500.51 x (m)φ4(x) (m) (d)40:52105rad=s 00.20.40.60.8 1−1−0.500.51 x (m)φ5(x) (m) (e)50:67105rad=s 00.20.40.60.8 1−1−0.500.51 x (m)φ6(x) (m) (f)60:83105rad=s Figure 2: The ﬁrst six orthogonal modes shape and the corresponding (approximated) natural frequencies of a ﬁxed-mass-spring bar with m= 1:5kg. 2.5 Galerkin Formulation In order of approximate the solution of Eqs.(4), (10) and (11) the Galerkin method (Hughes, 2000) is employed. Therefore, the displacement ﬁeld uis approximated by a linear combinationMecánica Computacional Vol XXXI, págs. 2673-2683 (2012) 2677 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar of the form uN(x;t) =NX n=1un(t) n(x); (19) where the basis functions nare the orthogonal modes shape of the ﬁxed-mass-spring bar, exempliﬁed in the end of section 2.4, and the coefﬁcients unare time-dependent functions. For a reason that will be clear soon, deﬁne u(t)ofRNas the vector in which the n-th component is un(t). SinceuNis not a solution of Eq.(4), when the ﬁeld uis approximated by uNa residual function is obtained. This residual function is orthogonally projected into the vector space spanned by the functions f ngN n=1in order to minimize the error incurred by the approximation (Hughes, 2000). This procedure results in the following NNset of nonlinear ordinary differential equations [M]u(t) + [C]_u(t) + [K]u(t) =f(t) +fNL u(t) ; (20) supplemented by the following pair of initial conditions u(0) = u0 and _ u(0) = _ u0: (21) where [M]is the mass matrix, [C]is the damping matrix, [K]is the stiffness matrix, and the upper dot again denotes the time derivative. Also, f(t),fNL u(t) ,u0, and _ u0are vectors of RN, which respectively represent the external force, the nonlinear force, the initial position, and the initial velocity. The initial value problem of Eqs.(20) and (21) has its solution approximated by Newmark method (Newmark, 1959). The reader interested in details about this integration scheme is encouraged to see Hughes (2000). 3 STOCHASTIC APPROACH 3.1 Probabilistic Model Consider a probability space ( ;A;P), where is sample space, Ais a-ﬁeld over and Pis a probability measure. In this probabilistic space, the elastic modulus is assumed to be a random variable E: !Rthat associates to each event !2 a real number E(!). Con- sequently, the displacement of the bar is the random ﬁeld U:[0;L][0;T] !R, which evolves according the following stochastic partial differential equation A@2U @t2(x;t;! ) +c@U @t(x;t;! ) =@ @x E(!)A@U @x(x;t;! ) +f(x;t); (22) being the partial derivatives now deﬁned in the mean square sense (Papoulis and Pillai, 2002). This problem has boundary and initial conditions similar to those deﬁned in Eqs.(2) and (3), by changinguforUonly. 3.2 Elastic Modulus Distribution The elastic modulus cannot be negative, so it is reasonable to assume the support of random variableEas the interval (0;+1). Therefore, the probability density function (PDF) of Eis aA. BARBOSA DA CUNHA JR, R. SAMPAIO 2678 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar nonnegative function pE:(0;+1)!R, which respects the following normalization condition Z+1 0pE()d= 1: (23) Also, the mean value of Eis known real number E, i.e., E[E] =E; (24) where the expected value operator of Eis deﬁned as E[E] =Z+1 0E()pE()d (25) Finally, one also wants that Ehas ﬁnite variance, i.e., Eh (EE)2i <+1; (26) which is possible (Soize, 2000), for example, if E ln (E) <+1: (27) Following the suggestion of Soize (2000), the maximum entropy principle (Shannon, 1948; Jaynes, 1957a,b) is employed in order to consistently specify pE. This methodology chooses forEthe PDF which maximizes the differential entropy function, deﬁned by S[pE] =Z+1 0pE() ln pE() d; (28) subjected to (23), (24), and (27), the restrictions that effectively deﬁne the known information aboutE. Respecting the constraints imposed by (23), (24), and (27), the PDF that maximizes Eq.(28) is given by pE() =1(0;+1)1 E 1 2 E!0 @1 2 E1 A 1 (1=2 E) E0 @1 2 E11 A exp  2 EE! ; (29) where 1(0;+1)denotes the indicator function of the interval (0;+1),Eis the dispersion factor ofE, and indicates the gamma function. This PDF is a gamma distribution. 3.3 Stochastic Solver: Monte Carlo Method Uncertainty propagation in the stochastic dynamics of the continuous–discrete system under study is computed by Monte Carlo (MC) method (Metropolis and Ulam, 1949). This stochastic solver uses a Mersenne twister pseudorandom number generator (Matsumoto and Nishimura, 1998), to obtain many realizations of the random variable E. Each one of these realizations deﬁnes a new Eq.(4), so that a new variational problem is obtained. After that, these new variational problems are solved deterministically, such as in section 2.5. All the MC simulations reported in this work use 45samples to access the random system. Further details about MC method can be seen in Liu (2001); Shonkwiler and Mendivil (2009); Robert and Casella (2010).Mecánica Computacional Vol XXXI, págs. 2673-2683 (2012) 2679 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 4 NUMERICAL EXPERIMENTS The numerical experiments presented in this section adopt the following deterministic pa- rameters for the studied system: = 7900kg=m3,c= 10kNs=m ,A= 625 mm2, k= 650N=m ,kNL= 6501013N=m3,L= 1m, andT= 8ms. Besides that, four values for the lumped mass are considered: m= 1:5;7:5;15;and75kg. The random variable E, is characterized by E= 203GPa andE= 10% . The initial conditions for displacement and velocity are respectively given by u0= 1 3(x) + 2x; and _u0= 0; (30) where 1= 0:1mm and 2= 0:5103. Note that u0reaches the maximum value at x=L. This function is used to “activate" the spring cubic nonlinearity, which depends on the displacement at x=L. The time-dependent external force acting on the system has the form of a sine wave with circular frequency equal to the ﬁrst natural frequency f(x;t) = 1(x) sin (1t); (31) where the external force amplitude is = 1N. 0 2 4 6 8−600−400−2000200400600 t(ms)displacement at x=L(µm) mean value 98% prob. (a)AL=m10 0 2 4 6 8−600−400−2000200400600 t(ms)displacement at x=L(µm) mean value 98% prob. (b)AL=m2 0 2 4 6 8−600−400−2000200400600 t(ms)displacement at x=L(µm) mean value 98% prob. (c)AL=m1 0 2 4 6 8−600−400−2000200400600 t(ms)displacement at x=L(µm) mean value 98% prob. (d)AL=m0:2 Figure 3: This ﬁgure illustrates the mean value (blue line) and a 98% of probability interval of conﬁdence (grey shadow) for the random process U(L;;), for several values of the continuous–discrete mass ratio.A. BARBOSA DA CUNHA JR, R. SAMPAIO 2680 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar 4.1 Analysis of Random System Envelope of Reliability and Phase Space The mean value of U(L;;)and an envelope of reliability, wherein a realization of the stochastic system has 98% of probability of being contained, are shown, for different values of the continuous–discrete mass ratio AL=m , in Figure 3. By observing this ﬁgure one can note that, as the value of lumped mass increases, the decay of the system displacement amplitude de- creases signiﬁcantly. This indicates that this continuous–discrete system is not much inﬂuenced by damping for small values of the continuous–discrete mass ratio. The mean phase space of the ﬁxed-mass-spring bar at x=Lis shown, for different values of the continuous–discrete mass ratio, in Figure 4. The observation made in the previous paragraph can be conﬁrmed by analyzing this ﬁgure, since the system mean orbit tends from a stable focus to an ellipse as the continuous–discrete mass ratio decreases. In other words the limiting behavior of the system when AL=m!0+is a mass-spring system. −600−400−200 0200400600−10−5051015 displacement at x=L(µm)velocity at x=L(m/s ) (a)AL=m10 −600−400−200 0200400600−10−5051015 displacement at x=L(µm)velocity at x=L(m/s ) (b)AL=m2 −600−400−200 0200400600−10−5051015 displacement at x=L(µm)velocity at x=L(m/s ) (c)AL=m1 −600−400−200 0200400600−10−5051015 displacement at x=L(µm)velocity at x=L(m/s ) (d)AL=m0:2 Figure 4: This ﬁgure illustrates the mean value of the ﬁxed-mass-spring bar phase space at x=L, for several values of the continuous–discrete mass ratio. 4.2 Analysis of Random System PDF The difference between the system dynamical behavior, for different values of m, is even clearer if one looks to the PDF estimations of the (normalized) random variable U(L;T;), which are presented in Figure 5. For large values of the continuous–discrete mass ratio, the PDF ofU(L;T;)displays bimodal shape, which tends to a unimodal shape as the lumped mass grows, i.e., the continuous–discrete mass ratio decreases. Furthermore, it can be notedMecánica Computacional Vol XXXI, págs. 2673-2683 (2012) 2681 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar that whenAL=m1the greatest probability occurs around the mean value of U(L;T;). −4 −2 0 2 400.20.40.60.81 normalized displacement at x=Landt=Tprobability density function (a)AL=m10 −4 −2 0 2 400.20.40.60.81 normalized displacement at x=Landt=Tprobability density function (b)AL=m2 −4 −2 0 2 400.20.40.60.81 normalized displacement at x=Landt=Tprobability density function (c)AL=m1 −4 −2 0 2 400.20.40.60.81 normalized displacement at x=Landt=Tprobability density function (d)AL=m0:2 Figure 5: This ﬁgure illustrates estimations to the PDF of the (normalized) random variable U(L;T;), for several values of the continuous–discrete mass ratio. 5 CONCLUDING REMARKS A model to describe the dynamics of ﬁxed-mass-spring bar with a random elastic modulus is presented. The aleatory parameter is modeled as a random variable with gamma distribution, being the probability distribution of this parameter obtained by the principle of maximum en- tropy. The paper analyzes some conﬁgurations of the model to order to characterize the effect of the lumped mass in the overall behavior of this dynamical system. This analysis shows that the dynamics of the random system is signiﬁcantly altered when the values of the lumped mass are varied. 6 ACKNOWLEDGMENTS The authors are indebted to Brazilian Council for Scientiﬁc and Technological Development (CNPq), Coordination of Improvement of Higher Education Personnel (CAPES), and Founda- tion for Research Support in Rio de Janeiro State (FAPERJ) for the ﬁnancial support.A. BARBOSA DA CUNHA JR, R. SAMPAIO 2682 Copyright © 2012 Asociación Argentina de Mecánica Computacional http://www.amcaonline.org.ar