revue roumaine des sciences techniques

ROMANIAN JOURNAL OF TECHNICAL SCIENCES APPLIED MECHANICS formerly REVUE ROUMAINE DES SCIENCES TECHNIQUES SÉRIE DE MÉCANIQUE APPLIQUÉE

September − December 2013

Volume 58, Nº 3

CONTENTS RĂZVAN-VLAD VASIU, OCTAVIAN MELINTE, VICTOR VLĂDĂREANU, DAN DUMITRIU, On the response of the car from road disturbances ............................. IULIAN GIRIP, ŞTEFANIA DONESCU, MIHAELA POIENARIU, LIGIA MUNTEANU, On the behaviour of the hysteretic wire-rope isolators under random excitation...... EMIL-ALEXANDRU BRUJAN, YOICHIRO MATSUMOTO, Giant response of microbubble contrast agents oscillation in a high-intensity focused ultrasound field............................................................................................................................ EMIL-ALEXANDRU BRUJAN, Behaviour of plasmonic nanoparticle-generated cavitation bubbles....................................................................................................................... SANDA CLEJA-TIGOIU, NADIA ELENA STOICUTA, OLIMPIU STOICUTA, Numerical algorithms for solving the elasto-plastic problem with mixed hardening................................................................................................................... STEFAN STAICU, CONSTANTIN OCNARESCU, LIVIU UNGUREANU, Dynamics modelling of a parallel flight simulator..................................................................... ABDELMADJID TADJADIT, BOUALEM TILIOUINE, Formulation analytique des efforts de réduction au droit de l’interface fluide-structure des barrages rigides à fruits composés sous excitations sismiques...............................................................

Ro. J. Techn. Sci. − Appl. Mechanics, Volume 58, Nº 3, P., Bucharest, 2013

14p 14p 8p

10p 33p

13p 12p

ON THE RESPONSE OF THE CAR FROM ROAD DISTURBANCES RĂZVAN-VLAD VASIU 1, OCTAVIAN MELINTE 2, VICTOR VLĂDĂREANU 2, DAN DUMITRIU 2

Abstract. The tire/road contact dynamics is investigated in this paper by assuming a rocky road. The vehicle model to be used in simulation is a seven degree of freedom full-car model. The contact domain between the tire and the road has a shape defined by the Lamé curve, and it is identified by checking the minimum distance between bodies in contact. The nonlinear response from road disturbances is analyzed by assuming explicit relationships for contact and friction forces, respectively. Since the road is rocky, a bristle model was chosen to take into account the effect of the road irregularities on the tire’s action. Key words: full-car vehicle model, contact force, friction force, rocky road, contact dynamics.

1. INTRODUCTION

The impact and the frictional slip are met together when they simultaneously develop in a contact interface, such as the interface between the tire and the road. The distribution of forces can be manifested in many different forms and the friction model changes the characteristics of the vibro-impact phenomenon in terms of duration, dissipation of energy, accelerations and decelerations [1]-[4]. The irregularities of the road (bumps, corners, peaks, shallows, etc.) can not be taken into account without introducing the contact and friction forces in the contact interfaces. We adopt in this paper a continuum approach for the vibrocontact dynamics by admitting explicit relationship between contact force and deformation. In this continuous approach, no difference is made between impact and contact, therefore the methods of non-impact dynamics can be used to solve the problem [4]. The modeling of the contact friction in the interfaces has been performed in a series of relevant papers [5]-[10]. The contact force depends on the deformation and it is defined by the interference distance or penetration. Impacts between bodies are generally defined by the condition of impenetrability [11]. The continuum modeling of tire/road vibro-contact dynamics is developed in this paper for a rocky road and the contact domain is modeled by checking the minimum distance between bodies. The model takes as inputs the function of 1

Technical University of Cluj-Napoca, Faculty of Mechanics, Cluj-Napoca, Romania

2

Institute of Solid Mechanics of the Romanian Academy, Bucharest, Romania

Rom. J. Techn. Sci. − Appl. Mechanics, Vol. 58, N° 3, P. x−x, Bucharest, 2013

Răzvan-Vlad Vasiu, Octavian Melinte, Victor Vlădăreanu, Dan Dumitriu

2

bristle displacement and the vertical tire force and produces, as outputs, the distribution of contact pressure in the interfaces between the tire and the road in a short interval of time. The novelty of the paper in relation to mentioned publications from the literature is referring to the modeling of the tire/road contact for a rocky road. Modelings for such roads we do not found in the literature. For such roads only the bristle model is able to take into account the effect of the irregularities on the tire’s action. Our modeling is based on a seven degree of freedom full-car model. The idea of taking the Lamé curve for modeling the contact domain between the tire and the road is new. This idea together with the explicit relationships for contact and friction forces in the contact interfaces, form a new strategy introduced in this paper for investigation the tire behavior on rocky roads. 2. CONTACT BETWEEN THE TIRE AND THE ROAD

The contact between the tire and the road can be identified by checking the minimum distance between the tire and the road [5] 1  min  (r1 − r2 )T (r1 − r2 )  , f1 (r1 ) ≤ 0, f 2 (r2 ) ≤ 0 , 2 

(1)

where r1 and r2 are the position vectors of two points belonging to the tire and the road, respectively, and f1 and f2 are bounding surface constraints. The interference distance is defined as min(− d ),

d d f1 (r1 ) ≤ − e1 , f 2 (r2 ) ≤ − e2 . 2 2

(2)

The tire load and velocity generate forces at the interface between the tire and the road. These forces act in three directions, namely the contact force Fc = Fcz acting in the z direction, the longitudinal component of the friction force Ftx acting in the x direction, and the lateral component of the friction force Fty acting in the y direction, respectively (Fig. 1). Fig. 2 shows the force and moments vectors acting on an arbitrary point P of the tire located at the central intersection point of the tire and road surface plane. The undeformed tire has the radius R and rotates at constant angular speed Ω about the wheel axis. In the tire plane x and y axes are parallel to the longitudinal axis of the tire and the road surface normal, respectively. Fx is the longitudinal force acting in the direction of the forward motion of the tire in the plane of the road, Fy is the lateral force acting in the direction perpendicular to Fx in the plane of the road, and Fz is the vertical force acting in the direction of the vector normal to the road plane. Mx is the overturning moment acting about the axis parallel to Fx vector, My is the

On the response of the car from road disturbances

3

rolling resistance moment acting about the axis parallel to the Fy vector, and Mz is the aligning moment acting about the axis parallel to Fz vector.

Fig. 1 – The contact between the tire and the road.

The τ, | τ | < 90° is the slip angle between the travel direction and the direction in which the tire is oriented, γ is the camber angle, that means the inclination angle between the z-axis and the tire plane and ω is the angular velocity of the wheel measured at the centre of the wheel-tire assembly (the linear speed of the center of the wheel in the longitudinal direction is V ), and T is the axle torque. If we denote with px, py and pz, the longitudinal tire contact stress, the lateral tire contact stress and vertical tire contact pressure, the forces and moments, respectively in the point P ( x, y ) belonging to the contact area A are given by

∫

∫

∫

A

A

A

Fx ( P ) = − p x dA , Fy ( P ) = − p y dA , Fz ( P ) = − pz dA ,

∫

∫ pz ⋅ xdA ,

A

A

M x ( P) = − pz ⋅ ydA , M y (= P)

M z ( P= )

∫ ( px ⋅ y − p y ⋅ x)dA . A

(3)

(4)


4

From (1.3) and (1.4) we can determine the forces Fx ,Fy ,Fz and moments Mx ,My ,Mz acting at the center of the tire in the tire reference frame, by applying a rotation γ about x [12]  Fx ( P )    P)   Fy (=  F ( P)   z  0 0  M x  1    0 cos γ −= sin γ   M y      0 sin γ cos γ   M z 

0 0   Fx  1  0 cos γ − sin γ   F  ,   y   0 sin γ cos γ       Fz 

(5)

 M x ( P ) − Fy ( P ) R cos γ − Fz ( P ) R sin γ    M y ( P ) + Fx ( P ) R cos γ  .   M z ( P ) + Fx ( P ) R sin γ  

Fig. 2 – Forces and moments acting on the tire.

To shape of the unknown contact domain between the tire and the road is taken as a superellipse shape defined by the Lamé curve [8],[9] n

n

 x   y  1, n > 0 ,  a ( t )  +  b( t )  =    

(6)

where x and y define the envelope of the contact area, a is half of the contact length, b is half of the contact width (radii of the oval shape are depending of time), and n the power of the ellipsoid. For the contact area A area of (6), we have


5

n 1/ n

a

x 4b 1 −    A(t ) =  a   0

∫

 1 41−1/ n a (t )b(t ) πΓ 1 +   n , dx = 1 1 Γ +  2 n

(7)

where Γ is the Gamma function Γ( z ) = lim n→∞

n !n z , ( z ≠ 0, −1, −2,...) . z ( z + 1)...( z + n)

3. MODEL DESCRIPTION

The description of the car is carried out based on the papers [13]-[15]. An alternative approach can be found in [16] concerning only the vertical dynamics by using a 7 DOF car model. The dynamic contact between the wheels and the road cannot be studied by alone, but only considering the car as a unique dynamical system. The sprung mass m of the vehicle is defined as a rigid body that can translate and rotate, exhibiting angular oscillations around the mass centre G. The tires B j , j = 1, 2,3, 4 are represented as the tips of the tetrahedron Q. The trajectories of these tips are arbitrary in the space, but subjected to given geometric requirements. Let us consider in the following that the projections of the trajectories of the contact points tyre-road are parallel lines situated in the plane Q. In Fig. 3, H is the horizontal plane, O1A a segment on the intersection of planes Q and H, α is the angle between Q and H, n the unit vector defining the direction with the highest inclination of the plane Q, and P the projection of G in the plane Q, which describes the line O1P. A mobile reference frame T1 of an orthogonal set of vectors i, is attached to P ( i1 , i2 ⊂ Q , the sense of i1 is along the vehicle motion direction). A fixed reference frame T0 of an orthogonal set of vectors k, is attached to O1. The relation between i and k is given by  i1   − sin β cos α cos β sin α cos β   k1  i  = − cos β − cos α sin β − sin α sin β   k  ,  2   2  i3   0  − sin α cos α   k3 

(8)

where α, 0 ≤ α ≤ π / 2 , and β, −π ≤ β ≤ π , are the relative angular displacements of g having the directions i1 and i2. The roll and pinch angles are = n k2 cos α + k3 sin α and= i1 n cos β − k1 sin β , where n is the unit vector defining the direction with the highest inclination of the faces of tetrahedron Q with respect to H.


6

A new mobile frame T of base vectors u of the principal inertial directions of the sprung mass, is assigned fixed on the sprung mass (Fig. 4)  u1  u  =  2  u3 

1 0  θ2

0 −θ2   i1  1 θ  i2  , −θ1 1   i3 

(9)

γ sin θ , i2 = with i1 × u1 = u2 u1 sin ϕ , and θ, ϕ and ψ are the Euler × γ i1 sin ψ , γ ×= angles. The frame T rotates with respect to T1 to angular velocity ω defined as ω = ω1u1 + ω2u2 + ω3u3 ,

(10)

ω1 = ψ + ϕ = θ 1 , ω2 =θ 2 , ω3 =0 .

(11)

The motion equations of the sprung mass m with respect to T0 , with neglecting of the term ω× τω , are written under the form m

d2 dt 2

(OG ) = R , τc ω = M G ,

(12)

where τc is the moment of central principal inertia tensor, R the force acting on m, MG the moment with respect to G. The force R acting on the sprung mass is defined as (13) R =F + Fa + Fw + Pj .

∑ j

Here F = − mgk3 is the weight, = Fa + Fa3i3 , Fa2 = 0 the aerodynamic force applied in the pressure center P1 linked on the frame T1, on the position Fa1i1

vector ρ P = la1 i1 + la3i3 , la2 = 0 , depending only on the speed v of the vehicle, 1

Fw =

∑ Fwk ik

the wind force applied in the point P2 linked on the frame T1, on the

k

∑

position vector ρ P = lwk ik , and supposed to be a constant. The quantities Pj are 2

k

the constraint forces acting in points A j , j = 1, 2,3, 4 , that link the sprung mass m to four unsprung masses m j , j = 1, 2,3, 4 of the tires B j , j = 1, 2,3, 4 (Fig. 5). The variation ∆V j of the normal reaction forces in points A j is given by ∆V j =−

∑ c j ∆ j − ∑ k j ∆ j , j

where c j and k j , j = 1, 2,3, 4 , are stiffness and

j

damping coefficients of the suspensions associated to the static equilibrium point

7


∆j = 0 of A j . The quantities ∆ j represent the variation of the distance between

the mass centers of the tires B j and the points A j ∆ j = z − z j + αb j − βa j , j = 1, 2,3, 4 .

(14)

Fig. 3 – Geometry of the problem.

Fig. 4 – The Euler angles.

The initial conditions of the system motion are t = 0 , v = 0 , Fa = 0 , Pj = Pj0 , ϕ = ϕ0 , ψ = ψ 0 ,

θ = θ0 , ϕ0 + ψ 0 = α 0 , θ0 =β0 .

(15)

From (3.6) we have at t = 0 R0 =F + Fw + ∑ Pj0 =0 . j

(16)


8

Fig. 5 – The full car model with seven degrees of freedom.

The suspension of the vehicle with two symmetric axes are modeled by the following seven equations of motion for m and m j , j = 1, 2,3, 4 , on the direction i3 , and the second equation (12)

mz +

∑ c j ∆ j + ∑ k j ∆ j =Fa3 , j

m j z j + c′j z j + k ′j z j − c j ∆ j − k j ∆ j =g j , j = 1, 2,3, 4 ,

 + Aα

∑ b j c j ∆ j + ∑ b j k j ∆ j =M1 , j

 − Bβ

(18) (19)

j

∑ a j c j ∆ j − ∑ a j k j ∆ j = M 2 , j

(17)

j

j

(20)

where c′j and k ′j , j = 1, 2,3, 4 , are the stiffness and damping coefficients of the tires B j , M1 and M 2 are the moments acting on the mass m with respect to G,


9

and A, B the principal central inertial moments corresponding to the directions u1 and u2 . The forces g j , j = 1, 2,3, 4 , are given by = g j c′j s j + k ′j s j , j = 1, 2,3, 4 ,

(21)

with s j (t ) , j = 1, 2,3, 4 , the bristle displacement functions. The boundary conditions attached to the tires can be written as M j = z0 F fj , Fshj = Fcj − Fz , j = 1, 2,3, 4 ,

(22)

where M j are the bending moment of B j corresponding to the contact domain, Fz is the vertical load acting on B j , and Fsh is the shearing force in B j ,

respectively. The moments M j are equal to the couple moments created by the friction force F fj , j = 1, 2,3, 4 in B j . The friction force F fj is given by [18],[19] t  k fj s j (t0 ) + vtj dt , | s j |< s j max ,  t0 Ftj =  | Fcj | vtj  , otherwise,  µ k fj | vtj | 

∫

(23)

where k fj is the bristle stiffness, s j (t ) the bristle displacement functions, j = 1, 2,3, 4 , t0 is the start time of the last sticking at that contact point, vtj the

relative tangential velocity in B j and parameter s j max is the maximum allowable deflection of the bristle. The contact force Fcj are the contact forces in B j , given by [17] n p q Fcj = k δ j + bδ j δ j , j = 1, 2,3, 4 ,

(24)

with n j , p j , q j , j = 1, 2,3, 4 , are constants, the coefficient k depends on the tire’s material and the geometric properties of the tire, and b is defined with respect to the coefficient of restitution 0 ≤ e ≤ 1 , and the Ftj are friction forces in B j respectively, given by (23). In the following we consider that the stiffness and damping coefficients c′j and k ′j , j = 1, 2,3, 4 , of the tires are constants. The suspension force is considered to vary continuously between a soft suspension (to insulate against road disturbances), and a hard suspension (to insulate against load disturbances)


10

k2 (t ) . c= k1 (t ) , k= 1 c= 2 c1 (t ) , c= 3 c= 4 c2 (t ) , k= 1 k= 2 3 k= 4

(25)

csh , if csh < c s (hard suspension),  = c1,2  cs , if css < c s < csh ,  s s  cs , if cs < c s (soft suspension).

(26)

k sh , if k sh < k s (hard suspension),  = k1,2  k s , if k ss < k s < ksh ,  s s  k s , if k s < k s (soft suspension).

(27)

4. RESULTS

We select the following parameters for the model: m = 1500kg , = 40kg , ki′ = 250kN/m , ci′ = 0.6kN/m , i = 1, 2,3, 4 , m = = 40kg , m= 3 m 4 1 m 2 k h = 20kN/m , k s = 12kN/m , c h = 5kNs/m , c s = 0.5kNs/m . The identification of the contact patches in the interval T0 is performed by checking the minimum distance between bodies according to (1) and (2). As results, for a speed of the vehicle of 10 m/s, a number of more than 500 contact points were detected in the given interval of time in t ∈ T0 = [0; 5sec] . By defining a contact patch consisting from a minimum nearest 5 contact points, a number of more than 100 contact patches were identified. The functions of bristle-displacement s j (t ) , j = 1, 2,3, 4 ,

are represented in Fig. 6.

Fig. 6 – Functions of bristle-displacement s j (t ) .


11

Table 1 shows the characteristics of eight contact patches identified in T0 for four values of the vertical loads, 2000, 3000, 5000 and 8000 N, respectively. Table 1 Dimensions of the contact patches in the interval T0 Tire

Vertical load Fz [N]

a [mm]

b [mm]

n

B1

2000

35 38.3 45.2 50.2 67.8 72.9 90.1 91.8 35.3 38.1 45.0 51.0 68.0 73.6 89.5 92.5 35.8 38.5 49.2 61.4 52,7 63.4 97.8 101.0 35.9 38.5 53,1 63.8 71.9 79.2 98.3 100.1

46 51.2 53.9 54.8 63.4 65.6 71.8 74.7 45.6 51.5 53.7 54.4 63.2 65.5 72.3 74.9 45.9 51.7 56.8 64.5 57.0 61.2 75.1 76.4 45.8 51.8 57.1 61.4 68.0 67.6 75.0 76.3

2.02 2.18 2.41 2.55 2.88 2.92 3.18 3.95 2.02 2.17 2.40 2.53 2.86 2.91 3.15 3.93 2.03 2.16 2.42 2.54 2.67 2.71 4.03 4.13 2.03 2.15 2.66 2.70 2.78 3.00 4.04 4.11

3000 5000 8000 B2

2000 3000 5000 8000

B3

2000 3000 5000 8000

B4

2000 3000 5000 8000


12

Fig. 7 – The maximum contact pressure in the first contact patch.

The maximum value of the contact pressures for the first contact patch, with respect to a / a0 and b / b0 , are plotted in Fig. 7. The a0 and b0 are the reference radii of the oval shape of the contact domain. The time variation of the roll and pinch angles α and β, are shown in Fig. 8.

Fig. 8 – The time variation of the roll and pinch angles α and β.

5. CONCLUSIONS

Road profiles are central into investigation of the interaction between the tire and the road. This paper outlines a novel approach for the modeling of the tire/road contact dynamics by assuming a full-car vehicle model in interaction with a rocky road. The vehicle model to be used in simulation is a seven degree of freedom full-

13


car model. The response from road disturbances is analyzed by assuming an explicit relationship between the contact force and the deformation. An important aspect of this model is that the contact area increases with deformation and a plastic region can appear for larger indentation, i.e. the damping depends on the indentation. Another advantage of this model is that the contact force has no discontinuities at initial contact and separation, and it begins and finishes with the value of zero. For a rocky road, the bumpiness has to be taken into account. It is not unusual for automotive designers to test virtual models of cars on virtual models of bumpy roads. The bristle model for the friction force is chosen in this paper to take into account the effect of the road irregularities on the tire’s action. The suspension force is considered to vary continuously between a soft suspension (to insulate against road disturbances), and a hard suspension (to insulate against load disturbances). The contact domain between the tire and the road has a shape defined by the Lamé curve, and it is identified by checking the minimum distance between bodies in contact. Acknowledgements. This research was elaborated through the PN-II-PT-PCCA-2011-3.10190 Project of the National Authority for Scientific Research (ANCS, UEFISCDI), Romania. The authors acknowledge the similar and equal contributions to this article.

Received on May 7, 2013

REFERENCES 1. JALALI, H., AHMADIAN, H., POURAHMADIAN, F., Identification of micro-vibro-impacts at boundary condition of a nonlinear beam, Mechanical Systems and Signal Processing, 25, 3, pp. 1073-1085, 2011. 2. FERRI, A.A., Friction damping and isolation systems, Journal of mechanical Design, 117, B, pp. 196-206, 1995. 3. BERGER, E.J., Friction modeling for dynamic system simulation, Applied Mechanics Reviews, 55, 6, pp. 535-577, 2002. 4. GILARDI, G., SHARF, I., Literature survey of contact dynamics modelling, Mechanism and Machine Theory, 37, 10, pp. 1213-1239, 2002. 5. KARNOPP, D., Computer simulation of stick-slip friction in mechanical dynamic systems, Journal of Dynamic Systems, Measurement, and Control, 107, pp. 100-103, 1985. 6. MENQ, C.H., BIELAK, J., GRIFFIN, J.H., The influence of microslip on vibratory response, Part I: A new microslip model, Journal of Sound and Vibration, 107, 2, pp. 279-293, 1986. 7. BRIŞAN, C., VASIU, R.V., MUNTEANU, L., A modular road auto-generating algorithm for developing the road models for driving simulators, Transportation Research part C: Emerging Technologies, 26, pp. 269-284, 2013. 8. MUNTEANU, L., BRIŞAN, C., DUMITRIU, D., VASIU, R.V., CHIROIU, V., MELINTE, O., VLĂDĂREANU, V., On the modeling of the tire/road dynamic contact, Transportation Research part C: Emerging Technologies, 2013, in progress.


14

9. DUMITRIU, D., MUNTEANU, L., BRIŞAN, C., CHIROIU, V., VASIU, R.V., MELINTE, O., VLĂDĂREANU, V., On the continuum modeling of the tire/road dynamic contact, CMES: Computer Modeling in Engineering and Sciences, 94, 2, pp. 159-173, 2013. 10. MUNTEANU, L., BRIŞAN, C., CHIROIU, V., DONESCU, Şt., A 3D model for tire/road dynamic contact, Acta Technica Napocensis, Series: Applied Mathematics and mechanics, 55, 3, pp. 611-614, 2012. 11. KIM, S.W., Contact dynamics and force control of flexible multi-body systems, PhD Thesis, Department of Mechanical Engineering, McGill University, Montreal, 1999. 12. BRENDAN, J.-Y., CHAN, B. J-Y., Development of an off-road capable tire model for vehicle dynamics simulations, PhD Thesis, Virginia Polytechnic Institute and State University, 2008. 13. TRUICĂ, N., NICOLAE, V., A mathematical model for the study of oscillations of an automobile transmission and suspension, Symposium Mechanisms and Mechanical Transmissions, Reşiţa, Romania, pp. 786-804, 1976. 14. MAILAT, F., DONESCU, Şt., CHIROIU, V., MUNTEANU, L., On the nonlinear adaptive control of semi-active suspensions, AMSE: Advances in Modelling, Series C: Automatic control, theory and applications, 60, 5, pp. 15-28, 2005. 15. WANG, Fu-Chen, Design and synthesis of active and passive vehicle suspensions, PhD Thesis, University of Cambridge, Department of Engineering, 2001. 16. DUMITRIU, D., Car vertical dynamics 3D simulator using a 7 DOF model, Acta Technica Napocensis, Series: Applied Mathematics and Mechanics, 55, 3, pp. 647-650, 2012. 17. HUNT, K.H., CROSSLEY, F.R.E., Coefficient of restitution interpreted as damping in vibroimpact, Journal of Applied Mechanics, 42, 2, pp. 440-445, 1975. 18. HAESSIG, D,A., FRIEDLAND, B., On the modeling and simulation of friction, Journal of Dynamic Systems, Measurement, and Control, 113, 3, pp. 354-362, 1991. 19. MA, O., Contact dynamics modeling for the simulation of the space station manipulators handling payloads, IEEE International Conference on Robotics and Automation, vol. 2, pp. 1252-1258, 1995.

ON THE BEHAVIOUR OF THE HYSTERETIC WIRE-ROPE ISOLATORS UNDER RANDOM EXCITATION IULIAN GIRIP 1, ŞTEFANIA DONESCU 2, MIHAELA POIENARIU 3, LIGIA MUNTEANU1

Abstract. Difficulties in characterizing of behavior of the wire-rope isolators under random excitation is the motivation for investigation of the modified Bouc-Wen model in the spirit of Baber and Noori, in order to include the experimentally observed stiffness and strength degradation, and pinching, respectively. The prediction of the behavior of such structures with hysteretic shear behavior under random excitation is investigated in this paper. The model is verified for the case of harmonic excitation with a range of exciting frequencies and amplitude levels, on the base of reported data in literature. Key words: wire-rope isolator, Bouc-Wen model, stiffness and strength degradation, pinching.

1. INTRODUCTION

The wire-rope isolators are assemblies made of stranded wire ropes which are wrapped around a metallic or fibrous core (Fig. 1a). The section of a rope is shown in Fig. 1b. The rope is wound in the form of helix and held between metal retainers (Fig.1c) [1]-[3]. The diameter of the wire rope, the number of strands, the rope length, the cable twist or lay, the number of ropes per twist section and the fashion of metal retainer are the main parameters that characterize the assembly. The rubbing and sliding friction between the strands of wire rope define the damping and the deformation of isolators. For steady periodic excitation, the wire-rope isolators exhibit nonlinear hysteretic behavior due to the fact that the hysteretic loops depend on the vibration level, being almost independent of frequency [4]-[7]. Information from damage of the wire-rope isolators surveyed after realistic cyclic and arbitrary dynamic loadings and the experimentally observed characteristics, indicate that such structures exhibit a wide variety of hysteretic features including inelastic loaddisplacement law without distinct yield point, progressive loss of lateral stiffness in each loading cycle (stiffness degradation), degradation of strength when cyclically

1

Institute of Solid Mechanics of the Romanian Academy, Bucharest Technical University of Civil Engineering, Bucharest 3 University Spiru Haret, Bucharest 2


Iulian Girip, Ştefania Donescu, Mihaela Poienariu, Ligia Munteanu

2

load is done to the same displacement level (strength degradation) and pinching [8],[9]. The pinching behavior is due to slipping during force reversal.

Fig. 1 – a) The wire-rope geometry; b) Section A-A of the rope; c) Helical shape of isolators.

The smoothly varying Bow-Wen model is widely used for different loadings of hysteretic systems [10]-[18]. In this paper, a hysteretic model is developed and used to obtain the nonstationary response of the wire-rope isolators with hysteretic shear behavior under random excitation. Isolators without degradation effects are well described by the classical Bouc-Wen model. No such isolators there are in our attention due to the fact that most of the wire-rope isolators exhibit special characteristics such as the stiffness and strength degradation and pinching [3]. To predict such behavior, the classical Bouc-Wen model is modified in the spirit of Baber and Noori [19] in order to include the above mentioned observed characteristics. The new restoring force model is a modification of the Bouc-Wen model based on the endochronic theory of plasticity, and seems to be a versatile tool to model nonlinear, inelastic behavior, stiffness and strength degradation, the pinching and the memory. The law takes into account the observed dependence of the input and response at earlier time and memory. The model is verified for the case of harmonic excitation with a range of exciting frequencies and amplitude levels, on the base of reported data in literature [6],[17]. 2. MODEL DESCRIPTION

The helical wire-rope devices are used as vibration isolation systems. Fig. 2a shows the model of the wire-rope isolator arranged as a hanging shaking platform [6]. Two rigid plates are hung in parallel on a trestle, through frictionless hinges connected with four rigid steel tubes to form a double pendulum system. Additional guiding rollers prevent lateral motion of the system. The isolator is mounted and fixed with its aluminum retainer bars to the upper and lower plates for hysteretic shear behavior tests. The upper plate of the mass m on the top of the

3

On the behaviour of the hysteretic wire-rope isolators under random excitation

isolator is fixed horizontally to the trestle through a force transducer. The lower plate of mass m under the isolator is excited on one end and connected to the displacement transducer on the others. The lower mass is acted by the forcing force F(t). Fig. 2b presents the isolator model characterized by the linear spring k, the linear damping c and the hysteretic element, respectively. The appropriate requirements are taking into account to be compatible to the geometry shown in Fig. 2a. In this way, certain parameters (α,β,γ and n) are evaluated directly from the experiment. With the features added by the hysteretic element, we hope to be more closely to the observed hysteretic behavior of the wire-rope devices [24]. The relative displacement between plates is = u x2 − x1 . The motion equation written with respect to u is mu + cu + R (u , z , t ) = F (t ) ,

(1)

where cu is damping restoring force, with c the damping coefficient, and R (u , z , t ) is the non-damping restoring force R (u , z , t ) = αku + (1 − α)kz ,

(2)

composed by the linear restoring force αkz , and the hysteretic restoring force (1 − α)kz , where 0 < α < 1 is the rigidity ratio representing the relative participations of the linear and nonlinear terms. The hysteretic element is defined by a nonlinear first-order differential equation with a form similar to the endochronic theory of plasticity [18]. The function z(t) is the hysteretic auxiliary variable representing the hysteretic displacement function of the time history of u. It is related to u(t) through the constitutive law the force-displacement [18]

(

)

dz = η h( z ) A − ν (β sgn(u ) | z |n −1 z + γ | z |n ) , du

(3)

where h( z ) is the pinching function (for h = 1 the function is not pinch), A is a parameter that controls the tangent stiffness and ultimate hysteretic strength, β,γ,n are the hysteretic shape parameters and ν,η the strength and stiffness degradation functions (for ν = η = 1 the model is not degrading). These functions depend on the dissipated hysteretic energy. The law (3) extends the Bouc-Wen model by including the pinching function. By setting dz / du to zero in (3) and solving it for z, we obtain the ultimate hysteretic strength zu 1/ n

 A  zu =    ν(β + γ ) 

.

The pinching function h( z ) is taken under the form [16]

(4)


h( z ) = 1 − ζ1 exp(−( z sgn(u ) − zu ) 2 / ζ 22 ) ,

4

(5)

where ζ1 < 1 is a variable which controls the magnitude of initial drop in the slope dz / du , and ζ 2 a variable which controls the rate of change of the slope dz / du . By dividing the equation (1) by m, we obtain u + 2ζ 0 ω0u + αω02u + (1 − α)ω02 z = f (t ) ,

(6)

where f (t ) is the mass-normalized forcing function, ω0 = ki / m with ki the initial stiffness, ζ 0 = c / 2 ki m the linear stiffness. The identification procedure is based on the signification of each unknown parameter {α, A, β, γ, n, ζ1 , ζ 2 } and two unknown functions {ν, η} . We are interested to identify what parameters and functions of interest can be evaluated directly from the features of the experimental data. In the following we will see that α,β,γ and n can be directly evaluated from the experiment, i.e. the restoring force against displacement. The system properties are evaluated from the model. The first natural frequency is calculated as ki / m , where m is the estimated mass of the system and ki the initial stiffness. The value of the linear damping ratio ξ0 may be chosen within the range 0.01 and 0.05 [18]. The problem of the redundant parameters in (3) was discussed in [20]. It was proved that A is redundant and, as a consequence, it is assumed to be 1. The parameter α is calculated as the ration of the final tangent stiffness k f to initial stiffness ki α=

= ki

kf ki

,

dR dR = |z= u= 0 k , k f = |z = zu = αk . du du

(7) (8)

The equations for loading paths are obtained from (3) for ν = η = 1 and h = 1 dz = 1 − (β + γ ) z n , zu > 0 , du

(9)

dz = 1 − ( γ − β) z n , zu < 0 . du

(10)

and for unloading case

5


The parameters β and γ are determined from (3) for h = 1 and ν = η = 1 written as dz = 1 − (β + γ ) z n , z ≥ 0 , u ≥ 0 , du

(11)

dz = 1 − ( γ − β) z n , z ≥ 0 , u < 0 , du

(12)

dz = 1 + (−1) n +1 (β + γ ) z n , z < 0 , u < 0 , du

(13)

dz = 1 + (−1) n +1 ( γ − β) z n , z ≤ 0 , u ≥ 0 . du

(14)

Fig. 2 – The model of the wire-rope isolator arranged as a hanging shaking platform [6].

In what concerns the parameter n, we obtain from (4) and (9) for nondegrading case ν =1 n

 z  dz = 1−   . du  zu 

(15)


6

The shape parameters are chosen so that β + γ > 0 and γ − β ≤ 0 , with β > 0 in order to have a positive energy dissipation. In what concerns the stiffness and strength degradation functions ν and η, they are related to the dissipated hysteretic energy W [18] ν(W ) = 1 + δνW , η(W ) = 1 + δηW ,

(16)

where W

 = ω02 1 − 

kf   ki 

tf

∫ zudt ,

t0

0
0 , are unknown parameters. The stiffness degradation occurs when the elastic stiffness degrades with increasing ductility, as typically shown in Fig. 3 to left. This behavior occurs due to geometric effects. The strength degradation is described by reducing the capacity in the backbone curve, as shown in Fig. 3 to right. Fig. 4 shows the typical pinching behavior in the diagram dz / du against z / zu , where zu is the ultimate hysteretic strength given by (4). The pinching is the typical behavior of structures that buckle when subjected to compressive loads. This behavior usually is the result of cracks or slips. To identify the remaining parameters δν , δη , ζ1 and ζ 2 , that cannot be directly evaluated from the experiment, a genetic algorithm can be applied in the same manners as in [12], by using experimental data.

Fig. 3 – Left: Typical stiffness degradation behavior; Right: Typical strength degradation behavior [18].


7

Fig. 4 – Typical pinching behavior for h( z ) ≠ 1 [18].

3. RESULTS

The model is verified firstly for the case of harmonic excitation. The numerical simulation was conducted using experimental data reported in [3,4] and [12]. The excitation signal is supplied by sine harmonic excitation. The restoring hysteretic force, acceleration and displacements signals are measured for a number of 11 frequencies ranging from 5Hz to 50Hz, with at least 5 different amplitudes levels in each case, recorded in a synchronous manner on a tape recorder and observed with a digital signal analyzer. Fig. 5 shows the experimental loop of the hysteretic restoring force against the displacement, for the case of steady periodic excitation without filtering [6]. After filtering, the loop is shown in Fig. 6. The necessity of filtering is due to the non-synchronism which may appear between the restoring force and displacement during recording. This non-synchronism can distort the shape and area of the loop. Therefore, these signals must be simultaneously measured at each recording time. The parameters α,β,γ and n are evaluated directly from the experiment (Fig.6), while the parameters δν , δη , ζ1 and ζ 2 are obtained from an identification procedure based on the genetic algorithm. The parameters are summarized in Table 1. Table 1 The parameter values used in the numerical simulation for the case of harmonic excitation.

α

β

γ

n

ζ1

ζ2

δν

δη

0.26

0.33

0,88

1.89

0.42

0.33

0.11

0.10


8

Results for harmonic excitation agree very well with results presented in [6]. The obtained hysteretic loops show typical nonlinearity (for example Fig. 7 shows the restoring force against the displacement for 6 Hz).

Fig. 5 – Experimental hysteretic loop of restoring force against displacement without filtering [6].

Fig. 6 – Experimental hysteretic loop of restoring force against displacement loop with filtering [6].

9


Fig. 7 – Restoring force against displacement at 6 Hz.

Next, the case of the random excitation is analyzed. The excitation signal is represented as the sum of a number of sinusoids of random phase and frequency 0.25l + π /100 , l = 20,..., 40 [20]. The frequency range is approximately 5-10Hz. The amplitude of the excitation is chosen so that the strong nonlinearities and sliding of the isolator to be manifested. Methods for design of excitation signals can be found in [22, 23]. Fig. 8 shows the excitation signal computed by the periodic random (grey line) and by inverse repeated random (black line), respectively. As written by other researchers [14], the genetic algorithm gives best results with a few cycles of the non-linear response. The first results are presented in Figs. 9 and 10. The numerical simulation of the hysteretic loop of the excitation against the displacement computed for the periodic random excitation signal, is shown in Fig.9. Although this diagram looks like to experimental hysteretic loop without filtering shown in Fig. 5, it has nothing to do with the synchronism between the excitation force and displacement. Its character is done by the random feature of the loading. A similar behavior was obtained in [24] for friction pendulum systems under severe dynamic loading such as earthquakes.


10

Fig. 8 – Random excitation signal: periodic random (grey line) and inverse repeated random (black line).

The parameters are summarized in Table 2. Table 2 The parameter values used in the numerical simulation in the case of random excitation.

α

β

γ

n

ζ1

ζ2

δν

δη

0.47

0.87

0,90

2.15

0.94

0.72

0.21

0.22

Effects of degradation parameters upon the time variation of u is shown in Fig. 11 at 5 Hz. The numerical simulation of the evolution of the hysteretic displacement function z (t ) with respect to the relative displacement between plates, u is shown in Fig. 12, for 5Hz. Reducing the capacity in the backbone of the hysteretic curve and the increasing ductility in material leads to large displacement in both situations. Effect of pinching in the hysteretic loop of the restoring force against displacement 5 Hz is presented in Fig. 12. The double-valued of the restoring force for all times except when passes through the origin leads to pinching of the curve in origin.

11


Fig. 9 – Theoretical hysteretic loop of the excitation against the displacement.

Fig. 10 – Theoretical hysteretic loop z (t ) against the relative displacement u.


12

Fig. 11 – Left: Effect of the stiffness degradation δv parameter on u at 5 Hz; Right: Effect of the strength degradation parameter δη on u at 5 Hz.

Fig. 12 – Effect of pinching in origin of the hysteretic loop of restoring force against displacement at 5 Hz.

4. CONCLUSIONS

This paper tries to predict the behavior of the wire-rope isolators with hysteretic shear behavior under random excitation. To do this, the Bouc-Wen model is modified in the spirit of Baber and Noori in order to include the experimentally observed characteristics such as the stiffness and strength degradation and pinching, respectively. The hysteretic law is based on the

13


endochronic theory of plasticity. The model is verified for the case of harmonic excitation with a range of exciting frequencies and amplitude levels, on the base of reported data in literature. For random excitation, the wire-rope isolator was subjected to large deformations. It was observed that the isolator exhibits reducing of its capacity in the backbone of the hysteretic loop and the increasing ductility, respectively, as effects on the degradation. The hysteretic loop of restoring force against displacement is pinched at the origin for the random excitation. Acknowledgements. This research was elaborated through the PN-II-PT-PCCA-2011-3.10190 Project nr. 149/2012 of the National Authority for Scientific Research (ANCS, UEFISCDI), Romania. The authors acknowledge the similar and equal contributions to this article.

Received on April 15, 2013 REFERENCES 1. GILBERT, C., LEKUCH, H., Severe service requires strenuous shock/vibration testing, Military Electronics/Countermeasures, ICMD of North America Corp., 8, 11, pp 34-41, 1982. 2. TINKER, M.L., Modeling of nonlinear vibration isolators using advanced continuous simulation language (ACSL), report 1-10, 2000. 3. DAVIS, L.P., WILSON, J.F., JEWELL, R.E., RODEN, J.J., Hubble space telescope reaction wheel assembly vibration isolation system, In: Damping Proceedings 1986 (Las Vegas, Nevada), AFWAL-TR-86-3059, Air Force Wright Aeronautical Labs, Wright-Patterson AFB, Ohio, BA-1-BA-22, 1986. 4. WONG, C.W., NI, Y.Q., LAU, S.L., Steady-State Oscillation of Hysteretic Differential Model I: Response Analysis, Journal of Engineering Mechanics (ASCE), 120, 11, pp. 2271-2298, 1994. 5. WONG, C.W., NI, Y.Q., KO, J.M., Steady-State Oscillation of Hysteretic Differential Model II: Performance Analysis, Journal of Engineering Mechanics (ASCE), 120, 11, pp. 2299-2325, 1994. 6. KO, J.M., NI, Y.Q., TIAN, Q.L., Hysteretic Behavior and Empirical Modeling of a Wire-Cable Vibration Isolator, Int. J. Anal. Exp. Modal Anal., 7, 2, pp. 111-127, 1992. 7. NI, Y.Q., KO, J.M., WONG, C.W., Identification of non-linear hysteretic isolators from periodic vibration tests, Journal of Sound and Vibration, 217, 4, pp. 737-756, 1998. 8. DEMETRIADES, G.F., CONSTANTINOU, M.C., REINHORN, A.M., Study of wire rope systems for seismic protection of equipment in buildings, Engineering Structures, 15, 5, pp. 321-334, 1993. 9. FUJITA, T., SASAKI, Y., FUJIMOTO, S., TSURUYA, C., Seismic isolation of industrial facilities using lead-rubber bearings, JSME International Journal: Series 3 – Vibration, Control Engineering, Engineering for Industry, 33, 3, pp. 427-434, 1990. 10. BOUC, R., Forced Vibration of Mechanical System with Hysteresis, Proceedings of 4th Conference on Nonlinear Oscillation, Prague, 1967. 11. WEN, Y.K., Method for Random Vibration of Hysteretic Systems, Journal of the Engineering Mechanics Division (ASCE), 102, 2, pp. 249-263, 1976. 12. SIRETEANU, T., GIUCLEA, M., MITU, A.M., Identification of an extended Bouc-Wen model with application to seismic protection through hysteretic devices, Computational Mechanics, 45, 5, pp. 431-441, 2010. 13. GIRIP, I., ŞTIUCĂ, P., MUNTEANU, L., On a shape identification problem, Revue Roumaine des Sciences Techniques– Série de Mécanique Appliquée, 55, 3, pp. 211-218, 2010.


14

14. POIENARIU, M., IONESCU, M.F., GIRIP, I., MUNTEANU, L., CHIROIU, V., On the damped beams with hysteresis, World Journal of Mechanics (WJM), 1, 1, pp. 6-14, 2011. 15. MUNTEANU, L., MOŞNEGUŢU, V., GIRIP, I., POPESCU, M.A., On the influence of density on the wave propagation in fluids, Revue Roumaine des Sciences Techniques– Série de Mécanique Appliquée, 57, 1, pp. 63-70, 2012, ISSN 0035-4074. 16. GIRIP, I., MUNTEANU, L. GLIOZZI, A., Inverse problems in the hysteresis modeling, chap. 5 In: Inverse Problems and Computational Mechanics – vol.1 (eds. L.Marin, L.Munteanu, V.Chiroiu), pp. 125-146, Romanian Academy Publishing House, 2011. 17. WONG, C.W., KOAND, J.M., NI, Y.Q., Modeling the hysteresis curves of wire-cable isolators, Report Dept. of Civil and Structural Eng., Hong Kong Polytechnic, pp. 644-648, 1994. 18. FOLIENTE, G.C., Stochastic dynamic response of wood structural systems, PhD Thesis, Virginia Polytechnic Institute and State University, Blacksburg, 1993. 19. BABER, T.T., NOORI, M.N., Modelling General Hysteresis Behaviour and Random Vibration Application, Journal of Vibration, Acoustics, Stress, and Reliability in Design – Transactions of the ASME, 108, pp. 411-420, 1986. 20. MA, F., ZHANG, H., BOCKSTEDTE, A., FOLIENTE, G.C., PAEVERE, P., Parameter Analysis of the Differential Model of Hysterezis, J. of Appl. Mech. ASME, 71, pp. 342-349, 2004. 21. STAMMERS, C.W., SIRETEANU, T., Control of building seismic response by means of three semi-active friction dampers, Journal of Sound and Vibration, 273, 5, pp. 745-759, 2000. 22. JÁVORSZKY, G.B., BOYD, S., KOLLÁR, I., VANDENBERGHE, L., WU, S-P., Optimal excitation signal design for frequency domain systemn identification using semidefinite programming, Report Technical University of Budapest, 2000. 23. SCHOUKENS, J., SWEVERS, J., PINTELON, R., VAN DER AUWERAER, H., Excitation design for FRF measurements in the presence of non-lieanr distorsions, Mechanical Systems and Signal Processing, 18, 4, pp. 727-738, 2004. 24. DIMIZAS, Panos C., KOUMOUSIS, Vlasis K., System identification of non-linear hysteretic systems with application toi friction pendulum isolation systems, 5th GRACM International Congress on Computational Mechanics, Limassol, Cyprus, 29 June – 1 July, 2005.

GIANT RESPONSE OF MICROBUBBLE CONTRAST AGENTS OSCILLATION IN A HIGH-INTENSITY FOCUSED ULTRASOUND FIELD EMIL-ALEXANDRU BRUJAN *, YOICHIRO MATSUMOTO †

Abstract: The interaction of microbubble contrast agents (BR14) with high-intensity focused ultrasound is investigated by high-speed photography. The experiments were conducted at an ultrasound frequency of 1.08 MHz with a peak negative pressure of 1.8 MPa. The microbubble oscillations run on a giant response with a long expansion phase followed by a violent collapse. The implications in high-intensity focused ultrasound surgery are also discussed. Key words: cavitation, microbubble contrast agents, giant oscillations.

1. INTRODUCTION

There is an increasing interest in understanding the interaction of microbubble contrast agents with ultrasound waves. Hypotheses of potential benefits from these interactions suggested that microbubble contrast agents loaded with therapeutic substances could be targeted for destruction with ultrasound and thus enhance diffusion-mediated delivery by increasing localized concentration of the substances. The ability to increase tissue permeability and concomitantly augment localized drug concentrations through targeted microbubble destruction has fueled interest in developing efficient methods for delivering drugs and genetic material. Most investigators who have used ultrasound contrast agents for therapeutic applications worked with perfluorocarbon bubbles stabilized by an albumin or lipid shell [1]. The main advantage of this type of contrast agents is their fragility when exposed to ultrasound. Several studies have been conducted on the behaviour of microbubble contrast agents. Chomas et al. [2], de Jong et al. [3], Bouakaz et al. [4], and Doinikov et al [5] describe the different effects of ultrasound on microbubbles and demonstrate these effects by experiments using high-speed photography. Depending on the applied ultrasound amplitude and frequency, effects such as stable oscillation of microbubbles, inertial cavitation, *

University Politehnica Bucharest, Department of Hydraulics, Spl. Independentei 313, 060042 Bucharest, Romania † University of Tokyo, Department of Mechanical Engineering, 7-3-1 Hongo, Bunkyo, Tokyo 1138656, Japan Rom. J. Techn. Sci. − Appl. Mechanics, Vol. 58, N° 3, P. x−x, Bucharest, 2013

Emil-Alexandru Brujan, Yoichiro Matsumoto

2

coalescence, fragmentation, ultrasound induced damage of the shell causing gas to escape from microbubbles, and jetting are described. The effects of various factors including the ultrasound driving frequency, pulse length, peak negative pressure, bubble size and shell properties on the fragmentation of microbubbles were also investigated by Chomas et al. [5] and Bloch et al. [6]. Microjet formation during collapse of microbubbles in the vicinity of a boundary was experimentally observed by Prentice et al. [7]. They also noted that the jetting and the microbubble translation towards the boundary are dependent on the relative distance between microbubble and boundary. In this article we describe experimental investigations on the interaction of microbubble contrast agents with high-intensity focused ultrasound. Previous in vivo investigations have indicated that microbubble contrast agents increase the ablation efficiency of high-intensity focused ultrasound [8, 9]. We found that the microbubble oscillations run on a giant response with a long expansion phase followed by a violent collapse. In the expansion phase the maximum microbubble radius grows rapidly to several ten times their original size while the rapid compression yields fragmentation of the microbubble. 2. METHODS AND MATERIALS

All experimental investigations were performed with a third-generation contrast agent, BR14 (Bracco Diagnostics, Geneva, Switzerland). BR14 consists of perfluorocarbon-containing microbubbles stabilized by a phospholipid monolayer. The mean diameter of the microbubbles is between 2.5 to 3.0 µm. The diluted contrast agent solution is pumped through 100 μm syringe needle with a microinjector. High-intensity continuous wave ultrasound is generated by a 2 mm thick concave piezo-ceramic transducer with an outer diameter of 80 mm and a focal length of 80 mm (see Fig. 1). A linear amplifier (E&I 2100L) takes a sinusoidal signal from a function generator (NF Corporation WF 1974) and drives the transducer at a frequency of 1.08 MHz. In the degassed water (2 ppm in O2 concentration), ultrasound is gathered at the focus of the transducer where an aluminium wall is placed. Then, in the localized focal volume near the wall surface, high amplitude pressure fluctuations (peak negative pressure 1.8 MPa) are obtained. The peak negative pressure was measured by using a PVDF needle hydrophone (Imotec Messtechnik, type 80-0.5-4.0) with a rise time of 25 ns, and a sensitivity of 10.53 MPa/V (calibrated by the manufacturer up to a frequency of 10 MHz). We also note that, for this value, no cavitation activity was observed in the absence of BR14 microbubbles. The microbubble dynamics were recorded with a high-speed image converter camera (DRS Hadland, Imacon 200) using parallel illumination provided by a cw Nd:Yag laser (λ = 532 nm). Independent of the interframing time, the exposure time on the fluorescent screen of the image converter camera is 5 ns. The image on the fluorescent screen was recorded with an

3

Microbubble contrast agents in high-intensity focused ultrasound

intensified scan ICCD camera system (Photometrics AT200A) with a 1360×1024 pixel resolution per frame. The signal from the ICCD camera was then digitized with 12-bit resolution (256 grey levels) and passed to a computer. A Nikon lens (Ai-Micro Nikkor, 105 mm, F2.8S) with 10 close up rings enabled a field view of 2.56 mm × 2.09 mm.

Fig. 1 – Experimental arrangement for the investigation of the behaviour of a microbubble contrast agent in high-intensity focused ultrasound field.

3. RESULTS AND DISCUSSION

Figure 2 shows the oscillation of a microbubble situated at a distance of 770 μm from the wall. The microbubble motion is characterized by a long expansion phase, that last for at least 10 μs, followed by a rapid compression. During the collapse phase, when the kinetic energy of the bubble surpasses its surface energy, the microbubbles fragment into a number of smaller cavities. The images obtained in this figure do not clearly show the compression of the bubble prior to fragmentation due to the limited temporal resolution, although it clearly shows the resulting fragments. Up to three fragments can be seen in frame 8. No translation of the microbubble towards or away from the wall was observed. The experimental microbubble radius was calculated using the long and short radii of the elliptically shaped bubble, measured from the photographic frames, R = [(Rlong(Rshort)2]1/3. The maximum bubble radius is 132 ± 10 μm (frame 3), i.e. 50 times larger than the equilibrium radius of the microbubbles.


4

Fig. 2 - Photographic series of the interaction of a BR14 microbubble contrast agent with highintensity focused ultrasound. Frame interval 2 μs, frame size 2.56 mm × 2.09 mm.

Figure 3 illustrates the behaviour of three microbubbles situated at various distances from the wall. Microbubble “a” is situated at a distance of 725 μm from the wall while the distance between microbubble “b” and the wall is 170 μm. The third microbubble is attached to the wall. As in the previous case, a long expansion phase followed by a violent collapse are the main features of the microbubble dynamics. The maximum radius of the microbubble situated farthest from the rigid wall is 120 ± 10 μm (approximately 40 times larger than the equilibrium radius of the microbubbles). The microbubble “b” expands to a smaller value of the maximum radius, probably because it is outside the focal volume of the ultrasound. There is no migration of the microbubbles situated far from the wall and no clear evidence for the formation of a microjet during microbubble collapse was found in these cases. Only the attached bubble shows interaction with the rigid boundary, including jetting (frames 12 and 13). Jetting behaviour has been previously noticed in micrometer-sized shelled and unshelled bubbles situated close to a rigid boundary [7,10]. The oscillation time of both microbubbles situated far from the wall is about 24 μs. The microbubble attached to the wall has a larger oscillation time which is most likely a consequence of the interaction with a rigid boundary. It is well known from studies with laser generated bubbles that a rigid wall in close proximity to a bubble results in a lengthening of the oscillation period with about 20% [11]. The maximum collapse velocity was obtained in the case of the microbubble attached to the rigid wall and is about 30 m/s (mediated between frames 12 and 13). We note, however, that this value is only a lower estimate of the collapse velocity because the temporal resolution is not sufficient to resolve the final collapse phase where the microbubble achieves the minimum radius. Figure 4 shows the oscillation of two microbubbles that are attached to the rigid boundary. The maximum radius of these bubbles are 120 ± 10 µm (frame 9), and fragmentation of both microbubbles can be observed starting with frame 10. The shape of the lower bubble in frame 2 and the shape of the upper bubble in frame 6 suggest the formation of a liquid jet directed towards the boundary. We also note that no shock wave emission from the collapsing microbubbles have been

5


observed, probably as a consequence of microbubble fragmentation. Shock wave emission from un-encapsulated microbubbles were previously observed using the same experimental set-up and similar ultrasound parameters [12].

Fig. 3 - Photographic series of the interaction of three BR14 microbubble contrast agents with highintensity focused ultrasound. Frame interval 2 μs, frame size 2.56 mm × 2.09 mm.

Fig. 4 - Photographic series of the interaction of three BR14 microbubble contrast agents with highintensity focused ultrasound. Frame interval 2 μs, frame size 2.56 mm × 2.09 mm.


6

When unencapsulated bubbles of different sizes are subject to only a very small acoustic pressure amplitude, they will respond with small oscillations about their equilibrium radius. Going up in the driving amplitude will bring out the effects of nonlinearity manifesting themselves in a non-sinusoidal bubble oscillation and in the occurrence of several resonances. At even larger oscillation amplitudes, the bubble can be elongated from its equilibrium radius to arbitrarily large radius values and can be compressed down to near zero radius [13]. The present experimental investigations show a dramatic increase in the response of BR14 microbubbles in a high-intensity focused ultrasound field. Physically, this giant response comes about through the instability of the bubble, when the static pressure is overcome by the sound pressure during part of the sound cycle. Bloch et al. [6] also investigated the behaviour of BR14 microbubbles at values of the peak negative pressure similar to that used in the present experiment. The microbubbles were exposed to 2.25 MHz, two-cycle pulses at peak negative pressures ranging from 180 kPa to 1.42 MPa. The maximum bubble radius obtained at the peak negative pressure of 1.42 MPa was about 20 μm, i.e. six times smaller than in this experiment. We note, however, that the giant response strongly depends on frequency. The peak reduces with increasing frequency because the duration of the inertial instability gets smaller in proportion to the smaller period of the driving field [13-15]. Thus, low frequency more easily leads to a large bubble expansion and strong collapse, high frequencies to more collapses in a fixed amount of time [14,15]. Pressure-wave-excited contrast agent bubbles (Levovist) in the vicinity of rat kidney fibroblast cells were investigated by Wolfrum et al. [16]. Their results showed that even at moderate peak negative pressure amplitudes of less than 2 MPa, the contrast agent bubbles expand to more than 30 times their original radius. In a recent paper, McDannold and co-workers [17] investigated, in rabbits, whether the combination of ultrasound contrast microbubbles and high-intensity focused ultrasound could induce lesions in the brain. The lowest values for the production of lesions with the microbubbles were 1.2 W for 10-second pulsed sonications and 0.6 W for 20-second pulsed sonications, almost 12 times lower than a previously determined threshold for lesion creation. The 50% probability for inducing tissue necrosis required a temperature increase of 5.9°C, approximately half the 11.4°C threshold established in prior experiments. Interestingly, the temperature measured with the microbubbles present in the ultrasound focus appeared to be below the threshold for thermal damage, indicating that the lesions induced during high-intensity focused ultrasound with contrast microbubbles are caused by other cavitational bioeffects. The present results indicate that the giant response of microbubble oscillations in a high-intensity focused ultrasound field may be a potential source of tissue damage because tissue close to the microbubble wall can be compressed and stretched beyond the damage threshold.

7


4. CONCLUSION

The interaction of microbubble contrast agents (BR14) with high-intensity focused ultrasound is investigated by high-speed photography. The experiments were conducted at an ultrasound frequency of 1.08 MHz with a peak negative pressure of 1.8 MPa. The present results indicate that the microbubble oscillations run on a giant response with a long expansion phase followed by a violent collapse. In the expansion phase the maximum microbubble radius grows rapidly to several ten times their original size while the rapid compression yields fragmentation of the microbubble. The giant response of microbubble oscillations in a high-intensity focused ultrasound field may be a potential source of tissue damage because tissue close to the microbubble wall can be compressed and stretched beyond the damage threshold.

Acknowledgement. This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS – UEFISCDI, project number PN-II-ID-PCE-2011-3-0079.

Received on September 29, 2013

REFERENCES 1. BEKEREDJIAN, R., GRAYBURN, P.A., SHOHET, R.V., Use of ultrasound contrast agents for gene or drug delivery in cardiovascular medicine, J. Am. Coll. Cardiol., 45, pp. 329-335, 2005. 2. CHOMAS, J.E., DAYTON, P.A., MAY, D., ALLEN, J., KLIBANOV, A., FERRARA, K., Optical observation of contrast agent destruction, Appl. Phys. Lett., 77, pp. 1056-1058, 2005. 3. DE JONG, N., FRINKING, P.J.A., BOUKAZ, A., GOORDEN, M., SCHOURMANS, T., JINGPING, X., MASTIK, F., Optical imaging of contrast agent microbubbles in an ultrasound field with a 100-MHz camera, Ultrasound Med. Biol., 26, pp. 487-492, 2000. 4. BOUKAZ, A., VERSLUIS, M., DE JONG, N., High-speed optical obervations of contrast agent destruction, Ultrasound Med. Biol., 31, pp. 391-399, 2005. 5. CHOMAS, J.E., DAYTON, P.A., MAY, D., ALLEN, J., KLIBANOV, A., FERRRARA, K., Threshold of fragmentation for ultrasonic contrast agents, J. Biomed. Opt., 6, pp. 141-150, 2001. 6. BLOCH, S.H., WAN, M., DAYTON, P.A., FERRARA, K.W., Optical observation of lipid- and polymer-shelled ultrasound microbubble contrast agents, Appl. Phys. Lett., 84, pp. 631-633, 2004. 7. PRENTICE, P., CUSCHIERI, A., DHOLAKIA, K., PRAUSNITZ, M., CAMPBELL, P., Membrane disruption by optically controlled microbubble cavitation, Nature Phys., 1, pp. 107-110, 2005. 8. HANAJIRI, K., MARUYAMA, T., KANEKO, Y., MITSUI, H., WATANABE, S., SATA, M., NAGAI, R., KASHIMA, T., SHIBAHARA, J., OMATA, M., MATSUMOTO, Y., Microbubble-induced increase in ablation of liver tumors by high-intensity focused ultrasound, Hepatol. Res., 36, pp. 308-314, 2006.


8

9. YU, T., FAN, X., XIONG, S., HU, K., WANG, Z., Microbubbles assist goat liver ablation by high intensity focused ultrasound, Eur. Radiol., 16, pp. 1557-1563, 2006. 10. BRUJAN, E.A., The role of cavitation microjets in the therapeutic applications of ultrasound, Ultrasound Med. Biol., 30, pp. 381-387, 2004. 11. BRUJAN, E.A., NAHEN, K., SCHMIDT, P., VOGEL, A., Dynamics of laser-induced cavitation bubbles near elastic boundaries: Influence of the elastic modulus, J. Fluid Mech,. 433, pp. 283-314, 2001. 12. BRUJAN, E.A., IKEDA, T., MATSUMOTO, Y., On the pressure of cavitation bubbles, Exp. Therm. Fluid Sci., 32, pp. 1188-1191, 2008. 13. BRUJAN, E.A., The effect of polymer concentration on the non-linear oscillation of a bubble in a sound-irradiated liquid, J. Sound Vib., 173, pp. 329-342, 1994. 14. CRUM, L.A., GAITAN, D.F., Sonoluminescence and its application to medical ultrasound risk assessment, Proc. SPIE, 1161, pp. 125-135, 1989. 15. GAITAN, D.F., CRUM, L.A., CHURCH, C.C., ROY, R.A., Sonoluminescence and bubble dynamics for a single, stable, cavitation bubble, J. Acoust. Soc. Am., 91, pp. 3166-3183, 1992. 16. WOLFRUM, B., METTIN, R., KURZ, T., LAUTERBORN, W., Observations of pressure-waveexcited contrast agent bubblesin the vicinity of cells, Appl. Phys. Lett., 81, pp. 5060-5062, 2002. 17. McDANNOLD, N.J., VYKHODTSEVA, N.I., HYNYNEN, K., Microbubble contrast agent with focused ultrasound to create brain lesions at low power levels: MR imaging and histologic study in rabbits, Radiology, 241, pp. 95-106, 2006.

BEHAVIOUR OF PLASMONIC NANOPARTICLE-GENERATED CAVITATION BUBBLES EMIL-ALEXANDRU BRUJAN *

Abstract: The transient dynamics of cavitation nanobubbles generated through the irradiation of gold nanoparticles, via the particle's plasmonic resonance, very close to the threshold laser fluence for bubble formation, is investigated numerically. Both for short and long laser pulses, the maximum bubble radius is in the nanometer range. The maximum bubble radius is almost proportional to the square root of the laser fluence. For long laser pulses, and very close to the threshold fluence for bubble formation, the energy dependence of the bubble size is stronger. The oscillation time of cavitation nanobubbles is smaller than that predicted by the Rayleigh formula. Excellent agreement was found between previous experimental results and the present calculations. The results indicate that the mechanical effects associated with bubble formation are compatible with intracellular nanosurgery and are also of interest for a gentle transfection of genes into biological cells. Key words: bubble, nanoparticle, nonlinear dynamics, modeling and simulation.

1. INTRODUCTION

It was recently demonstrated that the mechanical properties of cavitation bubbles, generated through the irradiation of plasmonic gold nanoparticles, can be used to produce a surgical device. A cavitation bubble is induced when a plasmonic gold nanoparticle is overheated with a laser pulse. As a result, the nanoparticle evaporates a very thin volume of the surrounding medium, thus generating a cavitation bubble that grows and then collapses within a very short period of time. The fast expansion of the bubble generates a localized mechanical impact on the environment with only minor thermal contributions [1]. It was also showed that plasmonic nanoparticle-generated bubbles can be efficiently used at the cellular and tissue level in vitro and in vivo as therapeutic, delivery and theranostic agents [2,3]. Cavitation bubbles can also aid diagnostics, for instance during imaging, as they may act as scattering centers for the probe. Although the thermalization of nanoparticles has been extensively studied [4], the dynamics of plasmonic nanoparticle-generated bubbles are less explored. Several theoretical studies have been conducted to investigate the initial phase of bubble formation around laser-irradiated nanoparticles and microparticles [5,6]. A *

University Politehnica Bucharest, Department of Hydraulics, Spl. Independentei 313, 060042 Bucharest, Romania Rom. J. Techn. Sci. − Appl. Mechanics, Vol. 58, N° 3, P. x−x, Bucharest, 2013

Emil-Alexandru Brujan

2

recent publication [7] describes the formation of a cavitation bubble around a plasmonic nanoparticle in which an expression for the acoustic wave pressure as a function of laser pulse energy is given. In any applications of plasmonic nanoparticles it is important to understand and control the spatial and temporal behaviour of cavitation bubbles in the liquid around nanoparticles. In a previous study [8] we have described the dynamics of cavitation bubbles at values of the laser fluence much larger than the threshold value for bubble formation. The present study describes numerical investigations of the dynamics of cavitation bubbles generated around gold nanoparticles irradiated with laser pulses at fluence values very close to the threshold value for bubble formation. We used a spherical model of bubble dynamics valid to first-order in the bubble wall Mach number. The effects of nanoparticle radius, laser fluence, and pulse duration on the maximum radius and oscillation time of the nanobubbles are discussed. Model predictions are then compared with the available experimental results. 2. MATHEMATICAL FORMULATION

Consider a spherical gold nanoparticle of radius Rnp suspended in a liquid of infinite extent and irradiated by a laser pulse with duration τL. When a certain laser fluence threshold is reached, the temperature of the nanoparticle exceeds the boiling point of the host liquid, and a vapour layer is formed on the surface of the nanoparticle. The large pressure in the vapour layer surrounding the bubble leads to a very rapid expansion of the layer with subsequent formation of a cavitation bubble. A spherical model of bubble dynamics is used to calculate the temporal development of the bubble radius [9]. The model considers the compressibility of the liquid, surface tension and the liquid viscosity. The bubble dynamics is described by the equation

(

)

3 1  + 6 RRR   + 2 R 3 = RR + R 2 − R2 R H, 2 c∞

(1)

where R is the time-dependent bubble radius, dots denote a time derivative, c∞ is the undisturbed speed of sound in the medium and H is the enthalpy change between bubble wall and infinity, n ( p0 + B )  P + B   = H  ( n − 1) ρ0  p0 + B 

( n −1) / n

 − 1 .  

(2)

The constants B and n relate to the Tait equation of state for water that was used to derive Eq. (2). The pressure P at the bubble interface is

3

Behaviour of plasmonic nanoparticle-generated cavitation bubbles

3 3  2σ   R0 − Rnp  P=  p0 +  3 R0   R3 − Rnp  

κ

 2σ R  − − 4η .  R R 

(3)

Here R0 denotes the initial radius of the bubble, κ the ratio of the specific heat at constant pressure and volume, p0 the initial pressure inside the bubble, σ the surface tension and η the medium viscosity. Assuming that the excess light energy absorbed by the particle during the laser pulse is spent in boiling at the critical point of water, the initial radius of the bubble can be estimated as [7]   3  ( F − Fc ) σabs  3  R0  =   + Rnp  Ecl  4πρcl   

1/ 3

.

(4)

In the above equation, Ecl is the internal energy of the liquid at the critical point, ρcl is the critical density of the liquid, σabs represents the absorption cross section, F is the laser fluence, and Fc is the critical laser fluence required to heat the nanoparticle to the boiling temperature Tboil , which in the case of short laser pulses 2 , where χl is the thermal diffusivity of the liquid) is given by [7] ( χl τ L > Rnp ), the expression for the critical laser fluence

becomes [7] Fc = 4πRnp cl ρl χl τ LTboil / σabs .

(6)

In the numerical calculations we take B = 3049.13 bar and n = 7.15 in the equation of state of the liquid. Other pertinent numerical values are κ = 1.4, Ecl = 2000 kJ/kg, p0 = 101325 Pa, ρl = 998 kg/m3, ρcl = 322 kg/m3, and σ= 0.07 N/m, cnp = 0.13 kJ/(kg⋅K), cl = 4.18 kJ/(kg⋅K), ρnp = 19300 kg/m3, and χl =1.4⋅10-7 m2/s. The temperature and the pressure inside the initial bubble are those at the critical point of water (Tcl = 647.1 K, pcl = 22.06 MPa). 3. RESULTS AND DISCUSSION

Figure 1 illustrates the temporal variation of the radius of cavitation bubbles generated from a spherical gold nanoparticle with a radius of 20 nm. The most prominent feature of the transient cavitation bubbles is their small size and short


4

Fig. 1 – Variation with time of the radius of cavitation bubbles generated after laser irradiation of a spherical gold nanoparticle with radius Rnp = 20 nm, for various values of the ratio F/ Fc . (a) bubbles generated from short laser pulses, (b) bubbles generated from a long laser pulse with duration τL = 10 ns. Solid line: F / Fc = 1; dotted line: F / Fc = 4; dashed line with one point: F / Fc = 7; dashed line with two points: F/ Fc = 10. The inset shows the bubble radius as a function of time for F / Fc = 1.1.

lifetime. For example, the maximum radius, Rmax , of the bubble generated from a nanoparticle irradiated with a laser pulse of 10 ns at F/ Fc =10 amounts to only 395 nm in water, and will be even smaller in a viscoelastic medium such as the cytoplasm. Such small bubbles are similar to those generated by focusing femtosecond laser pulses in the bulk of a liquid [10]. They are, however, three order of magnitude smaller than the bubbles generated by nanosecond optical breakdown in water (Rmax > 0.5 mm for τ L = 6 ns) [11]. For the smallest value of the laser fluence used in the present calculations ( F/ Fc =1.1) the oscillation becomes strongly damped. For example, in the case of long laser pulses, the maximum bubble radius is 38.4 nm and only two oscillation cycles can be seen. For short laser pulses, the maximum bubble radius is smaller with a factor of 5 for F/ Fc =10 (Rmax = 85 nm), and with a factor of 2 for F/ Fc =1.1 (Rmax = 21.5 nm). This makes a dissection mechanism associated with bubble formation compatible with intracellular nanosurgery because such small bubbles are essential to dissect

5


or knock out subcellular structures with nanometer precision. They are also of great interest for a gentle transfection of genes and transfer of other substances into specific cell types [12]. On the other hand, the small size of these bubbles is consistent with the fact that the collective action of a large number of nanoparticles is required to produce cell lysis.

Fig. 2 – Maximum bubble radius as a function of the ratio, F/ Fc , for various values of the nanoparticle radius, Rnp . (a) bubbles generated from short laser pulses, (b) bubbles generated from a long laser pulse with duration τL = 10 ns. Dashed line: Rnp = 10 nm; dashed line with one point: Rnp = 20 nm; dashed line with two points: Rnp = 40 nm.

Figure 2 presents the maximal bubble radius as a function of the ratio F/ Fc. At equal laser fluence, the maximum radius of the cavitation bubble decreases with decreasing nanoparticle radius and pulse duration. For example, for short laser pulses and F/ Fc = 4, the maximum radius of the bubble generated from a nanoparticle with Rnp = 10 nm (Rmax = 20.3 nm) is with a factor of 6 smaller than for Rnp = 40 nm (Rmax = 121.2 nm). The corresponding values for a long laser pulse ( τ L = 10 ns) are Rmax = 162 nm, for Rnp = 10 nm, and Rmax = 274.5 nm, for Rnp = 40 nm, respectively. For F/ Fc >2, the Rmax (F/ Fc) dependence is generally not very strong, and even for the largest nanoparticle (Rnp = 40 nm) irradiated at the


6

largest laser fluence ( F = 10 Fc ), Rmax is much smaller than the size of a biological cell. For example, the largest value of the maximal bubble radius is 514 nm and was obtained for Rnp = 40 nm and long laser pulses. Bubbles will be even smaller in cells where they are confined by the cytoskeleton [13]. The scaling law for the bubble size is the same for both short and long laser pulses and all values of Rnp , namely, the maximum bubble radius is proportional to the square root of the laser fluence. This is in contrast to the case of large bubbles (Rmax > 1 mm), generated by nanosecond optical breakdown in water, where the maximum bubble radius is proportional to the cube root of the laser pulse energy [11]. For long laser pulses ( τ L = 10 ns), this scaling law applies, however, only for F/ Fc >2. Closer to the threshold, the energy dependence of the bubble size is stronger. A similar trend was experimentally observed by Siems et al. [14] for 60-nm gold nanoparticles irradiated with nanosecond laser pulses.

Fig. 3 – Oscillation time of a bubble as a function of F / Fc , for various values of Rnp . (a) bubbles generated from short laser pulses, (b) bubbles generated from a long laser pulse with duration τL = 10 ns. Dashed line: Rnp = 10 nm; dashed line with one point: Rnp = 20 nm; dashed line with two points: Rnp = 40 nm.

7


It is also interesting to note that the oscillation time of cavitation nanobubbles is smaller than that predicted by the Rayleigh formula [15] 1/ 2

R Tosc

ρ  = 1.83  l   p0 

Rmax .

(7)

Figure 3 illustrates the collapse time of the nanobubbles in comparison to the Rayleigh collapse time, both for the case of short (Fig. 3a) and long (Fig. 3b) lasers pulses. For a nanoparticle radius Rnp = 20 nm irradiated with a laser pulse with duration τL = 10 ns at a fluence F/ Fc = 4, the oscillation time of the resulting nanobubble is 12.6 ns, which is less than half the value predicted by the Rayleigh formula (31.3 ns). The corresponding value obtained in the case of short laser pulses is 1.34 ns, with a factor of 6 smaller than the value predicted by the Rayleigh formula (7.4 ns). The differences become larger with decreasinglaser fluence and nanoparticle size. For large bubbles, as those generated by nanosecond or picosecond laser pulses in the bulk of a liquid, the expansion and collapse of cavitation bubbles are highly symmetrical and Rmax can be calculated from Tosc using Eq. (7) [11]. However, for bubbles with maximum radii smaller than a few micrometers the Rayleigh equation is not exact because it neglects surface tension and the viscosity of the medium surrounding the bubble. Both factors produce a pressure scaling with 1/R which adds to the hydrostatic pressure. This finding leads to the conclusion that, in the case of nanobubbles, the use of the Rayleigh oscillation time to determine the bubble size, as done in previous studies [16], leads to an underestimation of the maximum bubble size.

Fig. 4 – Comparison between simulated (dashed line) and experimental (filled symbols) results for the case of spherical gold nanoparticles with radius of 39 nm, irradiated with a short laser pulse τL = 100 fs, in water, at a delay of maximum bubble radius of 650 ps.


8

A comparison between simulated and experimental results obtained by Kotaidis et al. [17], for the case of spherical gold nanoparticles with radius of 39 nm, irradiated with a short laser pulse ( τ L = 100 fs), in water, is shown in Fig. 4. Very good agreement was found between experimental and numerical results. This implies, of course, that the neglected factors have little influence on the bubble dynamics or at least that some of the effects neutralize each other. It is interesting to note here that this model assumes that all the excess energy beyond heating the nanoparticle to the critical point is spent on boiling and used for forming the initial bubble. In the case of nanoparticles irradiated not under thermal confinement conditions, a large volume around the particle might be heated. Since the thermal gradient is very smooth, a significant heated volume might not reach the critical point. This heat cannot be used for forming the initial bubble and is lost. We also note that far above the threshold irradiation, Fc , the bubble can also nucleate during the laser pulse, which might lead to thermal insulation between the particle and the fluid by the vapour inside the bubble during the laser pulse as well as to optical scattering of the incident light. In this case, again, not all the applied laser pulse energy can be used to vaporize the liquid surrounding the nanoparticle. Therefore, the present model can only be used as an upper estimation limit for the bubble size. On the other hand, the surface tension of water is decreasing strongly with temperature until it vanishes at the critical point (see, for example, IAPWS Release on Surface Tension of Ordinary Water Substance, available at www.iapws.org), which can lead to larger bubbles than assumed in the model. Despite the desirable theoretical refinements that would be necessary to provide a complete description of bubble motion, the present formulation still gives a reasonable quantitative picture of bubble dynamics, and thus provides an essential practical connection between theory and experiment. Recent years have seen a continuous rise of interest in nanosurgery on a cellular and sub-cellular level. One important application is the separation of individual cells or other small amounts of biomaterial from heterogeneous tissue samples for subsequent genomic or proteomic analysis. According to the present results, when suitable laser parameters are used, the resulting cavitation bubble is in the nanometer range, and surgery can be performed at any desired location within a cell or a very small organism. It is worth noting here that time-resolved investigations of the behaviour of laser-induced nanobubbles around gold nanoparticles within cells are not yet available. However, Dayton et al. [19] investigated the oscillations of bubbles with a radius of 1.5 μm that were phagocytosed by leukocytes. By means of steak photography and high-speed photography with 100 million frames/s they observed that phagocytosed bubbles expanded about 20 - 40% less than free. The difference is due to the viscoelastic properties of the cytoplasm and cytoskeleton. This however remains a potentially important point that will form the object of a separate investigation.

9


4. CONCLUSION

The dynamics of cavitation nanobubbles, generated after irradiation of spherical gold nanoparticles with laser pulses close to the threshold value for bubble formation, was investigated by numerical calculations. The significant parameters of this study are the nanoparticle radius, Rnp , the laser fluence normalized by the threshold value for bubble formation, F/ Fc , and the laser pulse duration. The maximum bubble radius increases with increasing nanoparticle radius, pulse duration and laser fluence. For F/ Fc < 10, both for short and long laser pulses, the maximum bubble radius is in the nanometer range. For the smallest value of the laser fluence used in the present calculations (F/ Fc = 1.1) the oscillation becomes strongly damped. This makes a dissection mechanism associated with bubble formation compatible with intracellular nanosurgery and is also of great interest for a gentle transfection of genes and into specific cell types. For laser fluences larger than the threshold value for bubble formation, the maximum bubble radius scales with the square root of laser fluence. This is in contrast to the case of large bubbles (Rmax > 1 mm) generated by nanosecond optical breakdown in water where the maximum bubble radius is proportional to the cube root of the laser pulse energy. For long laser pulses and closer to the threshold, the energy dependence of the bubble size is stronger. The oscillation time of cavitation nanobubbles is smaller than that predicted by the Rayleigh formula. Thus, in the case of nanobubbles, the use of the Rayleigh oscillation time leads to an underestimation of the maximum bubble size. The Rayleigh equation cannot describe the oscillation time of bubbles smaller than a few micrometers because it neglects surface tension and viscosity that produce a pressure scaling with 1/R which adds to the hydrostatic pressure. The experimental data agree very well with the present numerical results.

Acknowledgement. This work was supported by a grant of the Romanian National Authority for Scientific Research, CNCS – UEFISCDI, project number PN-II-ID-PCE-2011-3-0079.

Received on September 29, 2013

REFERENCES 1. LAPOTKO, D., Pulsed photothermal heating of the media during bubble generation around gold nanoparticles, J. Heat Mass Transf., 52, pp. 1540-1543, 2009. 2. KHLEBTSOV, B.N., ZHAROV, V.P., MELNIKOV, A.G., TUCHIN, V.V., KHLEBTSOV, N.G., Optical amplification of photothermal therapy with gold nanoparticles and nanoclusters, Nanotechnology, 17, pp. 5167-5179, 2006. 3. LUKIANOVA -HLEB, E.Y., SAMANIEGO, A.P., WEN, J., METELITSA, L.S., CHANG, C.C., LAPOTKO, D.O., Selective gene transfection of individual cells in vitro with plasmonic nanobubbles, J. Control. Release, 152, pp. 286-293, 2011.


10

4. HARTLAND, G.V., Measurements of the material properties of metal nanoparticles by timeresolved spectroscopy, Phys. Chem. Chem. Phys., 6, pp. 5263-5274, 2004. 5. SCHMID, T., Photoacoustic spectroscopy for process analysis, Anal. Bioanal. Chem., 384, pp. 1071-1086, 2006. 6. VOLKOV, A.N., SEVILLA, C., ZHIGILEV, L.V., Numerical modelling of short pulse laser interaction with Au nanoparticle surrounded by water, Appl. Surf. Sci., 253, pp. 6394-6399, 2007. 7. EGEREV, S., ERMILOV, S., OVCHINNIKOV, O., FOKIN, A., GUZATOV, D., KLIMOV, V., KANAVIN, A., ORAEVSKY, A., Acoustic signals generated by laser-irradiated metal nanoparticles, Appl. Opt., 48, C38-C45, 2009. 8. BRUJAN, E.A., Numerical investigation on the dynamics of cavitation nanobubbles, Microfluid. Nanofluid., 11, pp. 511-517, 2011. 9. BRUJAN, E.A., Collapse of cavitation bubbles in blood, Europhys. Lett., 50, pp. 175-181, 2000. 10. VOGEL, A., NOACK, J., HUTTMAN, G., PALTAUF, G., Mechanisms of femtosecond laser nanosurgery of cells and tissues, Appl. Phys. B, 81, pp. 1015-1047, 2005. 11. BRUJAN, E.A., Dynamics of shock waves and cavitation bubbles in bilinear elastic-plastic media and the implications to short-pulsed laser surgery, Eur. Phys. J. Appl. Phys., 29, pp. 115-123, 2005. 12. ANDERSON, L.J.E., HANSEN, E., LUKIANOVA-HLEB, E.Y., HAFNER, J.H., LAPOTKO, D.O., Optically guided controlled release from liposomes with tunable plasmonic nanobubbles, J. Control. Release, 144, pp. 151-158, 2010. 13. HUTSON, M.S., MA, X., Plasma and cavitation dynamics during pulsed laser microsurgery in vivo. Phys. Rev. Lett., 99, 158104, 2007. 14. SIEMS, A., WEBER, S.A.L., BONEBERG, J., PLECH, A., Thermodynamics of nanosecond nanobubble formation at laser-excited metal nanoparticles, New J. Phys., 13, 043018, 2011. 15. RAYLEIGH, Lord, On the pressure developed in a liquid during the collapse of a spherical cavity, Phil. Mag., 34, pp. 94-98, 1917. 16. WAGNER, D.S., DELK, N.A., LUKIANOVA-HLEB, E.Y., HAFNER, J.H., FARACHCARSON, M.C., LAPOTKO, D.A., The in vivo performance of plasmonic nanobubbles as cell theranostic agents in zebrafish hosting prostate cancer xenografts, Biomater., 31, pp. 7567-7574, 2010. 17. KOTAIDIS, V., DAHMEN, G., VON PLESSEN, G., SPRINGER, F., PLECH, A., Excitation of nanoscale vapor bubbles at the surface of gold nanoparticles in water, J. Chem. Phys., 124, 184702, 2006. 18. DAYTON, P.A., CHOMAS, J.E., LUNN, A.F.H., ALLEN, J.S., LINDNER, J.R., SIMON, S.L., FERRARA, K.W., Optical and acoustical dynamics of microbubble contrast agents inside neutrophils, Biophys. J., 80, pp. 1547-1556, 2000.

NUMERICAL ALGORITHMS FOR SOLVING THE ELASTOPLASTIC PROBLEM WITH MIXED HARDENING SANDA CLEJA-TIGOIU *, NADIA ELENA STOICUTA †, OLIMPIU STOICUTA ‡

Abstract: The paper deals with elasto-plastic models for mixed hardening material, described by rate-type constitutive equations. In the first part of the paper update algorithms for solving rate type elastic-plastic models with mixed hardening mixed for small deformations have been emphasized, based on the radial return mapping algorithms proposed by Simo and Hughes [28]. We exemplify the elastic prediction and plastic corrector steps of the algorithms as well as the calculus of the appropriate algorithmic tangent moduli. We emphasize the role played by the algorithmic tangent moduli in solving the quasi-static, initial and boundary value problem, based on the methodology proposed and developed in the paper by Cleja-Tigoiu and Stoicuta [8]. Key words: small strain elasto-plastic model mixed hardening, rate type constitutive models, Radial Return Algorithm, finite element method, numerical application.

1. INTRODUCTION

The paper deals with the mathematical description of the elasto-plastic models with mixed hardening, within the constitutive framework of small elastoplastic strains. The rate independent type constitutive models were described and discussed by Cleja-Tigoiu and Cristescu [6], Kachanov [15], Khan and Huang [16], Paraschiv-Munteanu et al. [24], Simo and Hughes [28]. The variational formulation of the elasto-plastic problem is proposed and described by Han and Reddy [10], Johnson [14] and Hughes [12]. The general finite element method (FEM) is developed by Fish and Belytschko in [9], Belytschko et al. in [3], and it is applied to solve the elastic type problems, can see also Bath [2]. The assemble matrix of the internal forces and external forces vectors is given by Hughes [13]. In Section 2 we briefly introduce the rate-type constitutive model.

* University of Bucharest, Faculty of Mathematics and Computer Science, 14 Academiei, 010014 Bucharest, Romania † University of Petroşani, Department of Mathematics-Informatics, 20 Institutului, 332006 Petroşani, Romania ‡ University of Petroşani, Department of Control Engineering, Computers, Electrical Engineering and Power Engineering, 20 Institutului, 332006 Petroşani, Romania


Sanda Cleja-Tigoiu, Nadia Elena Stoicuta, Olimpiu Stoicuta

2

In the Section 3 of the paper we presented the update algorithm applied to the rate type independent elasto-plastic model with mixed hardening, which allow us to incrementally describe the elasto-plastic response of the model when the history of the total strain is given. We will focused mainly on the methodology proposed by Simo and Taylor [27], Simo and Hughes [28], which is based on the Radial Return Mapping Algorithm. The elastic prediction and plastic corrector are the steps of the algorithm. The algorithm has been widely studied by Maenchen and Sack [20], Nagtegaal [21] and Simo and Taylor [27]. Other generalized forms of Radial Return Method and their application in solving of some elastic-plastic or viscoplastic problems can be found in the following works: Schreye [26], Loret and Prevost [18] or Ortiz and Popov [23]. The classical Radial Return Method is proposed for the first time by Wilkins [31], Krieg and Key [17] and the development proposed by Simo and Hughes [28] is based on the implicit backward Euler method. The Radial Return Mapping Algorithm proposed by Simo and Hughes [28] provided as well as the calculus of the appropriate algorithmic tangent moduli. In Section 3 we discuss the algorithm proposed by Simo and Hughes [28] for the Prager hardening law, the modified algorithm for a in-plane stress state given by Cleja-Tigoiu and Stoicuta [8], and the algorithmic procedure applied by Lubarda and Benson [19] for Amstrong-Frederick hardening model. We emphasize the role played by the elastic trial stress state and the procedure to calculate the appropriate algorithmic tangent moduli. The update algorithm applied to the rate-type model is coupled with the equilibrium equation to be satisfy by the Cauchy stress at any time. The initial and boundary value problem for the rate-type elasto-plastic material consists of the equilibrium equation to be satisfies by the Cauchy stress, the elasto-plastic constitutive equations, which involves the Cauchy stress, at which the boundary and initial conditions have to be considered. In order to numerically solve the problem the discretized weak form of the equilibrium equation is coupled with the update algorithm. The quasi-static problem with the initial and boundary values associated with the rate-type elasto-plastic model with mixed hardening is formulated in the Section 4.1, following the numerical procedure developed by Cleja-Tigoiu and Stoicuta [8]. The discrete variational formulation of the problem and application of the finite element method is performed in Section 4.2. The nonlinear equation for the discrete displacement vector is solved by Newton-Raphson method, which requests the algorithmic tangent moduli, already calculated. Further we shall use the following notations:  the set of real numbers and  ≤0 = { x ∈  x ≤ 0} ; Lin the second order tensors; Sym ⊂ Lin the space of symmetric tensors; Sym2 the set of symmetric plane 1 tensors; ( e1 , e2 , e3 ) a Cartesian basis ; σ′= devσ= σ − ( trσ ) I 3 the deviatoric part 3

Numerical algorithms for solving the elasto-plastic problem with mixed hardening

3

of the stress tensor; I 3 the identity tensor in  3 ; I 2 = δij ei ⊗ e j the second-order identity tensor in  2 ; σ = ∂σ / ∂t the time derivative of σ ; ∂ C  the partial derivative with respect to C of the function  ;  the fourth order tensor; E the 0, if  < 0 the Young modulus; ν the Poisson ratio; x the norm of x;  ( ) =  1, if  ≥ 0 1 Heaviside function; 〈= z〉 ( z + z ) , ∀z ∈  the positive part of the real number 2 z ; a ⋅ b and A ⋅ B vector and tensor scalar products, respectively; ∇u( x, t ) =ui , j ei ⊗ e j the displacement gradient; H m (Ω) the Sobolev space; H 0m (Ω) closed space C0∞ (Ω) in H m (Ω) . 2. THE CONSTITUTIVE FRAMEWORK OF THE ELASTO-PLASTIC MODEL WITH MIXED HARDENING

We assume that the elasto-plastic body occupies a domain Ω ⊂  3 which is an open and bounded set with smooth boundary ∂Ω and closure Ω . Let us denote by x ( x1 , x 2 , x3 ) ∈ Ω a material point. The time interval is = denoted by = I [0, T ] ⊂  + . We introduce the following notation: • •

u : Ω × I →  3 , u( x, t ) ≡ ui ( x, t )ei the displacement vector; σ : Ω × I → Sym , σ ( x, t ) ≡ σij ei ⊗ e j the Cauchy stress tensor;

•

ε : Ω × I → Sym, ε= ( x, t )

•

tensor; α : Ω × I → Sym , α ( x, t ) ≡ α ij ei ⊗ e j the kinematic hardening variable;

•

k : Ω × I →  the isotropic hardening variable.

1  ∂ui ∂u j +  2  ∂x j ∂xi

  ei ⊗ e j the infinitesimal strain  

The set of the symmetric tensors is defined by Sym. We define the elastic domain: = E

{( σ, α, k ) ∈ Sym × Sym ×   ( σ, α, k ) ≤ 0}

(1)

The set of the symmetric tensors is defined by: Sym :=

{τ ∈ 

3×3

}

τT = τ, τij ∈ L2 (Ω)

(2)


4

where L2 (Ω) is the space of the square integrable functions defined on Ω . The yield function  : Ω ⊂ Sym × Sym ×  →  0 is dependent on the stress and hardening variables and is written in the form  ( σ= , α, k )

σ′ - α −

2 F (k ) 3

(3)

The variable α is called back-stress and characterizes the translational motion of the yield surface in the space of symmetrical tensor, Sym, while F gives the shape of the yield surface. We listed here some representations for the isotropic hardening function F, which is generally an increasing function of k : • the linear isotropic hardening ) Hk + σY F (k=

(4)

where H > 0 is the hardening constant and σY represents the initial yield constant, • the form proposed by Chaboche [5] F (k= ) Q(1 − e −bk ) + σY

(5)

where Q, b > 0, σY are the material constants. • the form introduced by Voce in [29] F (= k ) H k + ( K ∞ − K 0 )(1 − e −bk ) + σY

(6)

where K ∞ > K 0 > 0, b > 0 , H , σY are the material constants. The rate of the plastic strain tensor is described by the associated flow rule through: ε p = λ∂ σ  ( , σ α, k )

(7)

where the plastic factor λ is a function of the material state of material and is defined by the Kuhn-Tucker conditions

and the consistency condition

λ ≥ 0,  ≤ 0, λ =0

(8)

λ = 0

(9)

The elastic type constitutive equation is given by: = σ  (ε − ε p )

(10)

5


with ε e = ε − ε p the elastic strain tensor and  : Sym → Sym is linear mapping and is called the stiffness elastic tensor. The time variation of the isotropic hardening variable is given by the relation: 2 p p ε ⋅ ε 3

= k

(11)

The time variation of the kinematic hardening variable is given by Prager type evolution equation [25], see also Ziegler [30]: α = C ε p

(12)

or by the Armstrong-Frederick type [1] α =

2 p Cε − γk 3

(13)

with C and γ material constant, generally considered to be constant. As a peculiar case we mention the so-called linear hardening law with C = H ′(k ) which has been introduced by Simo and Hughes [28] and which will be considered together with γ =0 in the Section 4.1. Proposition 1. The plastic factor λ : Ω × I →  is calculated on the yield surface from the consistency condition (9) and the Kuhn Tucker conditions (8) by 〈β〉  ( ) h where β the complementary plastic factor is given by the expression λ=

(14)

β=  ∂ σ  ⋅ ε

(15)

and the hardening parameter h is given by h= ∂ σ  ⋅  ∂ σ  + C − γ∂ σ  ⋅ α +

2 F '(k ) 3

(16)

under the hypothesis that h > 0 . Here C , γ, σY are material constants. Proposition 2. The three-dimensional elasto-plastic model with hardening by the given Armstrong-Frederick law is described by the following rate-type constitutive equation


 β σ = ε −  ∂ σ  ( )  h   p β ∂ σ   ( ) ε=  h  β α = ( ) ( C ∂ σ  − γα )   h  k = β  ( )  h

6

(17)

with Q, b, C , γ, σY material constants and associated with the yield function given by (3). The expressions of the complementary plastic factor and the hardening parameter are given by relations (15) and (16). Remark 1. The set of constitutive equations allows us to define a differential system for the unknowns σ, ε p , α and k , when the deformation tensor ε is a given time dependent function. 3. UPDATE ALGORITHM APPLIED TO THE RATE TYPE INDEPENDENT ELASTO-PLASTIC MODEL WITH MIXED HARDENING

In this section we describe some algorithms for solving the differential system which describes the rate-type elasto-plastic constitutive model when the history of the strain tensor ε is given. We discuss the algorithm proposed by Simo and Hughes [28] for the Prager hardening law, the modified algorithm for a inplane stress state given by Cleja-Tigoiu and Stoicuta [8], and the algorithmic procedure applied by Lubarda and Benson [19] for Amstrong-Frederick hardening model. 3.1. The radial return algorithm proposed by Simo and Hughes [28]

The numerical algorithm for solving the rate-type elastic-plastic model with hardening proposed and developed by Simo and Hughes [28] is a predictorcorrector method, which applies the return mapping algorithm. Elastic Predictor step specifies for a given strain increment on the interval [t n , t n +1 ] the trial elastic state, which is defined for the irreversible freezing state. If the trial elastic state remains inside the appropriate yield condition, then an elastic state is realized. If not, namely if the trial state is outside the yield surface, then the state corresponds to plastic behavior. Plastic Corrector step and return mapping algorithm has the purpose to restore the consistency of the plastic state, which means that the corresponding plastic stress has to be on the appropriate yield surface.

7


Finally, the algorithmic elasto-plastic tangent modulus is built. Now we describe the Radial Return Algorithm proposed by Simo and Hughes in [28], for solving the elasto-plastic model with mixed hardening characterized by Prager type law. This elasto-plastic model is described by the following differential type system:

(

= σ  ε − ε p   p s−α ε = λ  F (k )  s−α  α= C λ F (k )  k = λ

) (18)

where s = σ′ is the deviatoric part of the stess tensor. If the time interval [ 0,T ] is discretyzed by the time sequence 0 = t0 < t1 < . . . < t N = T and the notation X (tn +1 ) = X n +1 is used to denote the the value of the field X at time tn +1 , the elastic type constitutive relation at time tn +1 becomes:

σ n +1  (ε n +1 − ε np+1 ) =

(19)

We apply the backward Euler method for the equation from (18), and the following algorithmic relationships are provided p ε np+= 1 ε n + ∆λ

s n +1 − α n +1 s n +1 − α n +1

kn += 1 kn + ∆λ α n += 1 α n + C ∆λ

(20) s n +1 − α n +1 s n +1 − α n +1

where ∆λ = λ∆t = λ (tn +1 − tn ) satisfy the discrete Kuhn-Tucker conditions (21) ∆λ ≥ 0,  ( s n +1 , α n +1 , kn +1 ) ≤ 0, ∆λ ( s n +1 , α n +1 , kn +1 ) = 0 The yield function is defined by:  ( s n +1 ,= α n +1 , kn +1 )

s n +1 − α n +1 −

2 F ( kn +1 ) 3

(22)

Predictor elastic. The trial elastic stress tensor and its deviatoric part, respectively, are defined by freezing the plastic strain at the previous time


(

σ ntr+1 = ka tr ε n +11 + 2µ(dev ε n +1 − ε np ); strn +1 = s n + 2µ dev ε n +1 − ε np

8

)

(23)

As ε np+1 , α n +1 and k n +1 are frozen at their previous values it follows that ε np+1 ≡ ε np , α n +1 ≡ α n , kn +1 ≡ kn

(24)

tr is defined by and consequently, the test function n+ 1

= ntr+1 :

s ntr+1 − α n −

2 F ( kn ) 3

(25)

Corrector plastic. Return algorithm. If the elastic trial state is outside the yield surface, then ∆λ cannot be zero and the following relations take place tr s n += 1 s n +1 − 2µ∆λ

s n +1 − α n +1 s n +1 − α n +1 ; α n += 1 α n + C ∆λ s n +1 − α n +1 s n +1 − α n +1

(26)

The unit normal n trn+1 is determined exclusively in terms of the trial elastic stress: ntrn +1 =

strn +1 − α n

(27)

strn+1 − α n

tr tr If we denote with ξ= n +1 s n +1 − α n +1 and ξ = n +1 s n +1 − α n . Then we have

tr tr ξ= n +1 s n +1 − α n − 2µ∆λ ( 2µ + C ) n n +1

(28)

Proposition 3. (The update algorithm). The following loading/unloading conditions hold: tr If n+ 1 ≤ 0 then ∆λ = 0 and the solution is elastic

(

σ ntr+1 = ka tr ε n +11 + 2µ(dev ε n +1 − ε np ); s ntr+1 = s n + 2µ dev ε n +1 − ε np ε np+1

≡ ε np ,

α n +1 ≡ α n , kn +1 ≡ kn

)

(29)

Precisely these relationships are taken to be the initial conditions for the solution of the plastic corrector problem. The consistency of the plastic condition is restored by the return-mapping algorithm. tr If n+ 1 > 0 then ∆λ ≠ 0 and the algorithmic expression of the plastic factor

λ n +1 is calculated from the consistency condition characterized by Kuhn-Tucker conditions. Hence, the plastic factor results

9


∆λ =

2 F ( kn ) 3 , 2 F ′(kn ) + 2µ + C 3

ξ ntr+1 −

(30)

and the update values of the state at time tn +1 tr s n += 1 s n +1 − 2µ∆λ

ε np+= 1

ε np

s n +1 − α n +1 s n +1 − α n +1 , α n += , 1 α n + C ∆λ s n +1 − α n +1 s n +1 − α n +1

s − α n +1 + ∆λ n +1 , s n +1 − α n +1

(31)

kn += 1 kn + ∆λ.

The radial return mapping algorithm is completely determined, if the algorithmic elastic-plastic tangent modulus is calculated. Proposition 4. (Simo and Hughes [28]) The algorithmic elasto-plastic tangent modulus can be approximated by the following expression:  ,  tr  0 nep+1 =  tr ka 1 ⊗ 1 + 2µ a [ θ1 ]n +1 I dev − 2µ a [ θ2 ]n +1 n n +1 ⊗ n n +1 ,  > 0

[θ1 ]n+1 =

1 − ∆λ

2µ a ξ ntr+1

;

[θ2 ]n+1 =

2µ a 2 F ′(kn ) + 2µ a + C 3

− ∆λ

2µ a ξ ntr+1

(32)

(33)

3.2. The algorithm proposed by Cleja-Tigoiu and Stoicuta [8]

In [8] a new numerical approach of the elasto-plastic problem with mixed hardening in the plane stress state, within the classical constitutive framework of small deformation, has been proposed. The numerical algorithm (which is an incremental one) consists of the Radial Return Mapping Algorithm, originally proposed by Simo and Hughes [28], adapted to the stress plane problem and coupled with a pseudo-elastic problem. The proposed algorithm has been called the revisited Simo-Hughes algorithm and it is applied to exemplify the mode of integration of the bidimensional problem. • We restrict ourselves to a plane stress state σ=

2

∑ σij ei ⊗ e j ,

i , j =1

σ ∈ Sym2 , i. e. σ ' ≡ σ (2)

(34)


10

since σ13 = σ23 = σ33 = 0 . Here is avoid the representation such as (σ11 , σ22 , 2σ12 ) is avoided, and for a plane symmetric tensor σ we have the representation σ = σ11e1 ⊗ e1 + σ22 e2 ⊗ e2 + σ12 e1 ⊗ e2 + σ21e2 ⊗ e1 , with σ12 = σ21

•

(35)

The kinematic hardening variable α is defined: α=

3

∑ αij ei ⊗ e j ,

α ∈ Sym

(36)

i , j =1

and α13 = α 23 = 0 , by trace null tr α =0 . • The generalized strain state is characterized by a strain tensor ε that has a non-zero component in the direction perpendicular to the shell, denoted by ε33 , which is added to the in plane strain ε(2) : ε = ε(2) + ε33e3 ⊗ e3 , ε(2) =

2

∑ εij ei ⊗ e j ,

ε(2) ∈ Sym2 .

(37)

i , j =1

•

The tensorial representations for the plain stress state and for the plane kinematic variable: σ ' = σ′(2) + σ′33e3 ⊗ e3 ; α = α (2) + α33e3 ⊗ e3 ; s = σ′(2) ; αˆ = α (2) ; σ′33 = − tr s; tr s = s11 + s22 ; α33 = − tr αˆ ; tr αˆ = αˆ 11 + αˆ 22

•

(38)

The passage from σ ≡ σ (2) ⊂ Sym2 to s = σ '(2) can be characterized by an invertible mapping P given by the matrix  2 / 3 −1/ 3  −1/ 3 2 / 3 = P   0 0  0  0

0 0 1 0

The plastic factor λ : Ω × I →   ( s, αˆ , k ) = 0 by: λ=

where

0  0 = , s P= σ , σ σ (2) . 0  1

(39)

is calculated on the yield surface

〈β〉  ( ) h

(40)


11

β=

E  ν E  ν 2E   ∂ trs   ε11 + ∂ trs   ε 22 + ∂ s ε12  ∂ s11  +  ∂ s22  + 1+ ν  1− ν 1+ ν  1− ν 1 + ν 12  

E (1 − 2ν)  2 2  E 2  E  = h  + C  ∂s +  +C + H  ∂ trs  + 2 3 1− ν  1+ ν  1+ ν

(41)

Here C is the kinematic hardening modulus. The variation in time of the component of the strain ε33 is defined by ε 33 = −

ν 1 − 2ν λ ∂ trs  ( ε11 + ε 22 ) − 1− ν 1− ν

(42)

Remark 2. The evolution equation for the strain component ε33 is derived from the consistency requirement to have a plane stress state, applying the procedure developed by Cleja-Tigoiu [7]. In the case of a plane stress state associated with the generalized plane strain, the elasto-plastic tangent modulus can be expressed in the following modified form 

ep

(2)  =  E(2) [Θ1 ] − [Θ2 ] n ⊗ E(2) P n + n ⊥ E(2) P I (2)

(

if  tr ≤ 0

)

if  tr > 0

(43)

where [Θ1 ] = I 4 −

[λ* ]n +1  s − αˆ tr  2 + (tr( s tr − αˆ tr ))2 tr

  1 = − [Θ2 ] (2)   2 E +C H+  3(1 − ν)  3

PE(2)

  [λ ]n +1  tr tr 2 tr tr 2   s − αˆ  + (tr( s − αˆ ))   *

(44)

In the relation above, the tensor E(2) : Sym2 × Sym2 is the elastic tensor and the normals n and n ⊥ have the following form n≡ n⊥ ≡

s − αˆ

s tr − αˆ tr =  s − αˆ 2 + (tr(s − αˆ )) 2  s tr − αˆ tr  2 + (tr( s tr − αˆ tr ))2 tr(s − αˆ )

tr(s tr − αˆ tr ) =  s − αˆ 2 + (tr(s − αˆ )) 2  s tr − αˆ tr  2 + (tr( s tr − αˆ tr ))2

(45)


12

Only when the plane stress state is considered to be associated with the plane strain, namely when the strain in the direction perpendicular to the plane is neglected, the elasto-plastic tangent moduli provided here can be reduced to the formulae that can be found in Simo and Hughes [28], namely the appropriate fourth order tensors are replaced by two scalars, written here in the formulae (32). 3.3. Algorithmic consistency condition proposed by Ludarda and Benson [19]

We exemplify the update algorithm proposed in [19] for the rate-type model, which is written in Proposition 2. The elasto-plastic model with mixed hardening characterized by the Armstrong-Frederick type rule is described by the following differential system: s = 2µe − 2µε p  ε p = λ∂ s   α C ε p − γ α k =   p k = ε

(46)

where β  λ = h c  β = ∂ s  ⋅ s  (s − α ) ⋅ α hc 2 F ′(k ) + C − γ =  2 F (k )

(47)

and which is associated with the yield surface given  ( s= , α, k ) 1/ 2 s − α − F (k )

(48)

To determine the algorithmic consistency condition, the authors applied the Euler forward method s n +1 =  kn += 1  = α  n +1  p n +1 ε=

where

s n + 2µde − 2µλ*n +1n n +1 kn + λ*n +1 α n + (dα ) n +1 ε np + λ*n +1n n +1

(49)


13

s n +1 − α n +1  n n +1 = s − α n +1 n +1   * (dα ) n +1 = λ n +1 (C n n +1 − γα n +1 )  * ∆t λ (tn +1 ) λ n +1 =  

(50)

In the method proposed in [19] in the second member of the equality written in (50)2 α n+1 is replaced with θ α n + (1 − θ)α n+1 . Thus the values at time tn +1 has been evaluated from the equality α n +1 − α n = λ*n +1 (Cnn +1 − γ (θα n + (1 − θ)α n +1 ))

(51)

where θ∈ [ 0,1] . In the case θ =0 the forward Euler rule takes place, θ =1 1 is used in the generalized 2 midpoint rule. We remark that a hybrid method has been employed to approximate the rate type equations which describe the model. Therefore we derived that dα= n +1 α n +1 − α n with

corresponds to the backward Euler rule and θ =

γ   dα n +1 =λ an +1 *n +1  n n +1 − α n  C  

(52)

where an +1 =

C

(53)

1 + γ (1 − θ)λ*n +1

The following equality results γ  * *  s n +1 − α n += αn  1 s n − α n + 2µ de − 2µλ n +1n n +1 − an +1λ n +1  n n +1 − C  

(54)

This can be finally written in the form s n +1 − α n +1 + (2µ + an +1 )λ*n +1n n +1= s n − α n + 2µde + bn +1λ*n +1α n

(55)

with bn +1 =

γ 1 + γ (1 − θ)λ*n +1

.

(56)


14

Remark 3. Let us denote by s *n+1 the trial elastic stress s *n +1 = s n + 2µde that corresponds to the elastic strain increment de. First the frozen plastic behavior has to be considered. Only if the trial elastic state is outside the appropriate trial yield condition at time tn +1 the approximate solution has to be characterized by the non-zero plastic factor. Following the method proposed by Simo and Hughes [28] the procedure to restore the consistency condition has to be pursued. If we take a trace product for the equation given above taking into account the fact that the plasticity function F ( , σ′ − α =2 F (k ) then σ α, k ) = 0 , namely it becomes: 2 F (kn + 2λ*n +1 ) + (2µ + an +1 )λ*n +1 =  =  2 F 2 (kn ) + 2µde + bn +1λ*n +1α n 

2

 + 2(s n − α n ) ⋅ (2µde + bn +1λ*n +1α n )  

1/ 2

(57)

That is the key equation of the numerical method proposed by Lubarda and Benson [19]. There is the algorithmic consistency condition for the considered hardening model with the Armstrong-Frederick evolution of the back stress. The nonlinear equation for the loading parameter λ *n+1 is proposed to be resolved using Newton’s iterative method. = G

2 F (kn + 2λ*n +1 ) + (2µ + an +1 )λ*n +1 −

 −  2 F 2 (kn ) + 2µde + bn +1λ*n +1α n 

2

 + 2(s n − α n ) ⋅ (2µde + bn +1λ*n +1α n ) 

1/ 2

(58)



The iterative procedure to find the zero value of G is then based on the approximate formula λ*,n +j1+1 = λ*,n +j1 −

G(λ*,n +j1 ) G ′(λ*,n +j1 )

(59)

where with G ′(λ*,n+j1 ) was denoted the derivative of the function G(λ*,n+j1 ) . The iterations are ended when a desired accuracy has been reached, i.e., when G(λ*,n+j1 ) falls into a prescribed error tolerance. Remark 4. The numerical procedure proposed in [19] has not been applied to solve any elasto-plastic non-homogeneous problem. The homogeneous simple shear deformation process has been considered and compared to the other models.


15

4. NUMERICAL ALGORITHM OF SOLVING THE ELASTO-PLASTIC PROBLEM WITH MIXED HARDENING

The quasi-static elasto-plastic problem with the initial and boundary values is formulated within constitutive framework of mixed hardening material with small strains. The discrete variational formulation of the problem is provided and the finite element method is performed in Section 4.2 in order to solve the problem, following the procedure applied by Cleja-Tigoiu and Stoicuta [8]. The rate type constitutive model considered here is the one introduced in Simo and Hughes [28]. The model has been called non-linear hardening, due to the presence of the scalar hardening variable k in the evolution equation of the Prager-type for the back stress., i.e. the constant C is replaced by H ′(k) . The Radial Return Algorithm is written to determine the algorithmic expression of the plastic factor λ from the consistency condition when the history t → ε (u (t )) is given in a fixed material particle and the tangent elasto-plastic moduli have also been calculated. The initial and boundary value problem will be formulated as in Cleja-Tigoiu and Stoicuta [8]. 4.1. Mathematical formulation of the elasto-plastic problem with mixed hardening

Problem P. Let be given the function f , g : Ω →  3 . Find the functions σ, ε p , α, k which are defined in Ω × [0, T ) , which satisfy the differential type constitutive equations:

= σ  ε(u ) −  λ ∂ σ   p ε = λ ∂ σ   2 α = λ H ′(k ) ∂ σ  3   2 = k λ 3 

where

(60)


 2 F (k )  = σ′ − α − 3   σ′ − α ∂ σ  =′ σ − α   β  λ = h H ( F ) c  β = ∂ σ  ⋅  ε(u )  2 2  F ′(k ) ∂ σ  ⋅  ∂ σ  + ∂ σ  ⋅ H ′( k )∂ σ  − hc = 3 3

16

(61)

The Cauchy stress σ satisfies quasi-static boundary value problem σ+b = div 0  σn f =   = u g

in Ω × I in Γ σ × I in Γu × I

(62)

at which the initial conditions have to be added , σ (0) σ= = u(0) u 0= 0 , ε (0) ε 0  p p = ε (0) ε= 0 , k (0) k0 0 , α (0) α=

(63)

To solve numerically the Problem P, we propose the following strategy: • we introduce the weak formulation associated with the equations in order to find the stress σ n+1 and the displacement vector u n+1 at the time t n +1 ; • •

for a given incremental strain history, an update algorithm for s, ε p , α, k is derived by using the constitutive equation; at every time t n +1 and for any material point in x ∈ Ω(tn +1 ), we use an iterative solution of the rate type given by the update algorithm, to constitutive equation boundary and initial value problem for the discretized equilibrium balance equation and coupled with the update algorithm.

4.2. Discretized variational formulation and application of the finite element method

We proceed to discretize the weak formulation of the problem to solve an associate non-linear elastic-like problem. We supposed that the boundary of the domain denoted with ∂Ω is decomposed into two parts Γ u and Γ σ , with Γ u ∪ Γ σ = ∂Ω and Γ u ∩ Γ σ = ∅ .

17


We define the set of kinematically admissible velocity field at every time t , is denoted by: ad =

 ∂wi 2 ∈ L2 (Ω), wi  wi ∈ L (Ω); ∂x j 

Γu

 = gi , i, j = 1, 2  ⊂ H 1 (Ω) 

(

)

2

(64)

where H 1 (Ω) is the Sobolev space. We replace the continuous time interval [ 0,T ] by the sequence of discrete times t0 , t1 ,..., t N with tn +1= tn + ∆t and we formulate at the moment tn +1 the pseudo-elastic problem Pne+1 . Problem Pne+1. Find the displacement field u n+1 , solution of the discretized weak formulation: L(w n +1 ) ∫  (εn+1 − εn+1 ) ⋅ ε(w n+1 ) dx = p

(65)

Ω

with L ( w n +1 ) given by L(w n +1= )

σ n +1n ⋅ g dΓ + ∫ b ⋅ w n +1dx, ∫ f ⋅ w n+1dΓ + ∫

Γσ

Γu

(66)

Ω

for all w n +1 ∈ adh , which is the finite-dimensional approximation to ad . At each iteration n + 1 one solves a so called pseudo-elastic problem Pne+1 , where w n+1 is the test function associated with the n + 1 iteration. The finite element method is utilized to solve the discrete variational formulation of the problem. Proposition 5. The discrete version of the weak formulation, namely of the problem Pne+1 reads: int Gˆ ( uˆ= Fnext ( uˆ n+1 ) − = +1 0. n +1 ) F

(67)

Proof. We proceed to discretize the weak form to give the associate nonlinear pseudo-elastic-like problem. Thus, the displacement field uˆ n+1 has be solution of the following equations


18

1 1    T T  e  e e e   = ˆ ˆ ( , ) , ( , ) d d w σ u B J ξ η ξ η ξ η ξ η ( )  n +1 n +1     e =1   −1 −1   1 1 nel    e T ( N (ξ, η) )T N (ξ, η)bˆ e J e (ξ, η)  dξdη + ˆ n +1 =  w n +1    e =1  −1 −1 

∑(

) ∫∫ (

nel

∑(

) (

)

) ∫∫

 ˆ en +1 +  w =i 1 =e 1   2 nel   e ˆ n +1 +  w  =i 1 =e 1 

 ( N (ξ, −2 ( i − 1) + 1) ) N (ξ, −2 ( i − 1) + 1)fˆ e ( J )e  dξ  +  1 1 + n i    −1  1  T ( N (−2 ( i − 1) + 1, η) )T N (−2 ( i − 1) + 1, η)fˆ e ( J )e  dη  1 2 n + i     −1  1

∑∑ (

) ∫

∑∑ (

) ∫

2 nel

T

(68)

T

T

e = where ( J1 ) = ( J )e ( J1 )e2  ; ( J 2 )e ( J 2 )1e  11

( J 2 )e2 

T

with

( J1 )1e

(( x ) − ( x ) ) + (( x ) − ( x ) ) ; ( J ) (( x ) − ( x ) ) + (( x ) − ( x ) ) =

( J 2 )1e

( x ) − ( x ) ) + (( x ) − ( x ) ) ( x ) − ( x ) ) + (( x ) − ( x ) ) ( (= ; (J )

e 1 3

e 2 1 4

e 2 3

e 2 2 4

e 1 2

2

e 1 2

e 2 1 3

e 2 2

e 1 1

e 2 2 3

e 2 1

e 2 2 2

2

e 1 1

e 2 2

2

e 2 1 2

e 2 1 4

e 2 1

e 2 2 4

2

(69) ˆ ∈ adh w

. for any In equality (68) the shape function matrix N (ξ, η) is by the form: N N (ξ, η) =  1 0

0 N1

N2 0

0 N2

N3 0

0 N3

N4 0

0  N 4 

(70)

with 1 (1 − ξ )(1 − η) ; N 2 ( ξ, η=) 4 1 N3 ( ξ, η = ) (1 + ξ )(1 + η) ; N 4 ( ξ, η=) 4 N1 ( ξ, η = )

and the strain-displacement matrix B e (ξ, η)

1 (1 + ξ )(1 − η) 4 1 (1 − ξ )(1 + η) 4

(71)

;


19

e  B11   0 B e (ξ, η) = e  B21  e  B12

0

e B12

0

e B13

0

e B14

e B21

0

e B22

0

e B23

0

e B11 e B11

e B22 e B22

e B12 e B21

e B23 e B32

e B13 e B31

e B24 e B42

0   e B24  e  B14  e  B41 

(72)

In the discrete form (68), for simplicity, we introduce the following notations:

(

σ ( uˆ ) ( ξ, η) J ) ∫ ∫ ( B (ξ, η) ) 1

 f int uˆ en +1  =   e

= f ext   1  n +1

(

)

e

(

)

e

1

T

e

−1 −1

1

1

∫−1 ∫−1 ( N (ξ, η) )

T

e n +1

e

 (ξ, η)  dξdη 

N (ξ, η)bˆ en +1 J e (ξ, η)  dξdη 

( N (ξ, −2 ( i − 1) + 1) )T N (ξ, −2 ( i − 1) + 1)fˆ e ( J )e  dξ n +1 1 i   −1  

 f ext =   2 i  n +1  f ext=   3 i  n +1

1

∫

∫−1 ( N (−2 ( i − 1) + 1, η) ) 1

T

(73)

e N (−2 ( i − 1) + 1, η)fˆne+1 ( J 2 )i  dη 

The integrals of the form (73) are approximated by the Gauss method. Given the above, the equality (68) can be written as

{

∑( nel

)

) } ∑{( wˆ )

(

T

nel

ˆ en +1  f int uˆ en +1  = w  

T e n +1

e

=e 1 =e 1 2 nel  e T  ext  e  2 nel  ˆ n +1 f + +  w   1 i  n +1     =i 1 =e 1 =i 1 =e 1 

∑∑ (

) (

)

}

 f3ext  +   n +1

∑∑ (

e  f ext    2 i   n +1 

) (

T ˆ en +1 w

)

(74)

In the following, the assembled matrix of internal forces and vectors of the external forces are denoted by = F int ( uˆ n +1 )

(

nel

)

(

)

nel

e

= A  f int uˆ en +1  ; F1ext A  f1ext  ;  n +1  n +1 e 1  =e 1  =  A   f 2ext  n +1 e =1   nel  = F3ext A   f3ext n +1 e =1   

(

)

(

)

= F2ext

nel

(

) 

e

(

)

e

1

) 

(

)

 +  f ext 1  n +1  3

We also use the following notation:

 ; 2 n +1  e   2  n +1 

(

+  f 2ext  n +1 

e

(75)


Fnext +1 =

∑ ( Fiext )n+1.

20

3

(76)

i =1

Then, discrete variational formulation (74) yields  ( uˆ , w  ˆ Tn +1  F int ( uˆ n +1 ) − Fnext ˆ ∈ adh G n +1 ˆ n +1 )= w +1 = 0, ∀w 

(77)

Consequently, the resulting system of non-linear equations (77) can be written as (67). 4.3. Radial Return Algorithm for the model with non-linear hardening

We use the Radial Return Mapping Algorithm which has been presented in the previous section 3.1, this time for the evolution equations written in (60). At the moment t n +1 we have the elastic type constitutive equation (19) together with s − α n +1 ε np+1= ε np + λ*n +1 n +1 ; kn +1= kn + λ*n +1; s n +1 − α n +1 (78) s n +1 − α n +1 2 α n += ( H (kn +1 ) - H (kn )) 1 αn + s n +1 − α n +1 3 where λ*n +1 = λ∆t = λ (tn +1 − tn ) satisfy the discrete Kuhn-Tucker conditions  (s n +1 , α n +1 , kn +1 ) ≤ 0, λ*n +1 ≥ 0,

0 λ*n +1 (s n +1 , α n +1 , kn +1 ) =

(79)

The yield function is defined by: , α n +1 , kn +1 )  ( s n +1=

s n +1 − α n +1 −

2 F (kn +1 ) 3

(80)

Elastic Predictor. We define trial elastic stress tensor like in (23) and in this stage ε np+1 , α n +1 and kn +1 are frozen as previous values like in (24). The test tr function n+ 1 is defined as in (25). The second step of the algorithm, Plastic corrector and return mapping algorithm is further applied. If the elastic trial state is outside the yield surface, then the plastic factor can not be zero and the procedure to restore the consistency condition follows.

21


tr * s= n +1 s n +1 − 2µλ n +1

s n +1 − α n +1 s n +1 − α n +1

(81)

s −α 2 α n += ( H (kn+1 ) - H (kn ) ) n+1 n+1 1 αn + s n +1 − α n +1 3

tr tr We introduce the notation β= n +1 s n +1 − α n +1 and β= n +1 s n +1 − α n and as a tr , i.e. direct consequence of the collinearity of β n+1 and β n+ 1

= n n +1

βn +1 = βn +1

βntr+1

(82)

βntr+1

In this conditions, from (81) we have:   2 tr * β= n +1 βn +1 −   2µλ n +1 + 3 ( H (kn +1 ) - H (kn ) )  nn +1  

(83)

Proposition 6. The following loading/unloading conditions take place tr * If n+ 0 and the solution is an elastic one 1 ≤ 0 then λ n+1 =

(

σ trn +1 = ka tr ε n +11 + 2µ(ε n +1 − ε np ); strn +1 = s n + 2µ dev ε n +1 − ε np

)

ε np+1 ≡ ε np , α n +1 ≡ α n , kn +1 ≡ kn

(84)

tr * If n+ 1 > 0 then λ n+1 ≠ 0 so the plastic factor λ n+1 is the solution of the nonlinear equation:

βntr+1 −

  2 2 2 * F (kn +1 ) −  2µλ*n +1 + λ n +1 ) - H (kn )   = 0  H (kn + 3 3 3   

(85)

Moreover, the plastic solution is described by the following relationships tr * s= n +1 s n +1 − 2µλ n +1

α n += 1 αn +

s n +1 − α n +1 , s n +1 − α n +1

s −α 2 ( H (kn+1 ) - H (kn ) ) n+1 n+1 , 3 s n +1 − α n +1

ε np+1= ε np + λ*n +1

s n +1 − α n +1 , s n +1 − α n +1

kn +1= kn +

2 * λ n +1 3

(86)


22

Proof. The plastic factor λ*n+1 is calculated from the discrete consistency condition. From the relation (82) we have: = βn +1

= βn +1 n n +1 and βntr+1

βntr+1 n n +1

(87)

and the equality (83) becomes: = βn +1 n n +1

  2 βntr+1 n n +1 −  2µλ*n +1 + ( H (kn+1 ) - H (kn ) )  n n+1 3  

(88)

or 2µλ*n +1 +

2 ( kn ) ) ( H (kn+1 ) − H = 3

(

βntr+1 − βn +1

(89)

)

On the yield surface  β n +1 , k n +1 = 0 the equality β n +1 = F (k n +1 ) holds, where β= n +1

s n +1 − α n +1 and − b( k n +  F (kn += 1 ) Q 1 − e 

2 / 3λ*n+1

)+σ  

Y

.

(90)

Consequently, the plastic factor λ *n+1 is determined from the following nonlinear equation βntr+1 −

  2 2 2 * F (kn +1 ) −  2µλ*n +1 + λ n +1 ) - H (kn )   = 0  H (kn + 3 3 3   

(91)

To solve the equation (91) we apply the Newton method, which consists in the following iterative procedure: λ*,n +j1+1 = λ*,n +j1 −

f (λ*,n +j1 ) f ′(λ*,n +j1 )

(92)

where the function f (λ*,n +j1 ) is given by j f (λ*n,+= 1)

βntr+1 −

  2 2 2 *, j F (knj+1 ) −  2µλ*,n +j1 + λ n +1 ) - H (kn )   (93)  H (kn + 3 3 3   

and f ′(λ*,n +j1 ) is the derivative of the function f (λ*,n +j1 ) :


23

f ′(λ*n +1 ) =−2µ −

 H ′(kn +1 ) + F ′(kn +1 )  2 ( H ′(kn+1 ) + F ′(kn+1 ) ) =−2µ 1 +  (94) 3 3µ  

where F ′(k n +1 ) = Qb e −b ( k n+1 ) in numerical application. The iterations are ended when a desired accuracy has been reached, i.e., when f (λ*,n +j1 ) falls into a prescribed error tolerance.

The Radial Return Algorithm is completely determined, if the algorithmic elastic-plastic tangent modulus  ep is calculated. Proposition 7. The algorithmic elasto-plastic tangent modulus can be approximated by the following expression:  ,  tr  0 e  ep  =    n +1 tr ka 1 ⊗ 1 + 2µ a [ θ1 ]n +1 Idev − 2µ a [ θ2 ]n +1 n n +1 ⊗ n n +1 ,  > 0

(95)

where

[θ1 ]n+1 =

1 − 2µλ*n +1 βntr+1

;

 H ′(kn +1 ) + F ′(kn +1 )  [θ2 ]n+1 = 1 +  3µ  

−1

(96)

(

− 1 − [ θ1 ]n +1

)

Proof. The algorithmic elasto-plastic tangent modulus is calculated from the relation:  ep =

∂σ n +1 ∂ε(uˆn +1 )

(97)

We differentiate the discrete constitutive equation (19) with respect to ε : dσ n +1 d * p  ka ( tr ε ( uˆn +1 ) )1 + 2µ(devε ( uˆn=  = +1 ) − ε n ) − 2µλ n +1n n +1  dε ( uˆn +1 ) dε ( uˆn +1 )  = ka 1 ⊗ 1 + 2µIdev

 ∂nn +1 ∂λ*n +1  − 2µ  λ*n +1 + n n +1 ⊗   ∂ε ( uˆn +1 ) ∂ε ( uˆn +1 )  

(98)

In the relation (98) we observe that the term ka 1 ⊗ 1 + 2µIdev = is the tensor of elastic moduli, where Idev =− I 1 1 ⊗ 1 is the second order deviator tensor. 3 We calculate in the following, the expression of the derivatives witch appearing in (98). Thus


∂n n +1 ∂n n +1 ∂βntr+1 = = ∂ε ( uˆn +1 ) ∂βntr+1 ∂ε ( uˆn +1 )

( I − n n+1 ⊗ n n+1 )

2µ βntr+1

24

Idev

(99)

The derivative of the plastic factor with respect to the total strain at time tn +1 can be derived from the consistency condition written in (9)    ∂  tr 2 2 2 *  βn +1 − F (k n +1 ) −  2µλ*n +1 + λ n +1 ) - H (kn )    = 0  H (kn + ∂ε n +1  3 3 3     ∂λ*n +1   H ′(kn +1 ) + F ′(kn +1 )   2µn n +1Idev  2µ 1 +  = ∂ε ( uˆn +1 )   3µ 

(100)

−1

∂λ*n +1  H ′(kn +1 ) + F ′(kn +1 )  = 1 +  n n +1 ∂ε ( uˆn +1 )  3µ 

Consequently, the following formula can be proved dσ n +1 = ka 1 ⊗ 1 + 2µ a [ θ1 ]n +1 Idev − 2µ a [ θ2 ]n +1 n n +1 ⊗ n n +1 dε ( uˆn +1 )

(101)

where

[θ1 ]n+1 =

1 − 2µλ*n +1 βntr+1

;

 H ′(kn +1 ) + F ′(kn +1 )  [θ2 ]n+1 = 1 +  3µ  

−1

(102)

(

− 1 − [ θ1 ]n +1

)

Newton Raphson algorithm The algorithmic elastic-plastic tangent modulus is involved in the Newton Raphson method which is applied for solving the nonlinear equation (67). The unknown uˆ n +1 is obtained based on the following iterative relation: −1

(

)

(103) uˆ nj ++11 = uˆ nj +1 − β j  K nj+1  Gˆ uˆ nj +1   where β j ∈ ( 0,1] and the Jacobian of the function Gˆ is given by the following expression:

25


∂Gˆ i ( uˆ n +1 )  K (= uˆ n +1 )  = i ,k ∂ ( uˆ n +1 )k =

 e ep ∫−1 ∫−1 ( B (ξ, η) ) (  )

nel  1

A e =1 

T

1

( )

  int e   ∂ f  i    n +1  A = e =1  ∂ ( u ˆ n +1 )k      nel 

(104)

  B e (ξ, η) J e (ξ, η)  dξdη  

e

Finally the Radial Return Algorithm can be assembled

{[ ] σ

e, j +1 e, j +1 ε p  , , n +1   n +1

[α ]en,+j1+1 , [ k ]en,+j1+1

{

}

= e

= Radial Return Algorithm u en,+j1 , [σ ]n +1 , ε p  , [α ]n +1 , [ k ]n  n e

e

e

}

j +1

1. compute the trial elastic stress for prescribed [ε ]n +1 = [ B ] uˆ nj +1 : σ trn +1 = ka tr ε n +11 + 2µ(ε n +1 − ε np );  tr : 2. check the yield condition  = n +1

(

strn +1 = s n + 2µ dev ε n +1 − ε np βntr+1 −

)

2 Q(1 − e −bkn ) + σY  :  3

tr If n+ 1 ≤ 0 , then: • elastic step:

σ n +1 = σ trn +1; s n +1 = strn +1; ε np+1 ≡ ε np ; α n +1 ≡ α n ; kn +1 ≡ kn ; nep+1 =  else • plastic step: proceed to step 3. endif

3. Radial Return Algorithm: • compute n n+1 and find λ*n+1 with Newton method: n n +1 =

βntr+1 βntr+1

Initialize: = λ*n,0+1 0= and kn0+1 0 Do until f (λ*,n +j1 ) < ERR , j ← j +1 Compute iterate


j f (λ*n,+= 1)

βntr+1 −

f ′(λ*n +1 ) =−2µ −

26

  2 2 2 *, j F (knj+1 ) −  2µλ*,n +j1 + λ n +1 ) - H (kn )    H (kn + 3 3 3     H ′(kn +1 ) + F ′(kn +1 )  2 ( H ′(kn+1 ) + F ′(kn+1 ) ) =−2µ 1 +  3 3µ  

λ*,n +j1+1 = λ*,n +j1 −

f (λ*,n +j1 ) f ′(λ*,n +j1 )

• update isotropic hardening variable: 2 *, j +1 . knj++11 = knj + λ 3 n +1 • update stress, hardening variable and compute the algorithmic elastoplastic tangent moduli: * tr s n += 1 s n +1 − 2µλ n +1n n +1 ; α n += 1 αn +

2 ( H (kn+1 ) - H (kn ) ) n n+1 ; 3

* p ε np+= 1 ε n + λ n +1n n +1 ;

nep+1 = ka 1 ⊗ 1 + 2µ a [ θ1 ]n +1 Idev − 2µ a [ θ2 ]n +1 n n +1 ⊗ n n +1 ; = [θ1 ]n+1

1 − 2µλ*n +1 ξtrn+1

;



H ′(k

) + F ′(k 3µ

)

−1

n +1 n +1 [θ2 ]n+1 = 1 +  − (1 − [ θ1 ]n +1 )





Return Remark 5. The discrete functions [ θ1 ]n+1 and [ θ2 ]n+1 are dependent on the evolution equation for the back stress, which is involved into the model, as can be seen from (96) and (33). 5. NUMERICAL APPLICATION

In order to exemplify the mode of integration of the problem, we consider the following bi-dimensional example. Numerical application presented here, aims to highlight the numerical algorithm of solving the Problem P. For this, we chose the trapezoidal panel domain (Figure 1), for which the vertical left edge is fixed, and the bottom and the right vertical edges are traction free, i.e. f = 0 . The traction f = 50sin t is applied on the top horizontal edge, while the internal forces are assumed to vanish b = 0 . Inside of this plates is inserted a hole with a radius of

27


20 mm. The application is realized for a time interval t ∈ [ 0,60] . The sheet is made of steel DP600, with the material parameters given by Broggiato et al. [4] = = C 17400 [ MPa= ; σY 394, 4 [ MPa= ]; H 1194[MPa]} ]; ν 0,3; {E 182000[ MPa ]= The Ω domain for this example is by the form:  x1   2 = Ω 2 ( x1 , x2 ) ⊂  0  x1  20 , 4  x2  100    Ω = Ω 2 \ Ω1 →  2 2 Ω = x , x ⊂  2 ( x1 − 80 ) + ( x2 − 60 )  20  1 ( 1 2 )

{

}

(105)

The numerical example put into evidence the behavior of the sheet, made of an elasto-plastic material with mixed hardening, when the Prager hardening law describes the evolution in time of the back stress.

Figure 1 – The geometry of the plate with the boundary

σij conditions.

In Cleja-Tigoiu and Stoicuta’s paper [7] the numerical procedure to solve the boundary value problem has been illustrated for the deformation in plane stress of the trapezoidal sheet. We exemplify here the behavior in the sheet through the graphs of the stress components , versus the corresponding strain components ij , produced in the four points of the trapezoidal plate with a hole. The cyclic behavior can be observed, as well as the peculiar behavior in compression and tension in the chosen material points. We specify also that for the meshing, the trapezoidal plate with hole was divided into 2500 elements, each element having four nodes. The distribution of the stress component in the direction of the applied force σ11 , and the intensity of the Cauchy stress, i.e. the so-called Mises stress, which is denoted by σMises , are plotted on the trapezoidal plate in Figures 2,3. The graphs are obtained for mixed hardening elasto-plastic materials, using the Cleja-Stoicuta solution for the in-plane stress and linear isotropic hardening


28

given by (4) and the Simo-Hudges algorithm, with the nonlinear isotropic hardening described by (5). The Figures 4-9 show the components of the stress versus the appropriate components of the strain. The components of the total strain p are plotted as function of time at the mentioned ε33 and the plastic strain ε33 material point in Figure 10.

Figure 2 – The distribution of the stress σ 11 on the trapezoidal plate with hole at t = 1,5 and

t = 4, 7 which correspond to the minimum and maximum of the applied force f.

Figure 3 – The distribution of the stress σ Mises on the trapezoidal plate with hole at t = 1,5 and

t = 4, 7 which correspond to the minimum and maximum of the applied force f given as in Figure 2.

29


Figure 4 – The component of the stress σ 11 versus the strain ε 11 at the material points given by node 9 obtained using the Cleja-Tigoiu & Stoicuta solution and the Simo-Hughes solution.

Figure 5 – The component of the stress σ 22 versus the strain ε 22 at the material points given by node 9, obtained using the Cleja-Tigoiu & Stoicuta solution and the Simo-Hughes solution.



30




31


p plotted as Figure 10 – The graphs of the components of the strain ε 33 and of the plastic strain ε 33 functions of time at the material points given by node 9, obtained using the Cleja-Tigoiu & Stoicuta solution and the Simo-Hughes solution.

6. CONCLUSIONS

Concerning the comparison of the numerical results using the Simo-Hughes algorithm and the in-plane algorithm, we note that the material responses are completely different on one side for the considered methods, and the other hand for the points located in the vicinity of the loading boundary, in Figures 4-6, and points close to the free boundary, Figures 7-9. Figures 4-6 show that the stresses take comparable values, while the strains are larger for the Simo-Hughes algorithm than in plane algorithm. Figures 7-9 show large differences between the components of the stress, especially for the shear and normal component in the perpendicular direction of the applied force. Finally we remark that the Chaboche isotropic hardening low, (5), leads to the existence of a limit cycle, which characterizes the evolution of the curves stress versus strain after a large number of periods of the applied force. The linear isotropic hardening low, (4) leads to a stress threshold and limit strain after a large numbers of period of the applied force.

Acknowledgement. The first author gratefully acknowledge the support from the University of Bucharest under Grant UB, Contract No. 339

Received on December 12, 2013


32

REFERENCES 1. ARMSTRONG, P.J., FREDERICK, C.O., A mathematical representation of the multiaxial Bauschinger effect, G.E.G.B. Report RD-B-N 731, 1996. 2. BATHE, J., Finite element procedure Vol. I, II, III, Prentice Hall, Engl. Cliffs, New Jersey, 1996. 3. BELYTSCHKO, T., WING, K.L., MORAN B., Nonlinear Finite Elements for Continua and Structures, British Library Cataloguing in Publication Data, Toronto, 2006. 4. BROGGIATO, G.B. CAMPANA, F. CORTESE, L., The Chaboche nonlinear kinematic hardening model: calibration methodology and validation, Meccanica, 43, 2, pp. 115-124, 2008. 5. CHABOCHE, J.L., Constitutive equations for cyclic plasticity and cyclic viscoplasticity, Int. J. Plasticity, 5, 3, pp. 247-302, 1989. 6. CLEJA-TIGOIU, S., CRISTESCU, N., Teoria plasticitatii cu aplicatii la prelucrarea materialelor, Ed. Univ. Bucuresti, 1985. 7. CLEJA-TIGOIU, S., Anisotropic Elasto-plastic Model for Large Metal Forming Deformation Processes, International Journal of Forming Processes, 10, 1, pp. 67-89, 2007. 8. CLEJA-TIGOIU, S., STOICUTA, N., Revisited Simo algorithm for the plane stress state, Applied Mathematics and Computation, (in press), 2013. 9. FISH, J., BELYTSCHKO, T., A First Course in Finite Elements, John Wiley&Sons Ltd., USA, 2007. 10. HAN, W., REDDY, B.D., Plasticity. Mathematical Theory and Numerical Analysis, Springer, New York, 1999. 11. HILL, R., The Mathematical Theory of Plasticity, Oxford University Press, Oxford, U.K., 1950. 12. HUGHES, T.J.R., Numerical Implementation of Constitutive Models: Rate- Independent Deviatoric Plasticity, chap.2 In: “Theoretical Foundations for Large Scale Computations of Nonlinear Material Behavior” (eds.: S. Nemat-Nasser, R. Asaro, G. Hegemier), pp. 29-57, Martinus Nijhoff Publishers, Dordrecht, The Netherlands, 1984. 13. HUGHES, T.J.R., The Finite Element Method. Linear Static and Dynamic Finite Element Analysis, Prentice-Hall, Englewood, New Jersey,1987. 14. JOHNSON, C., Numerical solution of partial differential equation by the finite element method, Cambridge University Press, Cambridge, 2002. 15. KACHANOV, L.M., Fundamentals of the Theory of Plasticity, Mir Publishers, Moscow, 1974. 16. KHAN, A.S., HUANG, S., Continuum theory of plasticity. Fundamentals of the Theory of Plasticity, John Wiley & Sons, INC., 1995 17. KRIEG, R.D., KEY, S.W., Implementation of a Time Dependent Plasticity Theory into Structural Computer Programs. Constitutive Equations in Visco-plasticity, Computational and Engineering Aspects, eds.: J.A. Stricklin and K.J. Saczalski, AMD-20, ASME, New York, 1976. 18. LORET, B., PREVOST, J.H., Accurate Numerical Solutions for Drucker- Prager Elastic-Plastic Models, Computer Methods in Applied Mechanics and Engineering, 54, pp. 259-277, 1986. 19. LUBARDA, V.A., BENSON, D.J., On the numerical algorithm for isotropic - kinematic hardening with the Armstrong - Frederick evolution of the back stress, Comput. Methods Appl. Mech. Engrg., 191, pp. 3583-3596, 2002. 20. MAENCHEN, G., SACK, S., The Tensor Code, In: “Methods in Computational Physics”, vol 3., ed.: B. ALDER, Academic Press, New York, pp. 211, 1964. 21. NAGTEGGAL, J.C., On the Implementation of Inelastic Constitutive Equations with Special Reference to Large Deformation problems, Comp. Meth. Appl. Mech. Eng., 33, pp. 469, 1982. 22. ORTEGA, J., RHEINBOLDTW, C., Iterative solution of nonlinear equations in several variables, Acad. Press, New York, 1970. 23. ORTIZ, M., POPOV, E.P., Accuracy and Stability of Integration Algorithms for Elastoplastic Constitutive Equations, Int. Journ. for Num. Meth. in Eng., 21, pp. 1561-1576, 1985.

33


24. PARASCHIV-MUNTEANU, I., CLEJA-TIGOIU, S., SOOS, E., Plasticitate cu aplicatii în geomecanică, Ed. Univ. Bucuresti, 2004. 25. PRAGER, W., A new method of analysing stress and strain in work-hardening plastic solid, Journal Appl. Mech., 23, pp.493-496, 1956. 26. SCHREYER, H.L., KULAK, R. L., KRAMER, J.M., Accurate Numerical solutions for Elastoplastic models, Journal of Pressure Vessel Technology, 101, pp. 226-334, 1979. 27. SIMO, J.C., TAYLOR, R.L., Return Mapping Algorithm for Plane Stress Elastoplasticity, International Journal for Numerical Methods in Engineering, 22, pp. 649-670, 1986. 28. SIMO, J.C., HUGHES, T.J.R., Computational Inelasticity, Springer, New-York, 1998. 29. VOCE, E., A Practical Strain Hardening Function, Metallurgica, 51, pp. 219-226, 1955. 30. ZIEGLER, H., A Modification of Prager’s Hardening Rule, Quarterly of Applied Mathematics, 17, pp. 55-65, 1959. 31. WILKINS, M.L., Calculation of Elastic Plastic Flow, In: “Methods in Computational Physics”, vol 3., ed.: B. ALDER, Academic Press, New York, pp. 211, 1964.

DYNAMICS MODELLING OF A PARALLEL FLIGHT SIMULATOR STEFAN STAICU *, CONSTANTIN OCNARESCU, LIVIU UNGUREANU

Abstract. Recursive matrix relations for dynamics analysis of a parallel manipulator, namely the spatial flight simulator, is established in this paper. Knowing the general motion of the platform, the inverse dynamics problem is solved using an approach based on explicit equations of parallel robots dynamics. Finally, some simulation graphs for the input forces and powers are obtained. Key words: connectivity relations, dynamics, flight simulator, parallel robot.

1. INTRODUCTION

Parallel robots are closed-loop structures presenting very good potential in terms of accuracy, stiffness and ability to manipulate large loads. One of the main bodies of the mechanism is fixed and is called the base, while the other is regarded as movable and hence is called the moving platform of the manipulator. The links of the robot are connected one to the other by spherical joints, universal joints, revolute joints or prismatic joints. Generally, the number of actuators is typically equal to the number of degrees of freedom and each leg is controlled at or near the fixed base [1-3]. These mechanisms can be found in practical applications, in which it is desired to orient a rigid body in space of high speed, such as flight simulators [4-6] and positional trackers [7,8]. A first spatial robot used in the amusement parks was built in 1928 by James E. Gwinnett [9] and was patented as a simulator by Ga1anov and Diaciun [10]. In the present paper, a recursive matrix method is applied to the inverse dynamics analysis of a spatial 3-DOF parallel mechanism, to prove that the number of equations and computational operations reduces significantly by using a set of matrices for dynamics modelling.

* “Politehnica” University of Bucharest, Department of Mechanics, Department of Mechanisms Theory and Robotics, Romania


Stefan Staicu, Constantin Ocnarescu, Liviu Ungureanu

2

2. KINEMATICS ANALYSIS

With a closed-loop structure, the parallel manipulator has a spatial construction and is used as a flying simulator in the amusement parks. Additionally, a closed cab with two doors is placed on the moving platform (Fig.1). The architecture of the robot consists of a fixed triangular base A1 B1C1 and an upper moving platform A4GB4C4 that is a rectangle with b and h the lengths of the geometric characteristics, two passive legs and three active extensible legs. The first active leg A is typically contained within the Ox0 y0 median vertical plane of the cab platform, whereas the remaining active legs B, C are starting from a perpendicular vertical plane, symmetrically placed relatively to the median plane. Provided with identical kinematical structure, these last active limbs connect the fixed base to the moving platform by two universal (U) joints interconnected through a prismatic (P) joint made up of a cylinder and a piston. Hydraulic or pneumatic systems can be used to vary the lengths of the prismatic joints and to control the location of the platform. There are three active prismatic joints A2 , B3 , C3 , four passive universal joints B1 , B4 , C1 , C4 , seven passive revolute joints A1 , A3 , A4 , D1 , D2 , E1 , E2 and one internal passive prismatic pair jointed at the point E2 .

Fig. 1 – Parallel flight simulator.

3

Dynamics modelling of a parallel flight simulator

For the purpose of analysis, we assign a fixed Cartesian coordinate system Ox0 y0 z0 (T0 ) at the point O of the fixed base platform and a mobile frame GxG yG zG attached on the mobile platform at its characteristic point G. Grübler mobility equation predicts that the device has certainly three degrees of freedom. The moving platform is initially located at a central configuration, where the platform is not rotated with respect to the fixed base and the point G is located at an elevation hG . We also consider that the three sliders A2 , B3 , C3 are initially starting from the positions A1 A2 = s A , B2 B3 = sB , C2C3 = sC (Fig. 2). To simplify the graphical image of the kinematical scheme of the mechanism, in what follows we will represent the intermediate reference systems by only two axes, so as is used in most of robotics papers [1,2,11]. It is noted that the relative rotation with angle ϕk ,k −1 or the relative translation of the body Tk with the displacement λ k ,k −1 must always be pointed along the direction of the zk axis.

Fig. 2 – Kinematical scheme of parallel mechanism.

The leg A consists of a fixed revolute joint, a moving cylinder of length l1 , mass m1A and tensor of inertia Jˆ1A , which has a rotation about z1A axis with the A A A A A 10 angle ϕ10 , the angular velocity ω10 and the angular acceleration ε10 . = ϕ 10 = ϕ


4

A prismatic joint is as well as a piston linked at the A2 x2A y2A z2A frame, having a A A A relative motion with the displacement λ 21 , the velocity v21 and the = λ 21 A A A  . It has the length l , mass m and tensor of inertia Jˆ A . A acceleration γ = λ 21

21

2

2

2

m3A

new revolute joint is introduced at a homogenous rod of length h, mass and tensor of inertia Jˆ3A , linked at the A3 x3A y3A z3A frame, having a perpendicular A A A rotation about z3A axis with the angle ϕ32 , the angular velocity ω32 and the = ϕ 32 A A  32 angular acceleration ε32 . Finally, the median axis of the moving platform of = ϕ mass m A and tensor of inertia Jˆ A rotates around the rod A G with the angle ϕ A , 4

4

A ω43

43

3

A = ϕ 43

A ε32

A  32 . = ϕ

the angular velocity and the angular acceleration Concerning other two active kinematical chains, the architecture of the leg B, for example, consists of the cross of a fixed Hooke joint characterised by a negligible mass, linked at the frame B1 x1B y1B z1B , which has the angle of rotation B B B , the absolute angular velocity ω10 and the angular acceleration ϕ10 = ϕ 10 B B 10 . This is connected at a moving cylinder B2 x2B y2B z2B of length l3 , mass ε10 = ϕ m2B and tensor of inertia Jˆ2B , having a relative rotation around B2 z2B axis with the

B B B B B  21 angle ϕ21 , so that ω21 , ε 21 . An actuated prismatic joint is as well as a = ϕ = ϕ 21 B piston of length l , mass m and tensor of inertia Jˆ B , linked to the B x B y B z B 4

3

3

B , λ32

3 3

3 3

B λ 32

frame, having a relative displacement the velocity and the = B  B . Finally, a second universal joint B is introduced at the acceleration γ 32 = λ 4 32 apex of the rectangular moving platform. Additionally, two rigid cranks of lengths l5 , l6 , masses m1D , m1E and tensors of inertia Jˆ D , Jˆ E , respectively, jointed at two fixed revolute pairs D , E , execute 1

B v32

1

1

oscillating rotations with the angles

D E , ϕ10 , ϕ10

the angular velocities

1 D E ω10 , ω10

and

D E . ω10 , ω10

the angular accelerations Finally, the internal prismatic pair E2 of negligible mass starts from an initial position D2 E2 = s and has the relative  . displacement λ , the velocity v = λ and the acceleration γ = λ Starting from the reference origin O and pursuing along five independent legs OA1 A2 A3 A4 , OB1B2 B3 , OC1C2C3 , OD1 , OE1 , we obtain following transformation matrices

5


ϕ α ϕ α ϕ T a10 = a10 a1 , a21 = θ , a32 = a32 a2 θ, a43 = a43 θ , ϕ ϕ α b10 = b10 , b21 = b21 a3 θ, b32 = θ,

(1)

ϕ ϕ α T ϕ α ϕ α c10 = c10 , c21 = c21 a3 θ , c32 = θ , d10 = d10 a4 , e10 = e10 a5 ,

where we denote the matrices: ϕ = aiα rot( z , αi ) , θ =rot( y, π ) , p= k , k −1 rot( z , ϕk , k −1 ) 2

( p = a , b, c , d , e ) .

(2)

A complete description of the absolute position and orientation of the platform requires four variables: the coordinates x0G , y0G of the point G and two Euler’s angles ψ, φ associated with two successive orthogonal rotations. Since all rotations take place successively about the moving coordinate axes, the general rotation matrix R = R21R10 is obtained by multiplying two known transformation matrices = R10 rot( z , ψ + α1 − α 2 ), = R21 rot( x, φ) .

(3)

Here, x0G , y0G and φ are chosen as the independent variables and ψ is a parameter of a parasitic motion. The parasitic motion from the four motions of the moving platform is permanently dependent on the independent variables. The conditions concerning the absolute orientation of the moving platform are given by the following matrix identity θT a40 = R,

(4)

from which we obtain significant relations between the angles of rotation A A A . ψ = ϕ10 − ϕ32 , φ = ϕ43

(5)

In the inverse geometric problem, the position of the mechanism is completely given through three variables x0G , y0G , φ . For the parallel simulator, the rotation angle φ is the only completely independent variable with other three pose parameters. Consider, for example, that during three seconds the rotation of the moving platform around its median axis and the motion of the characteristic point G along a rectilinear trajectory are expressed in the fixed frame Ox0 y0 z0 through the following analytical functions x0G − (d 2 + d3 + d 4 ) x0G∗

=

y0G − hG y0G∗

=

φ

π = 1 − cos t , 3 φ ∗

(6)


6

where the values 2 x0G* , 2 y0G* , 2φ∗ denote the final position and orientation of the moving platform. A

set

of

12

independent

variables

A A A B B B , ϕ10 , ϕ10 , λ 21 , ϕ32 , ϕ21 , λ32

C D E , ϕ10 ϕ10 , ϕC21 , λC32 , ϕ10 , λ will be determined by vector-loop equations

 r10A +

3

∑

 akT0 rkA+1,k

B 2 T B T  B4 T  GB4 r10 + bk 0 rk +1,k + b30 r3 − a40 r4 = =

∑

k =1

k =1

 = r10C +

T G + a40 r4

2





∑ ckT0 rkC+1,k + c30T r3C k =1

4

(7)  T  D2 T  D2G T  GC4 r10D + d10 r1 + a30 r3 = − a40 r4 =

  T  E2 T  E2G = r10E + e10 r1 + a30 r3 = r0G ,

where A  A A   A  G  A  )u1 , r32 r10 = 0, r21 = ( s A + λ 21 = l2u3 , r43 = 0, r4= hu3       1  B   B4 r10B = [d 2 + d3 + d 4 0 d5 ]T , r21B = 0, r32B = ( sB + λ32 )u1 , r3 = l4u3 , r4GB4 = − bu1 2  C        1 C C 4 4 r10= [d 2 + d3 + d 4 0 − d5 ]T , r21 = 0, r32 = ( sC + λC32 )u1 , r3C= l4u3 , r4GC= bu1 2 D   D2   D2G  (d 2 + d3 )u1 , r1 = (d1 − h)u1 r10 = l5u1 , r3 =       r10E= d 2u1 , r1E2= l6u1 , r3E2G= ( s + d1 − h + λ)u1

(8)

  u1 = [1 0 0]T , u3 = [0 0 1]T

From the vector equations (7) we obtain 12 analytical equations for the inverse geometric solution, giving the expressions of all 12 above independent variables. The motions of the component elements of the manipulator are characterized by the relative velocities of the joints   vk ,k −1 = λ k ,k −1u3 (9) or by the relative angular velocities and their associated skew-symmetric matrices    k ,k −1 = ωk .k −1 = ϕ k ,k −1u3 , ω ϕ k ,k −1u3 . (10) Deriving the relations (5) and the geometrical constraints (7), we obtain first the angular velocity ω A = φ and 12 matrix conditions of connectivity [12]: 43


7

  V = [Q]−1 P ,

(11)

where followings nonzero terms determines the contents of 12×12 invertible square  matrix [Q] and the column matrix P :  T   T  T T A T T T G T G = qiA1 uiT a10 u3{r21A + a21 r32 + a21 a32 a43 r4 } , qiA2 = uiT a10 u1 , qiA3 = uiT a30 u3 a43 r4     T T T T GB4 T T GB4 q Bj1 = −u Tj a10 u3 a21 a32 a43 r4 , q Bj3 = −u Tj a30 u3 a43 r4 ,

 T  T  T  T B T  B4 T  B4 q Bj4 = u Tj b10 u1 u3b21 {r32 + b32 r3 } , q Bj5 = u Tj b20 u3{r32B + b32 r3 } , q Bj6 = u Tj b20  T  T T T T  GC4 T  GC4 qCj1 = −u Tj a10 u3 a21 a32 a43 r4 , qCj3 = −u Tj a30 u3 a43 r4 ,  T  T  T  T C T  C4 T  C4 qCj7 = u Tj c10 u1 u3c21 {r32 + c32 r3 } , qCj8 = u Tj c20 u3{r32C + c32 r3 } , qCj9 = u Tj c20

 T  D2  T  T  D2G T T  D2G u3 r1 qiD1 = uiT a10 u3 a21 a32 r3 , qiD3 = uiT a30 u3 r3 , qiD10 = uiT d10  T  T  E2G  T  E2  T T T  E2G qiE1 = uiT a10 u3 a21 a32 r3 , qiE3 = uiT a30 u3 r1 , qiE12 = uiT a30 u3 r3 , qiE11 = uiT e10 u1    T  GB4    T  GC4   p Bj u Tj r0G + φ u Tj a40 u3 r4 , = pCj u Tj r0G + φ u Tj a40 u3 r4 piA = uiT r0G , =   piD = uiT r0G ,

  piE = uiT r0G

(i = 1, 2) ( j = 1, 2, 3)

(12)

Finally, the expressions of relative velocities are obtained from the column matrix  D A A A C C E B B B (13) V = [ω10 ω10 ωC21 v32 v21 ω32 ω10 ω10 v]T . v32 ω10 ω21 Considering some independent virtual motions of the spatial mechanism, virtual displacements and velocities should be compatible with the virtual motions imposed by all kinematical constraints and joints at a given instant in time. Let us assume that the robot has successively three virtual motions determined by following sets of velocities: Bv Cv Av v21 a = 1 , v32 a = 0 , v32 a = 0 Av Bv Cv v21 b = 0 , v32b = 1 , v32b = 0 Av v21 c

Bv = 0 , v32 c

Cv = 0 , v32 c

(14)

=1.

The characteristic virtual velocities are expressed as functions of the pose of the mechanism by the general kinematical equations (11). Expressions of relative accelerations are obtained from the column matrix  B B B C C D A A A E (15) Γ = [ε10 ε10 ε 21 γ 32 ε10 εC21 γ 32 γ 21 ε32 ε10 ε10 γ ]T .


8

using new conditions of connectivity:   Γ =[Q]−1 S ,

(16)

   where following terms determine the contents of column matrix S= P − [Q ]V :    A A T T T A T T T G A A T T T G ω10ui a10u3u3{r21A + a21 r32 + a21 a32 a43 r4 } − ω32 siA = uiT r0G − ω10 ω32ui a30u3u3a43 r4 −  A A T T T T T G A A T T −2ω10 v21ui a10u3u1 − 2ω10 ω32ui a10u3 a21a32u3 a43 r4 }   u T a T u r GB4 + ω A ω A u T a T u u a T a T a T r GB4 + s Bj = u Tj r0G + φ j 40 3 4 10 10 j 10 3 3 21 32 43 4    T  GB A A T T T +ω32 ω32u j a30u3u3 a43 r4 4 + φ 2u Tj a40 u3u3 r4GB4 +  A A T T T T T  GB4 T T T A  T T +2ω10 ω32u j a10u3 a21 a32u3 a43 r4 + 2ω10 φu j a10u3 a21 a32 a43u3 r4GB4 +  B B T T T B T  B4 T A  T T ω10u j b10 u3u3b21 r3 } − {r32 + b32 u3 r4GB4 − ω10 +2ω32 φui a30u3 a43   B B T T B B T T T  B4 T  B4 T −ω21 ω21u j b20 u3u3{r32B + b32 ω21u j b10 u3b21 r3 } − 2ω10 u3{r32B + b32 r3 } −     B B T T B B T T T −2ω10 v32u j b10 u3b21 v32u j b20 u3u1 u1 − 2ω21   u T a T u r GC4 + ω A ω A u T a T u u a T a T a T r GC4 + s Cj = u Tj r0G + φ j 40 3 4 10 10 j 10 3 3 21 32 43 4  T    GC T A A T T u3u3 r4GC4 + +ω32 ω32u j a30u3u3 a43 r4 4 + φ 2u Tj a40  A A T T T T T  GC4 T T T A  T T +2ω10 ω32u j a10u3 a21 a32u3 a43 r4 + 2ω10 φu j a10u3 a21 a32 a43u3 r4GC4 +  C C T T T T C T  C4 A  T T ω10u j c10 u3u3c21 +2ω32 φui a30u3 a43 u3 r4GC4 − ω10 {r32 + c32 r3 } −  T   C C T T T  C4 T  C4 T ω21u j c10 u3c21 −ωC21ωC21u Tj c20 u3u3{r32C + c32 r3 } − 2ω10 u3{r32C + c32 r3 } −     C C T T C T T T −2ω10 v32u j c10 u3c21 u j c20 u3u1 u1 − 2ωC21v32    T T  D2G A A T T A A T T −ω32 ω32ui a30u3u3 r3D2G − ω10ui a10u3u3 a21 a32 r3 siD = u Tj r0G − ω10   A A T T D D T T T T ω32ui a10u3 a21 a32u3 r3D2G −ω10 ω10ui d10u3u3 r1D2 − 2ω10    T T  E2G A A T T A A T T ω10ui a10u3u3 a21 a32 r3 −ω32 ω32ui a30u3u3 r3E2G − siE = u Tj r0G − ω10   A A T T T T E E T T −ω10 ω10ui e10u3u3 r1E2 − 2ω10 ω32ui a10u3 a21 a32u3 r3E2G −  A T T A T T T T  −2ω10 vui a10u3 a21 a32u1 − −2ω32 vui a30u3u1 (i = 1, 2) ( j = 1, 2, 3) .

(17)


9

3. INVERSE DYNAMICS MODEL

The dynamics analysis of parallel robots is complicated because the existence of a spatial kinematical structure, which possesses a large number of passive degrees of freedom, dominance of the inertial forces, frictional and gravitational components and by the problem linked to real-time control in the inverse dynamics. Considering all gravitational effects and neglecting the frictions forces, the relevant objective of the inverse dynamics is to determine the input torques or forces, which must be exerted by the actuators in order to produce a given trajectory of the end-effector. A lot of works have focused on the dynamics of Stewart platform. Dasgupta and Mruthyunjaya [13] used the Newton-Euler approach to develop closed-form dynamic equations of Stewart platform, considering all dynamic and gravity effects as well as viscous friction at joints. Tsai and Stamper [14] presented an algorithm to solve the inverse dynamics for the Delta translational robot, using Newton-Euler equations and Lagrange formalism. Three independent pneumatic or hydraulic systems A, B, C that generate three C     C input forces f 21A = f 21Au3 , f32B = f32B u3 , f32 = f32 u3 , which are oriented along the axes A2 z2A , B3 z3B , C3 z3C , control the motion of three moving pistons of the legs. The parallel robot can artificially be transformed in a set of three open chains Ci (i = A, B, C ) subject to the constraints. This is possible by cutting each joint for moving platform and taking its effect into account by introducing the corresponding constraint conditions. The force of inertia and the resulting moment of inertia forces of a rigid body TkA , for example,   A  A  kA0 ω  kA0 + ε kA0 rkCA  , FkinA 0 = − mk  γ k 0 + ω   inA   A CA  A A A A A  k 0 Jˆk ωkA0 ] M k 0 = − [mk rk γ k 0 + Jˆk εk 0 + ω

(

)

(18)

are determined with respect to the centre of joint Ak . On the other hand, the   wrench of two vectors Fk∗ A and M k∗ A evaluates the influence of the action of the  weight mkA g and of other external and internal forces applied to the same element TkA of the manipulator, for example:     Fk* A = mkA gak 0u3 , M k* A = mkA grkCA ak 0u3

(k = 1, 2, 3, 4) .

Pursuing the first leg A, two recursive relations generate the vectors A A      F= Fk 0 + akT+1,k FkA+1 , M kA = M kA0 + akT+1,k M kA+1 + rkA+1,k akT+1,k FkA+1 , k where one denoted

(19)

(20)


      − FkinA − Fk∗ A , M kA0 = FkA0 = − M kinA − M k∗ A .

As example, starting from (20), we develop a set of recursive relations:   A A  A A  A A  T T T A F4A = F40A , F= F30 + a43 F4A , F= F20 + a32 F3A , F= F10 + a21 F2 3 2 1  A  A   T T , M 3A = M 4A = M 40 M 30 + a43 M 4A + r43A a43 F4A  A    A   T T T T A M10 + a21 M 2A + r21A a21 F2 . M 2A = M 20 + a32 M 3A + r32A a32 F3A , M1A =

10

(21)

(22)

The fundamental principle of the virtual work states that a mechanism is under dynamic equilibrium if and only if the total virtual work developed by all external, internal and inertia forces vanish during any general virtual displacement, which is compatible with the constraints imposed on the mechanism. Applying the explicit form of equations of the parallel robots dynamics [15], compact matrix relations result for the input force of three prismatic actuators A A A B   Av Av Av Bv = f 21A u3T {F2A + ω10 M + ω M + ω M + ω M a 1 32 a 3 43a 4 10 a 1 +     E Bv B Cv C Cv C Dv D Ev + ω21 a M 2 + ω10 a M 1 + ω21a M 2 + ω10 a M 1 + ω10 a M 1 } B B C C   Bv Bv Cv Cv = f32B u3T {F3B + ω10 b M 1 + ω21b M 2 + ω10b M 1 + ω21b M 2 + (23) A A A D E Av Av Av Dv Ev + ω10 b M 1 + ω32 a M 3 + ω43b M 4 + ω10b M 1 + ω10 a M 1 } C C B B   C Cv Cv Bv Bv = f32 u3T {F3C + ω10 c M 1 + ω21c M 2 + ω10 c M 1 + ω21c M 2 + A A A D E Av Av Av Dv Ev + ω10 c M 1 + ω32 c M 3 + ω43c M 4 + ω10 c M 1 + ω10 c M 1 }. The relations (20)-(23) represent the inverse dynamics model of the parallel flight simulator. The various dynamical effects, including the Coriolis and centrifugal forces coupling and the gravitational actions are considered in this explicit equation. As an application let us consider a parallel simulator which has the following geometrical and architectural characteristics x0G* =0.05m , y0G* =0.1m , φ∗ =π /12 A3 D= d= 2 1 0.4 m, A1 E= 1 d= 2 0.95m, E1 D= 1 d= 3 0.65m D1F = d= = C2 F = d= 4 0.25m, B2 F 5 1.35m D1D2= l5= 1.3m, E1E2= l6= 1.2 m, B4C4= b= 1.2 m A4G = h = 1.4 m, FG = hG = 1.1m, ∆t = 3 s s= s= 0.33m A 0.29 m, s= B C

= l1 0.65m, = l2 0.85m, = l3 0.75m, = l4 1m

11


= m1A 1= kg , m2A 1.2= kg , m3A 1.5= kg , m4A 5kg B C B C D E m= m= m= 1kg , m= 1.2kg , m= 1.5kg , m= 1.4kg . 2 2 3 3 1 1

Using MATLAB software, a computer program was developed to solve the dynamics of the parallel simulator. To develop the algorithm, it is assumed that the platform starts at rest from a central configuration and moves pursuing successively some evolutions. Three examples are solved to illustrate the algorithm. In a first example, the platform’s point G moves along the vertical direction Oy0 with variable acceleration while all the other positional parameters are held equal to zero. As can be seen from Fig.3 and Fig.4 it is proved to be true that the active forces and powers of the actuators B and C are permanently equal to one another. For the case when the point G moves along a rectilinear trajectory without a rotation of the platform around its median axis, the graphs are illustrated in Fig.5 and Fig.6.

Fig. 3 – Input forces of three actuators.

Fig. 4 – Powers of three actuators.




12

Further on, we shall consider a general evolution of the platform, combining a rectilinear displacement of the point G and a rotation of the platform around its median axis. The input forces and powers of three prismatic actuators are graphically sketched in Fig.7 and Fig.8.



4. CONCLUSIONS

Some exact relations that give in real-time the position, velocity and acceleration of each element of a flight parallel simulator have been established in the present paper. The dynamics models take into consideration the masses and forces of inertia introduced by all component elements of the parallel mechanism. Based on explicit equations of parallel robots dynamics, the approach can eliminate all forces of internal joints and establishes a direct determination of the timehistory evolution of active forces and powers required by the actuators. The simulation certifies that one of the major advantages of the current matrix recursive formulation is the accuracy and a smaller processing time for the numerical computation. Choosing the appropriate serial kinematical circuits connecting many moving platforms, the present method can be easily applied in forward and inverse mechanics of various types of parallel mechanisms, complex manipulators of higher degrees of freedom and particularly hybrid structures, with increased number of components of the mechanisms.

Received on December 2, 2013

13


REFERENCES 1. TSAI, L-W., Robot analysis: the mechanics of serial and parallel manipulators, Wiley, 1999. 2. PLITEA, N., LESE, D., PISLA, D., VAIDA, C., Structural design and kinematics of a new parallel reconfigurable robot, Robotics and Computer-Integrated Manufacturing, 29, 1, pp. 219-235, 2013. 3. PISLA, D., ITUL, T., PISLA, A., GHERMAN, B., Dynamics of a Parallel Platform for Helicopter Flight Simulation Considering Friction, SYROM 2009, pp 365-378, Springer, 2009. 4. MERLET, J-P., Parallel robots, Kluwer Academic, 2000. 5. G. GOGU, T2R1-type parallel manipulators with bifurcated planar-spatial motion, European Journal of Mechanics, A/Solids, 33, pp. 1-11, 2012. 6. STEWART, D., A Platform with Six Degrees of Freedom, Proc. Inst. Mech. Eng., 180, 1, pp. 371386, 1965. 7. DI GREGORIO, R., PARENTI CASTELLI, V., Dynamics of a class of parallel wrists, ASME Journal of Mechanical Design, 126, 3, pp. 436-441, 2004. 8. PISLA, A., ITUL, T., PISLA, D., SZILAGHYI, A., Considerations upon the influence of manufacturing and assembly errors on the kinematic and dynamic behavior in a flight simulator Stewart-Gough platform, Mechanisms, Transmissions and Applications: Mechanisms and Machine Science, Volume 3, Springer, pp 215-223, 2012. 9. GWINNETT, J.E., Amusement devices, US Patent No. 1789680, 1931. 10. GALANOV, A.S., DIACIUN, V.C., Manipulator, URSS Patent No. 908589, 1982. 11. ANGELES, J., Fundamentals of Robotic Mechanical Systems: Theory, Methods and Algorithms, Springer, 2002. 12. STAICU, S., Méthodes matricielles en cinématique des mécanismes, UPB Scientific Bulletin, Series D: Mechanical Engineering, 62, 1, pp. 3-10, 2000. 13. DASGUPTA, B., MRUTHYUNJAYA, T.S., A Newton-Euler formulation for the inverse dynamics of the Stewart platform manipulator, Mechanism and Machine Theory, 34, pp. 711725, 1998. 14. TSAI, L-W., STAMPER, R., A parallel manipulator with only translational degrees of freedom, ASME Design Engineering Technical Conferences, Irvine, CA, 1996. 15. STAICU, S., LIU, X-J., LI, J., Explicit dynamics equations of the constrained robotic systems, Nonlinear Dynamics, 58, 1-2, pp. 217-235, 2009.

FORMULATION ANALYTIQUE DES EFFORTS DE REDUCTION AU DROIT DE L’INTERFACE FLUIDE-STRUCTURE DES BARRAGES RIGIDES A FRUITS COMPOSES SOUS EXCITATIONS SISMIQUES ABDELMADJID TADJADIT *, BOUALEM TILIOUINE*

Résumé. Les barrages, constituent des ouvrages à caractère stratégique important. L’évaluation de leur sécurité parasismique requiert, entre autres, une détermination précise des pressions hydrodynamiques et par suite de la distribution des efforts agissant au droit de l’interface fluide-structure. Cette tâche devient plus complexe en présence d’un parement amont à fruit irrégulier. Dans le présent article, une solution semi-analytique est présentée pour le calcul des efforts agissant sur des barrages rigides à parements irréguliers sous excitations sismiques. Les résultats numériques obtenus sont présentés pour différentes configurations géométriques des parements amont. A défaut de solutions exactes pour le cas d’un barrage rigide à double fruit, les résultats numériques obtenus sont comparés à ceux issus respectivement de l’application de la méthode exacte de Westergaard pour le cas particulier d’un barrage rigide à parement vertical, et à ceux de Chwang applicable uniquement pour le cas d’un fruit simple. Une très bonne concordance est observée et des conclusions d’intérêt pratique sont formulées. Mots clés: barrage rigide à double fruit, interface fluide-structure, excitation sismique, pressions hydrodynamiques, méthode semi-analytique, efforts de réduction.

1. INTRODUCTION

Les barrages de part leur importance, constituent des ouvrages stratégiques. Leur étude relève d’un processus complexe prenant en compte plusieurs phénomènes nécessitant la conjugaison de connaissances de plusieurs disciplines. L’évaluation de leur sécurité requiert un niveau élevé de précision car les dommages à encourir, particulièrement en présence de séismes, sont d’une valeur inestimable, non seulement en termes de vies humaines mais aussi en termes de retombées économiques. De ce point de vue, il est clair que la stabilité de tels ouvrages impose une investigation profonde pour une évaluation aussi réaliste que

*

Ecole Nationale Polytechnique, Laboratoire de Génie Sismique et de Dynamique des Structures. 10 Avenue Hassen Badi, El-Harrach 16200, Alger, Algérie. E-mails: [email protected], [email protected]


Abdelmadjid Tadjadit, Boualem Tiliouine

2

possible de l’ensemble des paramètres susceptible d’influencer directement ou indirectement sur le comportement du système global barrage-réservoir. Dans ce contexte, plusieurs travaux de recherches ont été effectués et Westergaard [13] fût le pionnier à avoir formulé une expression analytique permettant l’évaluation de la pression hydrodynamique pour le cas typique d’un barrage rigide vertical soumis à un mouvement harmonique du sol dans la direction horizontale. Au fil des années, l’étude classique de Westergaard fût enrichie davantage en intégrant plusieurs paramètres physiques, tels que la compressibilité du fluide dans le réservoir [2,5,11], la flexibilité du barrage [12] et l’absorption du fond du réservoir [4,6,8]. A côté des méthodes numériques, on retrouve également les méthodes analytiques. Leur domaine d’application se limite pratiquement pour des barrages présentant des interfaces fluide-structure à géométries simples. Toutefois, les solutions formulées, bien qu’approchées, sont souvent simples et pratiques à utiliser [9,10]. Pour l’ingénieur concepteur, le calcul de la résultante des efforts tranchants générés par les pressions hydrodynamiques et l’évaluation des moments de flexion associés constitue une tâche incontournable pour le dimensionnement de la structure du barrage. De nos jours, il existe plusieurs méthodes de calcul notamment celles relevant des calculs numériques de pointe. Toutefois, lorsque la situation le permet, les méthodes analytiques ou plutôt semi-analytiques restent d’actualité et sont d’un apport non négligeable. L’objectif de la présente étude est de formuler analytiquement les efforts de réduction agissant au droit de l’interface fluide-structure des barrages rigides à double fruit. A cet effet, un programme de calcul numérique sous langage Matlab7 [7] a été établi, ce qui a permis de valider plusieurs modèles de barrages rigides pour diverses configurations géométriques de la face amont. A défaut de solutions exactes pour le cas des barrages rigides à double fruit, les résultats numériques obtenus sont comparés à ceux issus respectivement de l’application de la méthode exacte de Westergaard [13] pour le cas particulier d’un barrage rigide à parement vertical et à ceux de la méthode de Chwang [3] applicable uniquement pour le cas d’un barrage à fruit simple. Après validation des résultats obtenus, l’étude est ensuite étendue aux cas des barrages rigides à double fruit. 2. FORMULATION THEORIQUE DU PROBLEME 2.1. Modèle Mathématique

On considère un barrage infiniment rigide à double fruit. L’origine des coordonnées est située à la base du réservoir et la surface libre d’eau est désignée par y = H (Fig.1).

3

Formulation analytique des efforts de réduction au droit de l’interface fluide-structure (...)

Y

Surface libre d'eau

s

3

s

2

n

Réservoir

s

4

H CH

Barage rigide

y

s

1

O

s

5

X

Fig. 1 – Barrage rigide à double fruit avec réservoir de longueur infinie soumis à un mouvement horizontal du sol.

Soit un point M (x, y) dans le repère (O, x, y) situé sur la partie amont du barrage à une élévation quelconque (y) du fond du réservoir, les coordonnées cartésiennes du point M sont définies par la relation suivante: x = (CH − y ) tg θ

(1)

Sur la Fig.1, on désigne par: - S1  S 2 : Contour délimitant la partie amont du barrage; - S3 : Contour définissant la surface libre de l’eau; - S4 : Contour définissant la limite de troncature du réservoir; - S5 : Contour définissant le fond du réservoir; - H : Hauteur du fluide dans le réservoir; - CH : Hauteur de la partie inclinée du barrage. Le système barrage-réservoir est bidimensionnel avec un réservoir de longueur infinie ayant un fond horizontal infiniment rigide. Le fluide dans le réservoir est considéré incompressible, non visqueux et irrotationnel, soumis à des excitations sismiques de courtes durées dans la direction horizontale. L’effet des ondes de surface éventuelles est considéré négligeable. 2.2. Formulation des équations gouvernantes

Sous les conditions précédentes, les pressions hydrodynamiques dans le réservoir sont gouvernées par l’équation des ondes de compression, qui prend, pour le cas particulier d’un fluide incompressible, la forme de l’équation de Laplace donnée par:


∇ 2 p( x, y )=

∂ 2 p( x, y ) ∂ 2 p( x, y ) + = 0 ∂x 2 ∂y 2

4

(2)

2.3. Conditions aux limites

Le réservoir est délimité par quatre frontières définies comme suit: 1. Sur la partie amont du barrage délimitée par le contour ( S1  S 2 ), on suppose que les particules fluides sont parfaitement solidaires avec les particules solides du barrage. Dès lors, sous l’effet d’un chargement sismique horizontal, le gradient des pressions hydrodynamiques dans la direction normale à la face amont du barrage et la force d’inertie développée dans la masse du liquide sont considérés en état d’équilibre, ce qui nous permet d’écrire : ∂p ( x, y ) = −ρ un ∂n S  S 1 2

(3)

où un représente le vecteur accélération des particules solides dans la direction normale à la face amont du barrage et ρ la masse volumique de l’eau. 2. Au niveau de la limite S3 , en absence de grandes fluctuations des pressions hydrodynamiques dans le réservoir, l’effet des ondes de surface est considéré négligeable, donc on prendra la pression atmosphérique p (x,H) = 0. 3. Le fluide étant considéré incompressible, la pression hydrodynamique générée par l’action sismique est considérée décroissante progressivement dans la direction x. Ainsi au droit de la limite S4 , on prendra p (y,∞) = 0. 4. Au niveau de la limite S5 (fond du réservoir), considérant que les particules fluides sont parfaitement solidaires avec les particules solides du barrage, sous chargement sismique horizontal, leurs vitesse dans la direction verticale est considérée nulle, ainsi le gradient des pressions hydrodynamiques correspondant est également nul, ce qui se traduit par: ∂p( x,0) =0 ∂n S 5

(4)

5


3. FORMULATION ANALYTIQUE DE LA RESULTANTE DES EFFORTS TRANCHANTS

Sous l’effet d’un chargement sismique horizontal, la pression hydro-sismique p pour le cas d’un barrage rigide à double fruit peut s’écrire sous la forme suivante: = p CS γ H C p

(5)

X g , γ est le poids volumique de l’eau, X g l’accélération horizontale du g sol, g l’accélération de la pesanteur et C p le coefficient des pressions

où CS =

hydrodynamiques. Le coefficient des pressions hydrodynamiques C p est, ici, approché par une série de fonctions de la variable réelle (y) [1] évoluant sur l’intervalle compact I = [0,H] donnée par la relation suivante: = C p ( y)

+∞

∑ Hi e −λ x cos λ i y , A

i

i =1

(2i − 1)π tel que λ i = , i = 1,2,...,+∞ 2H

(6)

Le calcul numérique des coefficients de pression C p ( y ) nécessite donc d’effectuer une troncature pour la série de fonctions donnée par l’équation (6). Pour une tolérance de calcul de l’ordre de 10-4, au-delà de 25 termes, on peut aisément démontrer par le biais du ‘critère de Cauchy uniforme’ la convergence uniforme de la série de fonctions C p ( y ) vers une fonction appelée « fonction limite de convergence uniforme ». Les coefficients inconnus Ai de la série de fonctions C p ( y ) sont obtenus après résolution du système d’équations linéaires donné par l’équation (7) ci-après à l’aide du programme de calcul numérique élaboré à cet effet.  F ji = { Ai }

{G j }

j, i {G j },=

1,2,..., +∞

(7)

Les éléments de la matrice symétrique  F ji  et ceux du vecteur colonne sont calculés suivant la formulation présentée dans la référence [1].

L’effort tranchant développé au droit du contour S1  S 2 , au-dessus d’une élévation quelconque (y) est défini comme suit: Ve ( y ) =

H

∫y

p ds

En substituant (5) et (6) dans (8), on obtient:

(8)


6

+∞ Ai −λ i x Ve ( y ) = CS γ H e cos λ i y ds H i =1 H

∫∑

(9)

y

Puisque C p ( y ) est une série de fonctions continues convergeant uniformément sur

∫

l’intervalle compact [0, H], l’interversion entre l’opérateur ( ) et l’opérateur (

∑ ) est donc permise. Ainsi, (9) prend la forme suivante: +∞

H

∑ ∫

Ai Ve ( y ) = CS γ H e −λ i x cos λ i y ds H i =1

(10)

y

En substituant x en fonction de y dans (10), on obtient après calcul intégral laborieux: Ve ( = y ) C S γ H 2C n ( y )

(11)

avec C n ( y) =

+∞

∑ λ i Hi 2 ( −e −λ (CH − y )tgθ sin(λ i y + θ) + sin(λ iCH + θ) − sin λ iCH + (−1) i+1 ) A

i

i =1

A partir de l’équation (11), on déduit les deux composantes de Ve ( y ) dans le repère (O, x, y) comme suit : Ve x ( y ) = Ve ( y ) cos θ ;

Ve y ( y ) = Ve ( y )sin θ

L’effort tranchant maximum est obtenu à la base du barrage, pour y=0: Ve max = Ve (0) = C S γ H 2C n max C n= max

+∞

∑ λ i Hi 2 ( −e −λ CH tgθ sin θ + sin(λ iCH + θ) − sin λ iCH + (−1) i+1 ) i =1

A

i

(12)

7


4. FORMULATION ANALYTIQUE DU MOMENT DE FLEXION

Le moment de flexion associé à Ve ( y ) au-dessus d’une élévation quelconque (y) est donné par: H

∫ V e ( y ) ds

M e ( y) =

(13)

y

Après calcul intégral (13) devient: M e(= y) C S γ H 2

+∞

−λ i (CH − y )tgθ

−e cos(λ i y + 2θ) + λ i (CH − y ) B +  ∑ λ i2 Hi 2  + cos( λ i CH + 2θ) + λ i H (1 − C )sin λ i H − cos λ i CH  A

(14)

i =1

avec B =

sin(θ + λ i CH ) − sin λ i CH + sin λ i H . cos θ

La valeur maximale du moment de flexion est obtenue à la base du barrage, pour y=0: = CS γ H 2 M e max

+∞

Ai  −e −λ i CH tgθ cos(2θ) + λ i CH B + cos(λ i CH + 2θ) +    λ 2 H 2  +λ i H (1 − C )sin λ i H − cos λ i CH  i =1 i

∑

A partir des relations générales précédentes (équations 11 et 14), on peut aisément en déduire les cas particuliers suivants : - Cas d’un barrage rigide à fruit simple (θ≠0, C=1) Ve ( y= ) CS γ H 2

+∞

∑ λ i Hi 2 ( −e −λ ( H − y )tgθ sin(λ i y + θ) + sin(λ i H + θ) ) A

i

i =1

M e(= y) C S γ H 2

+∞

 −e −λ i ( H − y )tgθ cos(λ i y + 2θ) + Ai   2 2   i =1 λ i H  +λ i ( H − y )sin λ i H + cos(λ i H + 2θ) 

∑

- Cas d’un barrage rigide à parement vertical (θ=0, C=0) Ve ( y )= C S γ H 2

+∞

∑ λ i Hi 2 ( − sin λ i y + (−1) i+1 ) A

i =1

M e ( y) = C S γ H 2

+∞

∑ λ i2 Hi 2 ( λ i ( H − y)sin λ i H − cos λ i y ) i =1

A


8

5. RESULTATS ET DISCUSSION

Maintenant que les efforts de réduction sont formulés analytiquement, dans ce qui suit on procédera à la validation des résultats numériques obtenus en comparant ces derniers, respectivement à ceux issus de l’application de la méthode exacte de Westergaard pour le cas d’un barrage rigide à parement vertical et à ceux de la méthode de Chwang lorsque le barrage présente un fruit simple. 5.1. Validation des résultats: Cas d’un barrage rigide à parement vertical (θ=0, C=0)

Dans ce qui suit les résultats numériques obtenus pour ce cas de figure sont comparés à ceux issus respectivement de l’application de la méthode exacte de Westergaard, et à ceux de Chwang. Il est à noter que pour ce cas particulier, la résolution du système d’équations linéaires (7) se réduit à un système diagonal d’équations. Le Tableau 1 ci-dessous présente les valeurs des erreurs relatives de p, Ve et M e , obtenues de l’application de la présente méthode et de celle de Westergaard pour diverses hauteurs du fluide de remplissage du réservoir. Table 1 Erreurs relatives de p, Ve et M e entre la présente méthode et la méthode exacte de Westergaard pour diverses hauteurs du fluide dans le réservoir.

p (%) Ve (%)

30.48 0.09 -0.03

M e (%)

0.21

Hauteur du fluide dans le réservoir H (m) 60.96 121.92 182.88 -0.92 -4.25 -8.19 -0.74 -3.28 -7.25 -0.69

-2.97

-6.63

243.84 -14.86 -13.35 -12.28

Pour des hauteurs faibles à moyennes (trois premières valeurs), on remarque que l’écart entre les résultats des deux méthodes est pratiquement négligeable. En revanche, lorsque la hauteur du barrage devient importante (deux dernières valeurs) l’effet de la compressibilité du fluide dans le réservoir se fait aussitôt ressentir, ce qui explique d’ailleurs, l’amplification de l’erreur relative entre les deux derniers résultats dans le tableau précédent. Nonobstant l’effet de la compressibilité de l’eau dans le réservoir, dans le Tableau 2, ci-après, on remarque que les résultats numériques obtenus par cette formulation analytique sont très proches des résultats issus de l’application de la méthode exacte de Chwang [3].

9


Table 2 Erreurs relatives sur p et Ve entre la présente méthode et la méthode exacte de Chwang.

Méthode analytique /Méthode Chwang p (%)

-0.02

Ve (%)

-0.05

Il est à noter que le signe (−) précédant certaines valeurs dans les deux tableaux précédents indique que les résultats de la présente méthode sont inferieurs respectivement à ceux de la méthode de Westergaard ainsi que ceux de la méthode exacte de Chwang. 5.2. Validation des résultats: Cas d’un barrage rigide à fruit simple ( θ≠ 0 , C=1 )

Pour ce deuxième cas de figure, la validation des résultats sera effectuée par analogie avec la méthode exacte de Chwang. L’équation (11) donne les coefficients C n de la résultante des forces hydrodynamiques dans la direction normale à l’interface barrage-réservoir. On en déduit les coefficients C x et C y comme suit: C y C n sin θ . C x C n cos θ ; = =

La Fig. 2 ci-après, montre en traits discontinus respectivement les courbes de C n , C x et C y obtenues par la présente méthode et en traits continus les courbes issues de l’application de la méthode exacte de Chwang. Un très bon accord est observé entre les deux résultats.

Fig. 2 – Comparaison entre les coefficients C n , C x

et C y respectivement donnés par la présente

méthode et la méthode exacte de Chwang (cas d’un barrage rigide à fruit simple).


10

5.3. Extension de l’étude pour le cas des barrages rigides à double fruit

La Fig. 3 ci-après, montre la distribution en élévation des efforts tranchants au-dessus d’une élévation quelconque y pour le cas d’un barrage rigide à double fruit. On remarque que plus la hauteur de la partie verticale augmente, plus la résultante des efforts tranchants devient importante. La partie inclinée du barrage permet de réduire sensiblement l’intensité des efforts tranchants. Le même raisonnement peut être adopté pour la distribution en élévation des moments de flexion (Fig. 4). Quel que soit l’inclinaison du parement amont, les valeurs maximales des efforts de réduction se produisent au fond du réservoir, contrairement au cas des pressions hydro-sismiques [10].

Fig. 3 – Distribution en élévation des efforts tranchants pour différents fruits composés.

Fig. 4 – Distribution en élévation du moment de flexion pour différents fruits composés.

11


6. CONCLUSIONS

Dans la présente recherche, des formules analytiques pour la détermination des efforts de réduction agissant au droit des barrages rigides à fruits irréguliers sous excitations sismiques sont développées. Les résultats numériques obtenus sont présentés pour différentes configurations géométriques des parements amont. A défaut de solutions exactes pour le cas des barrages rigides à double fruit, les résultats numériques obtenus sont comparés à ceux issus respectivement de l’application de la méthode exacte de Westergaard pour le cas particulier d’un barrage rigide à parement vertical, et à ceux de Chwang applicable uniquement pour le cas d’un fruit simple. Une très bonne concordance des résultats a été observée. Un programme de calcul numérique sous environnement MATLAB a été élaboré en vue de faciliter le calcul des efforts de réduction. L’intérêt pratique de cette approche réside dans le calcul des efforts sismiques en vue de l’analyse et le dimensionnement des barrages rigides à fruits composés. A partir des résultats obtenus, il est également possible de tirer les principales conclusions suivantes : 1- Pour le cas de barrages rigides présentant des hauteurs faibles à moyennes, l’utilisation de la présente méthode pour la détermination des efforts de réduction fournit des résultats très satisfaisants. En revanche, pour des hauteurs importantes, l’effet de la compressibilité du fluide dans le réservoir devrait être pris en considération. 2- Quelque soit la configuration géométrique de la partie amont des barrages à double fruit, les valeurs maximales des efforts de réduction se produisent au niveau de la base contrairement à celles des pressions hydro-sismiques [10]. 3- Le programme de calcul numérique élaboré dans le cadre de la présente étude devrait pouvoir être utilisé en vue de la détermination d’une configuration géométrique optimale des barrages à double fruit permettant une économie substantielle des quantités des matériaux de construction à utiliser. Received on August 9, 2013 NOTATIONS

Ai : Coefficients de la série de fonctions Cp(y). C : Rapport de la hauteur de la partie inclinée sur la hauteur du fluide dans le réservoir. CH : Hauteur de la partie inclinée du barrage. C p : Coefficient des pressions hydrodynamiques.

C S : Coefficient sismique. g : Accélération de la pesanteur (m/s2).


12

γ : Poids volumique de l’eau (N/m3). H : Hauteur du fluide dans le réservoir (m). i, j : Indices de boucle. M e : Moment fléchissant global de la force de pression p au-dessus d’une élévation y. n : La normale à l’interface barrage-réservoir. p : Pression hydro-sismique induite par un chargement sismique horizontal (N/m2). V e : Effort tranchant au-dessus d’une élévation quelconque y. x, y : Coordonnées rectangulaires. θ : Angle d’inclinaison du parement amont par rapport à la verticale. ρ : Masse volumique de l’eau (Kg/m3).

REFERENCES 1. AVILES, J., SANCHEZ-SESMA, F.J., Hydrodynamic pressures on dams with no vertical upstream face, Journal of Engineering Mechanics (ASCE), 112, 10, pp. 1054-1061, 1986. 2. CHOPRA, A.K., Hydrodynamic pressures on dams during earthquakes, Journal of Engineering Mechanics (ASCE), 93, pp. 205-223, 1967. 3. CHWANG, A.T., Hydrodynamic pressures on sloping dams during earthquakes, Part 2. Exact theory, Journal of Fluid Mechanics, 87, pp. 343-348, 1978. 4. FENVES, G., CHOPRA, A.K., Effects of reservoir bottom absorption and dam-water-foundation rock interaction on frequency response function for concrete gravity dams, Journal of Earthquake Engineering and Structural Dynamics, 13, 1, pp.13-32, 1985. 5. GOGOI, I., MAITY, D., A novel procedure for determination of hydrodynamic pressure along upstream face of dams due to earthquakes, Computers and Structures, 88, 9-10, pp. 539-548, 2010. 6. HATAMI, Kianoosh, Effect of reservoir bottom on earthquake response of concrete dams, Soil Dynamics and Earthquake Engineering, 16, 7-8, pp. 407-415, 1997. 7. LAPRESTE, J.T., Introduction à MATLAB, 2eme Edition avec MATLAB 7, Editions Ellipses, Paris, 2005. 8. LI, S.-M., LIANG, H., LI, A.-M., A semi-analytical solution for characteristics of dam-water reservoir system with absorptive reservoir bottom, Journal of Hydrodynamics, 20, 6, pp. 727734, 2008. 9. LIU, P. L.-F., Hydrodynamic pressures on rigid dams during earthquakes, Journal of Fluid Mechanics, 165, pp. 131-145, 1986. 10. TADJADIT, A., TILIOUINE, B., Pressions hydrodynamiques sur barrages rigides à fruits composés sous excitations sismiques, Colloque International sur la Réduction du Risque Sismique, Chlef, Algérie, 10-11 octobre 2012. 11.TILIOUINE, B., SEGHIR, A., A numerical model for time domain analysis of dams including fluid-structure interaction, CST 98 International Conference, Edinburgh, Scotland, August 18-20, 1998. 12. TILIOUINE, B., SEGHIR, A., Fluid-structure models for dynamic studies of dam-water systems, 11th European Conference on Earthquake Engineering, Paris, France, September 6-11, 1998. 13. WESTERGAARD, H. M., Water pressures on dams during earthquakes, Transactions, ASCE, 98, pp. 418-472, 1933.