M Guingand

J P De Vaujany

Y Icard

Fast Three-Dimensional Quasi-Static Analysis of Helical Gears Using the Finite Prism Method

This paper describes a model for analyzing external and internal helical gears. The finite prism method is used to study the elastic behavior of the structure. Contact deformations are also included in the model. Load sharing, pressure distribution, meshing stiffness and three-dimensional tooth fillet stresses are calculated at each instant. An experimental validation of the numerical model is also presented.

Introduction

Static-elastic analysis of the meshing of cylindrical gears is a challenging task as a number of modeling and numerical problems must be overcome. Two groups of problems are of primary interest. The first group is related to the computation of the contact pressure at the gear mesh. The pressures and the corresponding displacements within the contact area can be calculated by employing the theories of Hertz and Boussinesq. According to the actual geometry near the contact zone, either the results are given in the form of tabulated formulas, or a discretized contact area is used at the mesh where each grid square is subject to a certain contact pressure ͓1,2͔. Then, geometrical conditions of nonbody penetration allow formulation of a linear problem with consistency conditions ͓3͔. Once the pressures are known, the stress fields are calculated within the contact zone ͓4͔. When more than one pair of teeth transmit the load, the principle used for the calculation of load distribution is based on the analysis of generalized multicontacts of two bodies ͓3,5͔. Multitooth pair analyses of different studies subject show differences perhaps originating from the simplifying assumptions made ͓6-9͔.

The second group of problems of gear contact analysis are associated with bulk resistance. Stress calculations at the tooth root are made with the help of classical methods for continuous media. While many earlier studies used the material resistance method, general equations ͓10͔ are now used directly. Consequently, the finite element method ͑FEM͒ has been widely used in two-or three-dimensional ͑3D͒ modeling with a refined meshing ͓11-13͔. Therefore, it is necessary to develop computational tools that are refine enough to capture the major aspects of contact phenomena, but yet global enough to take into account the entire architecture of the gears.

The aim of this paper is to present a method for analyzing cylindrical gears in an efficient and precise manner. Parametrical analysis results using this method have been presented earlier previously ͓14,15͔, however, the fundamentals of this method and its experimental validation have not been published. The objective here is to formulate a elastic model of a cylindrical gear pair in such a way that reasonably accurate computation can be made within moderate calculation times, allowing the use of this formulation as a potential design tool. The finite prism method ͑FPM͒ is chosen to reduce the size of the numerical problem and accelerate the calculations. This method is able to consider the exact geometry of the teeth as well as the global axis 3D positions of the gears altogether. With this method, the bulk deformations can be modeled accurately as well as the stress distributions at the root of the teeth. In the following sections, FPM formulations are presented and are applied to gear pair problems.

The Finite Prism Method

Full analysis with the FEM is costly and even unnecessary for structures with a constant shape in one direction and simple boundary conditions. FPM can give results by using a much smaller number of input data ͓16͔. FPM can be considered as a special form of FEM, based on geometrical approximations concerning the shape of the displacements along a given direction. A standard FEM uses polynomial displacement functions in every direction while FPM uses simple polynomial displacements only along two directions where the meshing is carried out and serial functions for extrusion along the third direction. Thus, the 3D problem is approximated in a two-dimensional ͑2D͒ space, providing significant savings in computing time and memory requirements.

The general form of the displacement function is given as a product of polynomial functions and series. Therefore, for a linear prism where the z axis corresponds to the generating helical line, a displacement function is written in the following form:

Uϭ ͚ mϭ1 r N m ͑ x,y ͒F m ͑ z ͒ (1)
In this expression, the series are truncated to r terms. N m (x,y) is a polynomial expression including undetermined constants for the mth term. F m (z) is a series that satisfies the boundary conditions along the z-axis direction and also specifies the deflected shape of the structure in this direction.

The FPM used here is similar to the method proposed by Kantorovich ͓17͔, which is commonly used to reduce a partial differential equation to an ordinary differential equation. Cheung ͓18͔ developed the theory related to finite prisms. The displacement function U is often composed of two parts: the polynomial function N m (x,y), ruled by the shape of the transversal function and the series F m (z) satisfying the boundary conditions in the z direction. The transverse section of the prism is commonly modeled as an isoparametric eight-node element.

2.1 Choice of Displacement Functions. The global elastic line of a gear tooth along its face width is assumed to be the curve of a bending beam. Therefore, function F(z) is the solution of the spatial differential equation of a vibrating beam ͑see Fig. 1͒. The choice of this function is discussed by Olakorade ͓19͔. The differential equation is given in the form

d 4 F͑z ͒ dz 4 ϭ 4 b 4 z (2)
where b is the length of the prism in consideration, and is defined in Refs. ͓20͔ and ͓21͔. Therefore, the basic solution F(z) to Eq. ͑2͒ is given in the general form

F͑z ͒ϭC 1 sin ͩ b z ͪ ϩC 2 cos ͩ b z ͪ ϩC 3 sinh ͩ b z ͪ ϩC 4 cosh ͩ b z ͪ (3)
where coefficients C i are determined by applying the boundary conditions. They define the shape and amplitude of the displacements.

The length of the beam corresponds to the face width b of the gear teeth. In the particular case of ''free-free'' boundary condition, the prism is free at each end of the tooth

FЉ͑0 ͒ϭFٞ͑ 0 ͒ϭFЉ͑ b ͒ϭFٞ͑ b ͒ϭ0 (4)
where FЉ and Fٞ represent bending moment and transverse force terms, respectively and superscript ͑Ј͒ denotes differentiation with respect to z. These conditions express the fact that the bending moment and the transverse loads are zero at the free ends of the teeth. By substituting Eq. ͑4͒ into Eq. ͑3͒, four equations are obtained to allow the determination of C i of Eq. ͑3͒

F 1 ͑ z ͒ϭ1, F 2 ͑ z ͒ϭ1Ϫ 2 b z, . . . , (5)
F m ͑ z ͒ϭsin ͩ m b z ͪ ϩsinh ͩ m b z ͪ Ϫ␣ mͫ cos ͩ m b z ͪ ϩcosh ͩ m b z ͪͬ (6)
where

␣ m ϭ sin m Ϫsinh m cos m Ϫcosh m (7)
m Ϸ 1 2 (2mϪ3) (mϭ3,4,5, . . . ), and cos m cosh m ϭ1. In the other case of a tooth fixed at one end and free at the other end ͑case of a shaft pinion͒, the boundary conditions allow to write the conditions of a ''fixed-free'' prism

F͑0 ͒ϭFЈ͑ 0 ͒ϭFЉ͑ b ͒ϭFٞ͑ b ͒ϭ0 (8)
These equations show that the deflection and the rotation are zero at the fixed end and the bending moment and the transverse forces are zero at the free end. The corresponding serial functions for this case become

F m ͑ z ͒ϭsin ͩ m b z ͪ Ϫsinh ͩ m b z ͪ Ϫ␣ mͫ cos ͩ m b z ͪ Ϫcosh ͩ m b z ͪͬ (9)
where

␣ m ϭ sin m ϩsinh m cos m ϩcosh m (10) m Ϸ 1 2 (2mϪ1) (mϭ1,
F m u ϭF m , F m v ϭF m , F m w ϭF m Ј (12)
In the case of a spur gear pair, FPM can be used directly. For helical gears, a transformation of coordinates must be applied ͓5͔.

First, a curvilinear system (O,R,,Z) is used as shown in Fig. 2 such that

XϭR cos , Y ϭR sin , Rϭ͓X 2 ϩY 2 ͔ 1/2 , ϭArc tan ͫ Y X ͬ , Zϭz (13)
with RϾ0 and 0 рр f . Next, an oblique system (O,r,s,t) is applied in Fig. 2 as

Rϭr, ϭsϩt sin ␤, Zϭt cos ␤ (14)

Validation of FPM.

A pipe pressurized internally constitutes the first validation example. The pipe is fixed rigidly externally and an uniform pressure is applied inside as shown in Fig. 3. The serial function is defined along the axial direction. All the results obtained with a prism match well with the analytical predictions within an error nearly 4% on average. The larger differences are observed on the inner surface of the pipe.

In order to simulate the loading of the gear teeth, a clamped prismatic solid was considered as shown in Fig. 4 and the FPM results were compared with those obtained by using a classical 3D FEM. The main objective was to validate the results with different boundary conditions ͑fixed-free and free-free͒ and also to investigate the appropriate number of terms that must be considered in the series. The block is fixed rigidly in y direction for all z, the nodes are either free on the two small sides ͑case free-free͒, or fixed at boundary at zϭ0 ͑case of fixed-free͒. The stresses are analyzed on the dotted line aaЈ ͑see Fig. 4͒ close to the rigidly fixed side. The displacements are observed on the bbЈ line and the load is applied on top of the tooth. Figures 5 and6 show the results for the free-free boundary condition along the length of the beam. The displacements ͑Fig. 5͒ and the values of the principal stress ͑Fig. 6͒ are given for different numbers of terms in the series (T1 is one-term, T2 is two-term and soon͒. The results are compared to those obtained with the 3D finite element method.

For the meshing, hexahedron elements with 20 nodes are used in a commercially available finite elements package ͑two elements in X, four elements in Y, and ten elements in Z. Different meshing conditions are also considered to analyze the effect of the meshing smoothness on the results of both between FEM and FPM. The results show that, in this example case, eight terms are necessary to achieve acceptable accuracy in displacement predictions and five terms for stress computations.

In summary, the FPM uses a 2D mesh on a transversal section of an object, and an implicit approximation of the displacements is performed along the third dimension. This method gives satisfying results, compared to those obtained by other conventional methods. In addition, the calculation time is much shorter ͑80 times shorter according to Olakorede ͓5͔͒, and the memory space required is much smaller ͑nearly six times less͒ compared to conventional finite elements. Therefore, this method can be useful for solving gear design problems. ͑1͒ After determining the active profiles of the teeth, tooth contact analysis ͑TCA͒ simulation without load is performed for several meshing positions.

Application of FPM to Gear Analysis

͑2͒ This preliminary TCA allows the definition of the initial contact points and the unloaded transmission error caused, for example, by gear eccentricity, modified profile or undesirable profiles errors.

͑3͒ A 2D mesh of the gear is then produced to introduce the FPM method. Influence coefficients are used to calculate load sharing. A unit force is applied on the tooth contact surface and produces a displacement. Then, by considering unit forces at different points, a set of bulk influence coefficients can be deduced. A model of the contact calculates the surface effects between the teeth. At this stage of the process, the contact pressures and the static loaded transmission error are obtained. Finally, once the load sharing is known, the FPM is used again to calculate the 3D stresses at the root of the teeth at the last step of the process. inside the contact zone. The fixed-point algorithm is used to solve the system of equations. Coefficients C ik are the influence coefficients obtained with a uniform pressure P k applied on a small surface S k around point k of the contact zone of body i.

Load

Two terms define C ik , C ik ϭC ik f ϩC ik s . Term with superscripts f and s represents the bulk deformations and the surface deformations, respectively. C ik f is calculated using the FPM method and C ik s is found by using the wellknown Boussinesq theory ͓22͔. In order to implement these calculations, a local square mesh is produced in the tangent plane of the contact zone as shown in Fig. 9. The influence coefficients are obtained by applying unit loads at each node that are likely to be loaded. Two types of result are obtained: ͑i͒ load sharing and ͑ii͒ contact pressure distribution. Load sharing refers to the sharing of the torque among meshing tooth pairs. The contact pressure distribution is found along the face width for the different teeth in contact. 3.3 Model Results. In this section, results for an example helical gear pair are presented. The design parameters of the helical gear pair as listed in Table 1. The input torque is 2000 N m and the involute gear tooth profiles are not modified. Figure 10 presents potential contact lines on the teeth of the driving pinion for three different kinematic positions. For position 1, only teeth 2 and 3 are in contact. For an intermediate position ͑position 5͒, three pairs of teeth are in contact ͑teeth 2, 3, and 4͒. For tooth 3, the contact line on tooth 3 extent is along the entire tooth width. Tooth 2 experiences only a contact at the tip of tooth. At the last kinematic position, tooth 2 is not in contact while tooth number 5 enters the mesh.

The load sharing characteristics are presented in Fig. 11 as a result of computation of the elastic behavior. Tooth 1 is not loaded for this range of kinematic positions. Tooth 4 enters in the contact at position 2 and tooth 5 at position 8 while tooth 2 leaves the contact at position 6. The contact pressure distribution along the flanks of teeth 2 and 3 is presented in Fig. 12 for the kinematic position 1. Peaks of pressure on the edge of the teeth can be noted. Similarly, Fig. 13 presents the field of contact pressure on the square grid defined around the contact zone for several kinematic positions, from the beginning of a tooth meshing, to the exit of the tooth from the contact zone. The transmission error and the gear mesh stiffness are also predicted. An example of stiffness is shown in Fig. 14 as a function of the angular position of the driving pinion. The stiffness is calculated with the ratio of the input load to the transmission error ͑␣͒. The stresses in the tooth roots are computed by using the distribution of the load. Figure 15 presents the maximal principal stresses in four consecutive roots at kinematic position 5. Roots 2 and 3, located between loaded teeth experience compressive stresses.

Experimental Validation

Any numerical tool must be validated experimentally before it can be used in confidence as a design and analysis tool. This section aims at validating the gear analysis model presented in earlier sections experimentally.

Tests were performed on an existing back-to-back type fatigue test machine shown in Fig. 16. Table 2 specifies the parameters of the helical gears used in the experimental test. The gear manufacturing quality is ISO4 and tooth profiles are not modified. The schematic layout of the machine is shown in Fig. 17 were mounted on high-capacity tapered roller bearings. A rotating hydraulic torque device applies different torque values on the inner test loop and the system is driven by an electric motor.

During the tests, the angular rotation of the pinions is kept very low ͑5 rpm͒ so that a comparison with the quasi-static model is possible with no dynamic effects. The gears were instrumented with strain gages. A total of six gauges were mounted on each gear, as shown in Fig. 18 at a 120 deg tooth fillet position and two locations along the tooth fillet width. The gauges are located near the peak bending stress location. An optical encoder was also used to mark the angular position as gears rotate. The tests were performed for one rotation of the larger gears that corresponded to nearly five rotations of the pinion, so that the strain-stress variations could be monitored for each strain gage. When the strain gauged tooth fillet entered the meshing zone, the signal varied from positive to negative values ͑see Fig. 20͒. The results given for the circular positions 1, 1Ј or 1Љ are also plotted in Fig. 20.

The tested gear was modeled numerically. Figure 19 shows the FPM of the test gears and the boundary conditions used. The rim thickness is greater than 3.5 times the module, a five-tooth segment is sufficient to perform a reasonably good numerical simulation. The input torque value was set at 600 N m. In Fig. 20, example predictions and the corresponding stresses during meshing are compared. The measurements for the three gauges situated at an angular distance of 120 deg agree very well with the simulation results ͑solid line in Fig. 20͒. Differences between the numerical and measured results can be explained by secondary effects such as pitch error, shaft misalignments, and helix angle errors ͓23͔.

Conclusions

The present paper describes a numerical tool that aids the design of cylindrical gears. The elastic behavior of gear contact was studied with the finite prism method. This method has been shown to correlate to the results of conventional FEM well. It was shown that the numbers of terms in the series functions influence of the accuracy of the predictions. The integration of this method to calculate the gear bending coefficients and the application of Boussinesq theory allows to computation of the load sharing and the tooth root stresses. The FPM gives results in three dimensions with significantly less computational demand. An experimental test results were also shown to match the predictions of this model well.

  Fig. 1 Gear prism element

3. 1

 1 Description of the Process. A realistic geometry of the gear distinguishes three levels of definition: the architecture ͑web and rim͒, the position of the gear ͑misalignment, eccentricity, instantaneous meshing position͒, and the teeth ͑profile, crowning, surface waviness, pitch, and distortion errors͒. Obviously, for each meshing position, the load is shared among the meshing teeth, contact pressure distribution is established and a stress field spreads inside the considered elements. The complete process, integrated in the design software called PRINCE ͑French acronym for: Program Interactif de Calcul d'Engrenage͒, can be summed up as follows ͑Fig. 7͒.

Fig

  Fig. 4 A prismatic solid

  Analysis. The determination of tooth load distribution is a problem of contacts. Once potential contact zones are established, the load sharing in these areas must satisfy the conditions of displacement compatibility in every point k in Fig. 8. Inside the contact zone W 1k ϩW 2k ϩ⑀ k Ϫ␣ϭ0 and P k у0, (15) and outside the contact zone W 1k ϩW 2k ϩ⑀ k Ϫ␣у0 and P k ϭ0 (16) W 1k and W 2k are the normal displacements of bodies 1 and 2 at point k, respectively, P k is the contact pressure at point k, ⑀ k is the normal initial deviation at point k ͑computed in the unloaded TCA͒, and ␣ is the global adjustment between the gear bodies ͑transmission error expressed on the line of action͒. The normal displacement of body i at point k can be written as W ik ϭ ͚

Fig

  Fig. 8 Contact zone

  . The shafts

Fig. 13

 13 Fig. 13 Contact pressure during meshing on one helical tooth

Fig. 18

 18 Fig. 18 Location of the gauges

Table 1 Parameters of the example helical gear pair

 1

	Pinion	Gear

Table 2 Geometrical characteristics of the gears

 2

		Pinion	Gear
		͑driving and driven͒	͑driving and driven͒
	Number of teeth	159	32
	Working module ͑mm͒ Width ͑mm͒ Pitch diameter ͑mm͒ Helix angle ͑deg͒	3.56 86 566.039 19.967	3.56 80 113.92 19.967
	Helix hand	Left	Right
	Working pressure angle ͑deg͒ Tip diameter ͑mm͒ Root diameter ͑mm͒	25.347 572.1 556.1	25.347 121.2 105.1

Acknowledgments

The initial research program was financed by the EUROCOP-TER FRANCE Company mainly during the Ph.D of A. A. Olakorede. Additional works were performed in the framework of the Ph.D of H. C. Kim and J. P. de Vaujany. The authors thank D. Remond and J. Mahdoudh for their help during the experiments.


F Zara

T Redarce

P Poignet

Theme F "Medical Robotics for training and guidance": results and future works Thème F « apprentissage et assistance aux gestes médico-chirurgicaux » : bilan et perspectives

This paper presents the projects of the Theme F "Medical Robotics for training and guidance" inside the GdR STIC Santé. Three scientific meeting days have been organized during the period 2011-2012. They were devoted to physical simulators of behavior for gesture learning, command of hand prostheses by myoelectric signals or brain activity and the manipulation of objects by the artificial hand, and the last to the use of robots for medical gestures. The next event, scheduled for early 2013, will focus on the evaluation of gesture and especially "Evaluation of gesture -To do what ?".

Aims

The Theme F "Medical Robotics for training and guidance" inside the GdR STIC Santé aims to promote interaction between scientific and medical communities and to improve the visibility of the community. Its secondary goal is to develop mutual response of the ANR (French National Research Agency) to create new collaborations between these two communities. This research group focuses on three main themes: (1) modeling and visualization of organs; (2) learning tools; (3) interaction human / machine. Indeed, the final goal is to develop pertinent simulators and platforms demonstration for medical training and scientific research.

The scientific challenges concern (1) the definition of new biomechanical models, parallel algorithms and data structures to obtain real-time simulation of organs [START_REF] Buttin | A Biomechanical Model of the Female Reproductive System and the Fetus for the realization of a Childbirth Virtual Simulator[END_REF][START_REF] Baudet | Integrating Tensile Parameters in Hexahedral Mass-Spring System for Simulation[END_REF][START_REF] Saadé | A preliminary study for a biomechanical model of the respiratory system[END_REF][START_REF] Promayon | A 3D Discrete Model of the Diaphragm and Human Trunk[END_REF][START_REF] Nesme | Physically Realistic Interactive Simulation for Biological Soft Tissues[END_REF]; (2) the elaboration of pertinent scenarios of training and their evaluation [START_REF] Vadcard | Developpement de l'expérience en chirurgie orthopédique percutanée[END_REF][START_REF] Chieu | Student Modeling in Orthopedic Surgery Training: Exploiting Symbiosis between Temporal Bayesian Networks and Fine-grained Didactic Analysis[END_REF][START_REF] Tonetti | Assessment of a percutaneous iliosacral screw insertion simulator[END_REF][START_REF] Sclaverano | BiopSym: A Simulator for Enhanced Learning of Ultrasound-Guided Prostate Biopsy[END_REF][START_REF] Luengo | Novel Technology for Learning in Medicine[END_REF]; (3) and the proposition of new physical devices for the user to interact with the numerical simulations [START_REF] Dupuis | Does forceps training on a birth simulator allow obstetricians to improve forceps blade placement[END_REF][START_REF] Moreau | Simulation of an instrumental childbirth for the training of the forceps extraction: control algorithm and evaluation[END_REF][START_REF] Moreau | Assessment of forceps use in obstetrics during a simulated childbirth[END_REF][START_REF] Moreau | A method to analyze skill based on the curvatures of the position and the orientation: application to the evaluation of obstetric gestures[END_REF][START_REF] Moreau | Design of a new instrumented forceps: Application to a safe obstetrical forceps blades placement[END_REF][START_REF] Dupuis | A new obstetric forceps for the training of junior doctors: a comparison of the spatial dispersion of forceps blade trajectories between junior and senior obstetricians[END_REF][START_REF] Goulette | Computer Assisted Dental Implantology : A New Method and a Clinical Validation[END_REF][START_REF] Chendeb | Design and Realization of a Force Feedback System Dedicated to Virtual Surgery[END_REF][START_REF] Goulette | A Framework for Fast Computation of Hyper-Elastic Materials Deformations in Real-Time Simulation of Surgery[END_REF]. *Author to whom all correspondence should be addressed e-mail : florence.zara@univ-lyon1.fr

Results for the per iod 2011-2012

During the last period, three scientific meeting days were organized, each one being focused on a specific research topic, and mobilizing the different communities. Note that, some of these days are organized with research groups of the GdR Robotique.

More details about these meeting days are presented on the research group Web page: http://stic-sante.org/gdrpresentation/45-themef.

November 7, 2011 -The first meeting was devoted on the physical simulation of the behavior. Indeed, physical simulation is essential and coupling numerical simulation / physical simulation remains a scientific lock for learning and assistance actions of medical gestures. Moreover, following the meeting days organized in 2008 and 2009 by the research group, around the modeling and simulation of organs, which have dealt with aspects of numerical simulation and haptic interactions, we wanted to now focus our attention on the simulation of behavior.

November 25, 2011 -The second scientific meeting day was devoted to review the research activities in the field of prostheses active hand in France including commercial products available and developments in laboratories. Moreover, we are interested in the control of prostheses by myoelectric signals, in the developments of brain machine interfaces, and in the manipulation of objects by artificial hands, as in "full hand" as your fingertips.

May 30, 2012 -The third meeting was an open meeting day with several presentations: design of compliant mechanisms for surgical robotics; co-operation of an ultrasound probe for assisted biopsy of the prostate; handling magnetic gastrointestinal endoscopy; automatic recognition of high level and low level surgical tasks in the operating room from video images; robotics use in confocal endomicroscopy.

Following the different meeting days of the research group, a consortium has been created with LIRIS (Université Lyon 1), laboratory TIMC-IMAG (Université Grenoble 1), CAOR (Ecole des Mines ParisTech), laboratory Ampère (INSA de Lyon), LSE (Université Grenoble 2), Professor Olivier Dupuis (MD, HCL Lyon-Sud) and the society Didhaptic of Laval. This consortium has presented the SAGA project (Simulator for learning medical gestures of childbirth) to the ANR. It has been accepted (2013-2016).

Futur e plans

The next meeting will take place in March 2013 focused on the evaluation of the gesture. Moreover, we plan to organize in 2013 a school about the elaboration of medical simulators for training. This school will present issues of medical simulation (modeling of organs from medical data, biomechanical models, constitutive laws of organs, modeling interactions between organs and medical devices, etc.), and the coupling with suitable haptic interfaces. It will be open to students and researchers in this community.


J Bender

J Dequidt

C Duriez

G Zachmann

F Caillaud

X Faure

F Zara

F Jaillet

J-M Moreau

Workshop on Virtual Reality Interaction and Physical Simulation VRIPHYS (2013)

Keywords: Descriptors (according to ACM CCS): I, 3, 5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically based modeling I, 6, 8 [Simulation And Modeling]: Types of Simulation-Parallel

9:8)1(%)1" 7()"1; 7()%77 %;%)<* $%&'%()18"9:8)1(%)1" %4<% ")%" )"(% =%76 >:?

: Overview of our approach to dynamically adapt the simulation according to criteria.

Introduction

The accurate and interactive simulation of deformable objects has become a challenge in computer graphics [DCA99, SDR * 05]. As an alternative to the classical approach of the FEM, the Tensor-Mass solves the mechanical equations by directly computing forces applied on each node, instead of forces applied on each element. This approach offers a good compromise between accuracy, computation time and control.

A symbolic computation approach for this model [FZJM12] enables to multiply the possible combinations between geometrical elements (quads, triangles, tetrahedra, hexahedra with different integration nodes) and constitutive laws (from linear and non linear elasticity). This new possibility gives wide opportunities at low development cost, to mix several mechanical behaviors with different kinds of elements inside a same object.

In this paper, we explore the opportunity to improve the simulation in precision or time, by adaptively modifying the model according to criteria (see Fig. 1).

Background on the Tensor-Mass Model

Considering an object's element E with n nodes, the contribution of the forces applied on each node Pi for i ∈ [0, n -1] of this element is defined by: FE(Pi) = -∂WE(Pi) ∂Ui with WE the total strain energy of the element in 2D defined by:

WE = E W law (x, y) dx dy = ng -1 i=0 wi W law (ξi, ηi) det(J)
with J the transformation between spatial (x, y) and local (ξ, η) coordinates, ng the number of sampling points (ξi, ηi) of weights wi and W law the strain energy depending on the law (Hooke's or Saint Venant-Kirchhoff's law).

The total force F i applied on node Pi is then computed by summing the contributions from all the neighboring elements.

At time t, the Newton equation governing the object's movement involves the following formulation for the node Pi of mass mi and displacement U i :

mi Üi (t) = Fi(t).
This equation is linked to the differential equation system by assuming the sparse matrix of the masses to be diagonal, with the distribution of the object's mass on each node (mass-lumping).

Adaptive Tensor-Mass simulation

As stated previously, the possibilities offered by [FZJM12] allows us to modify an object or a scene during the simulation. The necessary actions are part of an adaptive strategy, and may concern the constitutive law (Hooke, Saint Venant-Kirchhoff, etc.), the geometry (quad, triangle, etc.), the refinement (addition or suppression of nodes), or the shape functions used for each element of the object (Fig. 1). Algorithm 1 summarizes our approach.

Algorithm 1 Adaptation at each step of the simulation according to the evaluation of the criteria. endfor 14: endfor Geometrical criteria The Hookean model is known to reach its limit for large displacements (detected by the edge length ratio RL) or for rotations (where unwanted area/volume increase is detected by RA). When this happens, a good solution is to switch to a non-linear law, leading to a more accurate simulation while controlling the computation time (see Fig. 2) † .

Other criteria like Aspect Ratio or Harmonic Mean of the length of the edges can be useful to consider the quality of the elements. This time, a geometric action like refinement is advocated (as a change in the constitutive law should have limited action). † The Open Source Framework SOFA was used for the implementation [ACF * 07].

Physical criteria. Besides classical geometric estimation errors, we can also consider a physical criterion noted RS linked to the stress applied on each element E of area A(t) at time t, with:

RS = 1 A(t) E FE(x, y)dx dy.
This physical criterion has to be discretized. But as the forces at the Gaussian point gi are not known, we also need to interpolate it according to the forces applied on the vertices. We obtain:

RS = 1 A(t) ng -1 i=0 wi n-1 j=0 Λj(gi) FE(Pj)
with FE(Pj) the force applied on the vertex Pj and wi the weighting function applied on the Gaussian point gi. Note that most of the terms have already be computed during simulation. Fig. 3 shows our results.

GPU implementation

A GPU implementation of the Tensor-Mass model has been proposed in [FZJM12]. But in adaptive simulation on the GPU, elements are handled by groups with the same association geometry/constitutive law. When an element is promoted to another group, this implies to keep up-to-date the list of elements' indices; and to load the new list for each group as well as all information linked to the considered elements. The GPU implementation involves: 1. low memory (50 times less than Abaqus TM ). Thus, for each group, we can set an array of the number of elements, meaning that there is no need to manage resizing at each adaptation step; 2. reduced reading and writing times to transfer data from RAM to GPU memory (14 % of the total time, in the worst case).

  1: {lc ← list criteria} 2: for each element e do 3: p ← position of element's nodes 4: u ← displacement of element's nodes

FFigure 2 :

 2 Figure 2: Simulation of a cantilever beam. The green-red range stands for max (RL,RA) ∈ [0, 1]; and blue when the quad element has shifted from Hooke to Saint Venant-Kirchhoff constitutive law.

Figure 3 :

 3 Figure 3: Gravitational force applied to a beam fixed on its top. Successive refinements (a-h) of the mesh according to the RS criterion. At each step, triangles overpassing 90 % of max(RS) are subdivided.

c The Eurographics Association 2013.


X Faure

F Zara

F Jaillet

J-M Moreau

Workshop on Virtual Reality Interaction and Physical Simulation VRIPHYS

Keywords: (according to ACM CCS): I, 3, 5 [Computer Graphics]: Computational Geometry and Object Modeling-Physically based modeling I, 3, 7 [Computer Graphics]: Three-Dimensional Graphics and Realism-Animation and Virtual Reality I, 6, 8 [Simulation And Modeling]: Types of Simulation-Parallel

The realistic and interactive simulation of deformable objects has become a challenge in Computer Graphics. In this paper, we propose a GPU implementation of the resolution of the mechanical equations, using a semi-implicit as well as an implicit integration scheme. At the contrary of the classical FEM approach, forces are directly computed at each node of the discretized objects, using the evaluation of the strain energy density of the elements. This approach allows to mix several mechanical behaviors in the same object. Results show a notable speedup of 30, especially in the case of complex scenes. Running times shows that this efficient implementation may contribute to make this model more popular for soft bodies simulations.

Introduction

Following the increasing demand of realism in Computer Graphics, physically-based simulation has become a very active research field over the last decade. This is particularly apparent in cloth and hair animation, medical simulation, interactive entertainment, and more generally in all Virtual Reality applications where animation, interaction or alteration of deformable objects is required in interactive time.

Several methods have been published to model the physical behavior of deformable objects [NMK ⇤ 06]. Among them, the Finite Element Method (FEM) allows the resolution of the differential equations systems governing the movements of objects [START_REF] Works | Finite element modeling for stress analysis[END_REF][START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. This method is based on the discretization of each object into elements (hexahedra or tetrahedra, for example) and provides the resolution of a linear system of the form KU = F, with K the stiffness matrix of the considered object, U the displacement of all nodes of the object, and F the external forces and boundary conditions applied on the object. Moreover, K may be dependent of U, judging from the mechanical behavior of the object. As the physical properties of objects are directly integrated in this mechanical formulation, this produces realistic, physically-based simulations of deformable objects.

As an alternative to the classic approach used in the FEM, the Tensor-Mass (TM) approach was introduced in Computer Graphics by Delingette, Cotin and Picinbono [DCA99, CDA00, PDA00, Pic03], and extended by Schwartz [SDR ⇤ 05]. This approach consists in another way to solve the mechanical equations of the objects, by directly computing the forces applied on each node considering the evaluation of the strain energy density.

To the difference of the classical FEM, the forces between elements do not need to be expressly computed, as only the forces applied on nodes are considered. This approach allows to consider objects discretized on elements with different geometrical models or mechanical properties. The main difficulty remains in the formulation of the equations involved by the forces computation, more specially using an implicit integration scheme like the Euler implicit one, which involves the computation of the differentials of the forces. In our work, this problem is addressed using formal computation.

Several formulations of the Tensor-Mass model have been presented in Computer Graphics to account for various mechanical behaviors: the linear Hookean model [START_REF] Cotin | A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation[END_REF]; the non-linear geometrical model based on Saint Venantc The Eurographics Association 2012.

Kirchhoff's elasticity model [PDA00], or anisotropic material [START_REF] Picinbono | Non-linear anisotropic elasticity for realtime surgery simulation[END_REF]. Moreover, Schwartz [SDR ⇤ 05] presented another extension for non-linear visco-elastic deformations, with some pre-computing of the various tensors to accelerate the process.

Concerning parallel physically-based animation on the GPU, several papers were presented: for non-linear FEM soft tissue modeling based on CUDA [CTAO08], FEM cloth simulation [RNSS ⇤ 06], implicit FEM solver for deformation simulation [ACF11], or Mass-Springs Systems [MHS05, GW05]. With respect to the Tensor-Mass model, Mosegaard [SM06] addressed its parallelization for the linear case with an explicit integration scheme, using texture memory.

But, as far as we know, the topic of the parallel implementation on the GPU of the Tensor-Mass model has not yet been addressed neither for the non-linear elasticity model, nor using an implicit integration scheme. An approach based on the texture memory as in [SM06] would be possible, but would not permit to consider topological changes in the simulated object's mesh.

In this work, we consider both linear and non-linear elastic models, with two of the most employed formulations in interactive simulations, namely Hookean and Saint Venant-Kirchhoff's models. Moreover, we present the main steps of the resolution of the mechanical equations using the Tensor-Mass approach, employing a semi-implicit as well as an implicit integration scheme to ensure unconditional stability for any time step. The topic of its parallelization on the GPU is also addressed. Obtained results exhibit notable speedup, enabling the interactive simulation of soft-bodies with linear or non-linear mechanical behaviors. Running times shows that this implementation will be valuable for relatively large scenes, contributing to make this model more popular in the domain of deformable object simulation.

Simulation of object deformations

Continuum mechanics deals with a continuous way to describe an object moving or being deformed under the action of stress. The deformation F may be formulated according to the displacement U(X) of a point X of the object by F(X) = X +U(X).

Different elasticity models exist depending on the expected mechanical behavior of the material. The straintensor e is deduced from the deformation gradient rF = I + rU, a quantity depending on rU. Therefore, e may be directly expressed as a function of rU. Moreover, the deformation energy W , involved by the deformation, naturally depends on the material's characteristics. Its derivative gives the corresponding opposite force causing the deformation.

Finally, the movement of the deformable object may be expressed by the following differential equations system:

M Ü + D U + KU = F, (1)
with U the displacement of the object and M, D, K respectively the mass, damping and stiffness matrices of the simulated object.

As for the classical FE approach, the TM formulation is based on the domain's discretization into several elements, and naturally, some steps are similar. However, its specificity is the discrete and local resolution of the mechanical equations. Details of the mains steps involved for each node of each element are now presented.

Discretization of the displacement in an element

The displacement of a point X(x, y, z) inside an element E of the object is defined by:

U E (X) ' n 1 Â j=0 L j (X) U j , (2)
with n the number of nodes P j = (P jx , P jy , P jz ) defining the element, U j the displacement of each node P j from its initial position, and L j (X) some interpolation functions. These functions are defined according to the type of element used for the discretization, with:

n 1 Â j=0 L j (X) = 1. (3)
One of the simplest and most used subdivision of the domain consists in P1 tetrahedral elements (with n = 4 interpolating nodes placed at each vertex). Indeed, this kind of element corresponds to a linear interpolation, that is valid for both considered cases in this paper, namely Hookean and Saint Venant-Kirchhoff's materials, as pointed out in [START_REF] Picinbono | Non-linear anisotropic elasticity for realtime surgery simulation[END_REF].

The linear interpolation functions are defined using the barycentric coordinates of X inside the element, leading, for a tetrahedron, to:

L j (X) = a j • X + b j (4) = a jx x + a jy y + a jz z + b j (5)
for j = 0..3 with (considering O the origin and OP j the position of P j from this origin):

⇢ a j = ( 1) j (OP j+2 OP j+1 ) ⇥ (OP j+3 OP j+1 ) b j = ( 1) j OP j+1 • (OP j+2 ⇥ OP j+3 ). (6)

Strain energy for an element

Let us define the strain-tensor and the energy of deformation according to the mechanical behavior of the object.

Hooke's law. For small deformations (or strains) -usually less than 10% of the size of the object -the elasticity is considered as linear with a linear relationship between stress and strain. To model such behaviors, only the linear part of the Green-Saint Venant strain-tensor is considered:

e l (X) = 1 2 ⇣ rU T (X) + rU(X) ⌘ , (7)
The energy density of deformation, measuring the strain energy per unit undeformed volume on an infinitesimal domain around the material point X, is then defined by

W l (X) = l 2 (tr e l (X)) 2 + µ tr e l (X) 2 , (9)
where l and µ are the Lamé coefficients characterizing material stiffness.

Saint-Venant Kirchhoff's constitutive law. For displacements larger than 10% of the object's size, the non-linear elasticity behavior must be considered. The simplest hyperelastic model is a direct extension of the linear elastic material. It is based on the Green-Saint Venant strain-tensor e nl , and its associated strain energy density W nl with:

e nl (X) = 1 2 ⇣ rU T (X) + rU(X) + rU T (X) rU(X) ⌘ W nl (X) = l 2 (tr e nl (X)) 2 + µ tr e nl (X) 2 . ( 10
)
Generalization. As rU and rU T are constant inside an element when considering P1 tetrahedral elements, the straintensor and the deformation energy density (noted e law and W law according to the chosen mechanical law -Hooke's or Saint Venant-Kirchhoff's) are also constant inside the tetrahedron (i.e., do not depend on X). Hence, the total strain energy W E of a P1 tetrahedron follows:

W E = Z E W law (X) dX = W law Z E dX = W law Vol E , (11)
with Vol E the initial volume of the considered element.

Computation of the forces and their differentials

Considering a given element E of the discretized object, the forces applied on any node P i for i 2 [0, n 1] of this element, is defined by:

F E (P i ) = ∂W E (P i ) ∂U i , (12)
with W E (P i ) the energy density of deformation of the considered element, evaluated at node P i .

Whereas the information of the forces is sufficient for an explicit or semi-implicit integration scheme (enabling to compute velocities and displacements according to accelerations and velocities, respectively), the computation of their differentials is mandatory when an implicit integration scheme is used [TPBF87, BW98]. Considering an element E, the differentials of the force on node P i is given by:

dF E (P i ) = h dF 0 E (P i ) . . . dF n 1 E (P i ) i , (13)
where

dF j E (P i ) = h ∂FE (Pi) ∂Uj i 33 for j 2 [0, n 1].
Finally, if we consider the whole object involving m nodes, the force applied on P I with I 2 [0, m 1] is computed by summing all the contributions from the neighboring elements, with:

F(P I ) = F I = Â E2N I F E (P i ) = Â E2N I ∂W E (P i ) ∂U i , (14)
where N I is the set of elements containing node P I , and with i 2 [0, n 1] the local node in E corresponding of the global node I of the object (in the global vertex indexation).

Considering the whole object involving m nodes, the differentials of the force on P I with I 2 [0, m 1] follows:

dF(P I ) = dF I = h dF 0 I . . . dF m 1 I i , (15)
with, for J 2 [0, m 1]

dF J I = Â E2N J I dF j E (P i ) = Â E2N J I  ∂F E (P i ) ∂U j 33 ,
where N J I is the set of elements containing global nodes P I and P J , and with i (resp. j) 2 [0, n 1] the local node in E corresponding of global node I (resp. J ) of the object.

Equation of movement for each node

Once the forces applied on each node have been computed, Newton's equation governing the movement of the object may be used. At time t, considering the node P i of mass m i , we have:

m i Üi (t) = F i (t). ( 16
)
This equation is linked to the differential equations system (1) by assuming the sparse matrix M to be diagonal, and distributing the object's mass (computed according to its density r) on each node of its discretization. This simplification, called mass-lumping, yields to the formulation of equation ( 16) for each node, which enables the computation of the acceleration according to the applied forces.

The next section will explain in particular our parallel implementation of Euler's semi-implicit and implicit integration schemes used to obtain velocity (from acceleration) and displacement (from velocity) of the nodes throughout simulation.

Parallel integration in SOFA

The Open Source Framework SOFA was used for our parallel implementation of the Tensor-Mass model using the OpenCL language [START_REF] Stone | Opencl: A parallel programming standard for heterogeneous computing systems[END_REF]. This framework, written in C++ and using the XML script language, enables comparisons between models and methods implemented by research groups in medical simulation (http://www.sofaframework.org) [ACF ⇤ 07].

c The Eurographics Association 2012.

The objects simulated in SOFA are defined by their topology, geometry and visual model (by loading OBJ files, for example), their mechanical state (position, velocity, acceleration and force stored as vectors), their mechanical behavior (constitutive laws), as well as the collision model. Then, the simulation loop is executed, involving the following main steps for each node of the discretized object: (1) computation of the forces, (2) computation of the acceleration according to Newton's second law of motion, (3) integration of the acceleration to obtain velocity, (4) integration of the velocity to obtain displacement.

To integrate the Tensor-Mass model in SOFA, only two functions, namely addForce() and addDforce(), should be implemented, enabling the computation of the forces and their differentials (only required for Euler's implicit integration scheme) for each node. But, as the initial GPU implementation of SOFA is based on CUDA [CTAO08], the various steps of the simulation loop have also been parallelized using the OpenCL language.

Forces computation on the GPU

Strategy of parallelization. The same parallelization strategy as the one presented in [ACF11] for an implicit Finite Element solver was adopted in this work, and tailored to our case. Thus, to obtain an efficient parallelization on the GPU, the object's mesh is not partitioned (as this would imply concurrent conflicts), but a two-step process is used instead, involving a parallel computation of forces divided into two tasks, each involving a set of kernels:

1. First, the forces applied on each node of a given element are computed and stored. This computation (kernel1) is done for each element of the object. 2. Second, for each node of the object (kernel2), these partial forces are summed to obtain the total forces applied on each node (involved by the different elements of the object).

Moreover, to ensure an efficient pooling of these kernels (i.e. an efficient use of cores), the number of nodes and elements may be chosen as a power of 2, involving grouping of kernels running in parallel of size a power of 2. Data structures. Specific data structures were defined for a discretization of the object into N elements involving m nodes, with n the number of nodes of each element and Nn the maximal number of neighbor elements for a node (see Fig. To be more efficient on the GPU, one may note that the float4 data type is used for the storage of partial and total forces. Moreover, basic arithmetic instructions are naturally coded into the hardware (using Arithmetic and Logic Units).

Communication. Fig. 1 illustrates how the global memory is accessed in reading for the initial volume's value of each element i (vol i in equation ( 11)). The float4 data type is again used. And to limit the number of accesses, only one kernel over 4 in a grouping of kernels (wrap) will read data from the global memory. In this example, kernel 0 reads initial volumes {vol 0 , vol 1 , vol 2 , vol 3 } and kernel 4 reads initial volumes {vol 4 , vol 5 , vol 6 , vol 7 }. Note that, the same strategy is used to read masses of nodes during the numerical integration scheme.

Parallel algorithm. Algorithm 1 presents the parallel algorithm for force computation -SOFA function addForce(). In this algorithm, function Force() enables the force computation F E (P i ) of each node P i in an element (equation ( 12)). It depends on constant values (Lamé coefficients l, µ; coefficients of interpolation functions a j , b j for j = 0..3 -equation ( 6)) and some variables (initial position and current displacements of nodes). This function is written using formal computation software and enables to easily generate equations (expressed as polynomials in function of U) of any kind of elasticity models. Consequently, we do not need to construct the matrix formulation of the problem of the form K U = F as for the classical FEM. • Consider node P 2 of element E 0 , and node P 1 of element

E 1 . We have index[0][2] = index[1][1] = 2.
Consequently, node P 2 in global indexation is a common vertex of elements E 0 and E 1 . • Then, the partial forces applied on this node involved by these two elements remain to be computed. Considering element E 0 (resp. E 1 ), ForceIndex[0][2] (resp.

Semi-implicit integration method on the GPU

Once the forces and the acceleration (according Newton's law) are computed, a numerical integration scheme is used to obtain the velocity (according to the acceleration) and displacement (according to the velocity) of the nodes.

Euler's semi-implicit scheme. Firstly, one of the simplest integration methods is considered, namely the Euler semiimplicit scheme (also called Euler symplectic), defined by:

V i (t + h) = V i (t) + h Vi (t) U i (t + h) = U i (t) + h V i (t + h) ( 17
)
where h is the time step and V i represents the velocity Ui .

Parallel algorithm. Algorithm 2 presents the parallel algorithm for the computation of the velocity and displacement for each node, using this numerical integration scheme (equation ( 17)). As previously for the computation of the forces, the float4 data type is used to store accelerations (Accel), velocities (V) and displacements (U) of nodes.

Algorithm 2 Computation of velocities and displacements on the GPU using Euler's semi-implicit method.

U[i] U[i] + h * V[i] ; 11: end for

The parallel implementation of this scheme is straightforward. However, to obtain a stable simulation, the time step h must be small, especially when stiffness increases. c The Eurographics Association 2012.

X. Faure & F. Zara & F. Jaillet & J-M. Moreau / An implicit Tensor-Mass solver on the GPU for soft bodies simulation 3.3. Implicit integration method on the GPU Euler's implicit scheme. To enable stability for larger time steps, an implicit integration scheme is mandatory like the Euler implicit scheme defined by:

V i (t + h) = V i (t) + h Vi (t + h), U i (t + h) = U i (t) + hV i (t + h). (18)
Noting M, F, V and U the mass matrix, forces, velocities and displacements vectors, respectively, this scheme may be reformulated in a linear form (Ax = b) [TPBF87, BW98]:

✓ M h ∂F ∂V h 2 ∂F ∂U ◆ | {z } A DV |{z} x = h F(t) + h 2 ∂F ∂U V (t) | {z } b ( 19
)
with DV = V (t + h) V (t) and ∂F ∂U , ∂F ∂V the Jacobian matrices encoding the variation of forces resulting from displacement and velocity changes. Consequently, we have to solve this linear system to obtain DV and next update the velocities and the displacements with:

V (t + h) = DV +V (t) U(t + h) = U(t) + h V (t + h). ( 20
)
Conjugate Gradient method. The Conjugate Gradient (CG) method, presented in Algorithm 3, is usually favored to solve the linear system (19) in a few iterations [BW98].

In this algorithm, the matrix ∂F ∂V is null since only internal forces (depending on displacement U) are considered. Thus in this algorithm, only ∂F ∂U is considered but never explicitly computed, directly computed h 2 ∂F ∂U V (t) or h ∂F ∂U DV . Algorithm 3 Conjugate Gradient algorithm to solve the Ax = b system from Euler's implicit integration scheme.

1: b h F(t) + h 2 ∂F ∂U V (t) 2: x 0 3: d r b 4: r 0 < r, r > 5: for i = 1 to n do 6: d f h 2 ∂F ∂U d 7:
A Md + d f Force() at line 9, to compute differential forces dF E (P i ) of each node P i in an element (equation ( 13)). This function depends on the same parameters, in addition to the time step h and the current nodes' velocity, to in fact directly compute h 2 dF E (P i ) V E (t) or h dF E (P i ) DV E (with V E (t) and DV E the elements of V (t) and DV corresponding to E).

For the parallel implementation of the CG method, we used the same parallel strategy as in [ACF11], but using the OpenCL language instead of CUDA. During the execution of this algorithm, only data a, b and r are transmitted between GPU and CPU to detect the end of the algorithm.

Finally, the velocities and displacements of nodes are updated in a similar way to Algorithm 2, but replacing the equations ( 17) by (20).

Results

For the simulations, a CPU Intel R Xeon R , 4 cores @3.07 GHz; and a GPU GeForce GTX 560, 2047 MB, 56 cores @1.620 GHz are used. Moreover, a 3D beam of size 5 ⇥ 1 ⇥ 1 m is considered, with a density r = 1, 000 Kg.m 3 and a Young modulus of 100 MPa.

Semi-implicit and implicit solvers

Precision. To test the precision of the two numerical integration schemes implemented, a traction test is performed on a Hookean beam discretized into 2,000 parallelepipeds, all split into 6 tetrahedra [MHS05]. A small deformation corresponding to 5 % of the initial beam size is imposed progressively with a displacement of nodes performed with a linear velocity. Tables 1 and2 show the running times to achieve the final state (directly linked to the given time step h) and the Young modulus computed according to the resulting strain. To ensure the stability using the semi-implicit scheme, the largest c The Eurographics Association 2012. time step is h = 0.0002, whereas in the implicit case, it is left unconstrained. However, to be complete, tests with the same time steps are presented is the first rows of Table 2.

The precision ratio (initial Young modulus compared to the average of the computed ones) reach 1.00611 (with a time step of 0.00012) and 1.00003 (with a time step of 0.093) for the semi-implicit and the implicit scheme, respectively. Therefore, despite its additional computation cost, an implicit scheme permits to achieve better precision at a reduced final time cost as the time steps may be considerably increased without loosing stability.

Memory. Table 3 shows the memory usage (linear according to the number of elements) and for each data structure, the percentage according to the total required memory. Here, a beam discretized into N = 60,000 tetrahedra is considered, involving m = 10,000 nodes, with a max number Nn = 24 of neighbor elements for any node. Moreover, an implicit scheme is employed, involving the additional vectors V or DV for the semi-implicit one. Note that P(t 0 ) and P(t) denote the initial and current positions of nodes in this table. Only 9.52 Mb are required for our similation, compared to the 560 Mb needed by the Abaqus FE software for the same scene. In the sequel, we will consider only the Euler implicit scheme, which provides very good results, both in precision and time.

Performances on the GPU

As previously, traction tests are preformed on beams (discretized into 12,000 tetrahedra) with an imposed 5 % deformation of their initial size.

Wrap's size. Fig. 4 presents running times (with time step h = 0.1) in function of the wrap's size, i. e. the number of cores involved at the same time (Fig. 1).

A stagnation in the execution time appears around 32 kernels in parallel for a total number of 56 cores. Indeed, the number of nodes and elements are chosen as a power of 2, involving an optimal wrap's size equal to a power of 2. In the remaining part of our analysis, we assume the wrap size to be 64, the nearest power of 2 greater than 56.

Number of elements per kernel. In Algorithm 1 for SOFA's function addForce(), we chose to execute in Task1 one kernel1 per element, and we used a similar strategy for SOFA's function addDForce(). A speedup of 25.5 is obtained for SOFA's corotational FEM implementation and of 29.5 for our TM, for a Saint Venant-Kirchoff's material beam composed of 307,200 elements.

Consequently, we only used 55 % of the GPU power (with 56 cores), due to synchronization between CPU and GPU in the Conjugate Gradient algorithm. Moreover, we may note that the running time on the CPU is linear according to the number of elements, whereas this time is affine on the GPU with a constant time cost for any number of elements. Therefore, our parallelization on the GPU is interesting for simulations with more than 1,000 elements due to the GPU's cost, judging from the elasticity model used (Hooke or Saint Venant-Kirchhoff).

Data transfer. The data transfer between local and global memory on the GPU is linear and corresponds for each element to constant values (Lamé coefficients l, µ; coefficients of interpolation functions a j , b j for j = 0..3 -equation ( 6)) and some variables (initial position and current displacements of nodes), requiring a total memory of 160 bytes per element.

Scaling study. Fig. 7 presents the computation time for Hookean and Saint Venant-Kirchhoff's elasticity models according to the number of tetrahedra in the mesh. The running times is nearly affine according to the number of elements.

Hookean and Saint Venant-Kirchhoff's behavior

Fig. 8 presents a rendered beam for Hookean and Saint Venant-Kirchhoff's material subjected to the gravitational force. Fig. 9 shows a rabbit interactively deformed with Saint Venant-Kirchhoff's elasticity model. Fig. 10 presents some results for torsion and bending of the beam considering both Hookean and Saint Venant-Kirchhoff's materials. In this article, a supplementary color MPEG file is also provided, which contains several views of the 3D simulation for better illustration of our results.

As expected, large deformations involved an increase of the volume for Hooke's elasticity model, whereas Saint Venant-Kirchhoff's material exhibits a much correct behavior.

Discussion

These results show that our parallelization of the Tensor-Mass formulation provided very positive results. Its efficiency is increased for large simulations including more than 1,000 elements, that is due to the implementation cost on the GPU.

We also note that the use of the Conjugate Gradient algorithm (needed in the implicit scheme), generated some synchronizations between the CPU and the GPU for data transfer. But, as the majority of the running time is dedicated to the addDForce() function, this does not really affect overall performances.

Consequently, to increase speedup, we have now to optimize the addDForce() function, and more specifically the DForce() function called inside it, which computes the differentials of the forces. As we ever said, this computation is based on formal computation to easily obtain their formulation and to generate the corresponding code. This is one of the reasons why a better speedup is obtained than for the corotational SOFA's FEM (see Fig. 6).

Conclusion

In this paper, we presented an implementation of the Tensor-Mass model on the GPU that considerably speeds up the simulation times. Comparisons between running times on the CPU and the GPU suggest that the parallel implementation of the model becomes interesting for increasingly complex computations.

Besides, we used a formal encoding for the equations of the forces (function Force()). However, this allows us to generalize the technique to almost any existing combination of constitutive law and mesh type, without extensive additional coding. This is particularly interesting as it will allow to implement and perform easily a variety of tests and comparisons for various materials on different solvers for an object within the same framework (SOFA, for example).

Moreover, the differentials of the forces are also be obtained formally (function DForce()), providing an easy way to implement Euler's implicit integration scheme, which requires these differentials. This computation would otherwise be very fastidious, as the energy equation must be derived twice; and finally, an implicit solver had never been implemented before for the Tensor-Mass model.

In addition, results are presented for tetrahedral meshes, but the derivation for other element types is straightforward, as long as interpolation functions L exist for each type (for example: triangle or quadrangle in 2D, and quadratic or cubic tetrahedron or hexahedron, prism, etc. in 3D). The only issue is that W E in equation ( 11) is no more constant over the element, necessitating computation that may be done formally as suggested above.

Consequently, we believe that a GPU version of the TM approach will help to make this model popular in the Computer Graphics community. Our implicit parallel solver for linear and non-linear behavior yields a stable version with high time reduction, in comparison with the explicit CPU implementation of the linear TM that may be found in the literature. Moreover, relatively similar simulation times are achieved in comparison to other existing FEM, but allowing better interaction control. This leads to an efficient alternative implementation for interactive physically-based simulation. Indeed, it combines advantages of discrete models (like Mass-Spring Systems) such as topological changes and interactions with the surrounding scene, and those of FEM, as it is derived as well from continuum mechanics.

  2-3 for an illustrated example): • index (of size N ⇥ n) stores the relationship between local and global indexation for each node of each element. For example, index[e][v] gives the global indexation for node v of element e. • PartialForce (of size Nn ⇥ m ⇥ 4) enables the storage of 3D coordinates of partial forces for each node considering its global indexation. Thus, PartialForce[][v] gives partial forces of node v of the object. • TotalForce (of size m ⇥ 4) stores the sum of the partial forces for each node considering its global indexation. Thus, TotalForce[v] gives the forces of node v of the object. • ForceIndex (of size N ⇥ n) stores the index of data structure PartialForce. For example, ForceIndex[e][v] indicates where the forces applied on vertex v involved by the element e are stored in Partial-Force.

Figure 1 :

 1 Figure1: Only kernels of a wrap (grouping of kernels) with a multiple of 4 as number will read data from global memory.

c

  The Eurographics Association 2012. Algorithm 1 addForce() function on the GPU. 1: {N : number of elements} 2: {n : number of nodes per element} 3: {m : total number of nodes} 4: {Nn : max number of neighbor elements for a node} 5: // Task 1 -Computation of partial forces 6: for e = 0 to N 1 do [e][v]][index[e][v]] Force(); 10: end for 11: end for 12: // Task 2 -Sum of partial forces 13: for i = 0 to m 1 do 14: // Execution of m kernel2 15: for j = 0 to Nn 1 do 16: TotalForce[i] TotalForce[i] + PartialForce[i][j]; 17: end for 18: end for Illustration. To illustrate this parallelization, let us consider a simple deformable object divided into 2 tetrahedra. Fig. 2 exhibits the local and global indexations of nodes. Fig. 3 presents the data structures resulting from this example.

  (a) Local vertex indexation. (b) Global vertex indexation.

Figure 2 :

 2 Figure 2: Simple object discretized into 2 tetrahedral elements (N = 2, n = 4, m = 5, Nn = 2).

ForceIndex[ 1 Figure 3 :

 13 Figure 3: Data structures for the GPU forces computation.

  1: {m : total number of nodes} 2: // Task 3 -Computation of velocities 3: for i = 0 to m 1 do ] + h * Accel[i] ; 6: end for 7: // Task 4 -Computation of displacements 8: for i = 0 to m 1 do 9: // Execution of m kernel4 10:

  The parallel computation of differential forces (SOFA's function addDForce()) is similar to the parallel computation of forces (SOFA's function addForce()) presented in Algorithm 1. The only difference resides in the use of function DForce() instead of

Fig. 5 Figure 4 :Figure 5 :

 545 Fig.5presents the variation of the number of elements treated in kernel1 and inside its equivalent kernel in function addDForce(), requiring additional loops inside them. The results validate our strategy, i.e. to treat only one element per kernel.

Figure 6 :

 6 Figure 6: Speedup between CPU and GPU for TM and SOFA's corotational FEM.

Figure 7 :

 7 Figure 7: Average time (in second) of one step as a function of the number of elements in the mesh.

c

  The Eurographics Association 2012.

Figure 8 :

 8 Figure 8: Rendered beam for Hookean and Saint Venant-Kirchhoff's material, and its initial state (from left to right).

  (a) Initial mesh (zoom) (b) Initial undeformed rabbit (c) Curved left ear (d) Little jump

Figure 9 :

 9 Figure9: Simulation of the deformation of a Saint Venant-Kirchhoff's material rabbit (initial 3D mesh courtesy of L. Stanculescu, with radius-edge ratio q < 1.0).

Figure 10 :

 10 Figure 10: Hookean versus Saint Venant-Kirchhoff's behavior for torsion and bending test an a beam.

Table 1 :

 1 Results for Euler's semi-implicit scheme.

	Time step	Computed Young's modulus (Pa)	Running time (s)
	0.00005	98,775,400	220
	0.00010	99,360,500	108
	0.00020	100,041,000	52
	0.00012	99,392,300	127 (average)
	Time step	Computed Young's modulus (Pa)	Running time (s)
	0.00005	96,177,200	841
	0.00010	97,839,800	418
	0.00020	98,910,300	205
	0.00012	97,642,430	488 (average)
	Time step	Computed Young's modulus (Pa)	Running time (s)
	0.040	99,993,400	5
	0.080	99,998,200	2
	0.160	99,998,500	1
	0.093	99,996,700	3 (average)

Table 2 :

 2 Results for Euler's implicit scheme.

Table 3 :

 3 Memory required for each data structure.

Table 4 :

 4 Table4presents running times (in second) for different parts of the simulation loop, in the case of Hookean and Saint Venant-Kirchhoff's (SVK) elasticity models. For each function, the percentage is given according to the total time. As expected, most of the execution time is consumed by the addDForce function. Execution times (in second) for Hookean and Saint Venant Kirchhoff's elasticity models.

		addForce	addDForce	CG	Total
	Hooke (s)	0.0150	0.1367	0.0059	0.1576
	Hooke (%)	9.53	86.75	3.72	100
	SVK (s)	0.0147	0.1865	0.0048	0.2060
	SVK (%)	7.14	90.53	2.33	100

c The Eurographics Association 2012.


Hu Zheng

Dong Wang

�� Jonathan Barés

��� Robert Behringer

Jamming by compressing a system of granular crosses

A disordered stress-free granular packing can be jammed, transformed into a mechanically rigid structure, by increasing the density of particles or by applying shear deformation. The jamming behavior of systems made of 2D circular discs has been investigated in detail, but very little is known about jamming for non-spherical particles, and particularly, non-convex particles. Here, we perform an experimental study on jamming by compression of a system of quasi-2D granular crosses made of photo-elastic crosses. We measure the pressure evolution during cyclic compression and decompression. The Jamming packing fraction of these quasi-2D granular crosses is φ J � 0.475, which is much smaller than the value φ J � 0.84 for-2D granular disks. The packing fraction shifts systematically to higher values under compressive cycling, corresponding to systematic shifts in the stress-strain response curves. Associated with these shifts are rotations of the crosses, with minimal changes in their centers of mass.

Introduction

Granular materials can be easily found in nature and in industrial applications. They have gas-like, liquid-like or solid-like properties depending on their density or shear stress [1][2][3][4]. Many industrial devices exploit these state properties, and a number of geo-hazards involve uncontrolled state changes. For example, in controlled flow in an industrial process, such as in a hopper, gas-like or liquidlike states are desirable [5]. On the other hand, most geo-hazards, like landslide or debris flow occur when a jammed solid state transitions to a liquid-like state [6,7]. Hence, it is very important to understand the granular phase properties, and the transitions between fluid and solid, i.e. unjammed or jammed, states.

Characterizing the jamming transition in disordered granular systems has attracted significant interest [2,[8][9][10][11][12]. Liu and Nagel [8] hypothesized a zero temperature jamming phase diagram. In a space of shear stress, τ and packing fraction, φ, they proposed jammed states for τ < τ Y (φ), where the yield stress curve, τ Y (φ), is a monotonically increasing curve that intersects the φ axis at φ J . This proposal was substantiated by O'Hern et al. using numerical simulations of frictionless discs. These studies used protocols that produced macroscopically isotropic states. Recently, Bi et al. [2] showed that for frictional disks, jammed states occur below the frictionless isotropic φ J . They used a protocol consisting of shear strain applied � e-mail: hz64@phy.duke.edu �� e-mail: dw148@phy.edu.ed ��� e-mail: jb@jonathan-bares.eu ���� e-mail: bob@phy.duke.edu to an initially unjammed (stress-free) isotropic state. In a range φ S < φ < φ J , they found that increasing shear strain led from unjammed states to fragile and then robustly shear jammed states, where the last two states are characterized by non-zero shear stress. The isotropic jamming point of a two dimensional frictionless disc granular system is φ � 0.84 [9,13].

However, most of the experiments and simulations of jamming have focused on simple geometry particles like discs and spheres or at least on convex particles. Hence the focus of the present experiments is on non-convex crossshaped particles. Role of particle shape has been recognized as significant in controlling the microstucture of various granular materials [1,[START_REF] Duran | Sands, powdersm and Grains: An Introduction to Physics of Granular Materials[END_REF][START_REF] Wood | Soil Behavior and Critical State Soil Mechanics[END_REF][START_REF] Azéma | [END_REF]. For instance, Athanassidis et al. [17] have shown that the macroscopic mechanical properties of a granular material can be controlled by tuning the grain shape via a 3D printing approach to make different shape particles including tetrahedra, octahedra, and dodecahedra [17]. And several authors have shown that 3D systems of rods [18] or of non-convex particles, such as staples, or star-shaped particles have interlocking or entangled states that are mechanically stable [19][20][21]. For instance, star-shaped non-convex particles, i.e. hexapods are interesting to both granular science and architecture design. These particles form stable granular aggregates by pouring [22]. Zhao et al. [21] conducted laboratory experiments to show the stability and structure of poured hexapod packings by varying the arm length and coefficient of friction between particles. However, experimental methods generally cannot detect the internal details of such packings, or the forces between particles, which The photo-elastic method, used here, is an effective experimental technique to study 2D granular materials [10,[23][24][25][26] because it provides access to the forces, and hence local pressure on each particle. We combine this technique with direct visualiation of particle displacements and rotations. The latter is particulary important for crosses as they gradually evolve to lower stress states at a given φ under cyclic compression/decompression. We apply strain with a bi-axial apparatus [2,27] which allows controlled strains in two directions, as in Fig. 1(c) and (d).

Experimental techniques and protocols 2.1 Experimental techniques

Specifically, the experiments are conducted using photoelastic crosses in a two dimensional bi-axial apparatus. We cyclically compress this system through the isotropic jamming point, and then decompress the system to the original boundary configuration. The geometry of the crosses, cut out of a 6.35mm thick photo-elastic sheet by water-jet cutting, is shown in Fig. 1 (a). The length of each arm is 18cm, from the center of cross to the edge. The end of each arm is a semicircle and the corner between two neighbor arms is arc-shaped to avoid force concentration. The radii of the semicircle and arc is 3cm. The friction coefficient between crosses is ∼ 0.7.

The particle photo-elastic property allows us to measure the forces acting on them via digital images of the fringe patterns, as typified by Fig. 1 (b). In particular, the concentrated intensity variation close to a contact point gives a measure of the force at a contact, and helps determine the existence of contacts. The density of fringes corresponds to the magnitude of the force: the more fringes, the bigger the force. Fig. 1 (b) shows a photo-elastic image of a cross particle subjected to compression in all four directions, as indicated by the arrows. We use the 'G 2 ' method [23] as measure of the fringe density, and hence of the forces acting on a particle. The G 2 value is defined as the pixel-wise intensity gradient in the photo-elastic digital image. Moreover, the actual contact force can also be calculated in principle based on a complex fringe patternfitting algorithm pioneered by Majmudar [27]. In this paper, we only use G 2 as a measure of the pressure of each particle, which has been demonstrated to be a one-to-one function of pressure at the particle scale [28].

Experimental apparatus

The experimental apparatus performed bi-axial and uniaxial compression. An illustration of the system is shown in Fig 2 (a). The photo-elastic cross particles are placed on a transparent Plexiglas table in a horizontal plane, such that gravity plays no role in the granular behavior, other than through weak friction with the base (0.36 ± 0.04). Four UV tube lights illuminate the UV-fluorescent ink on the top of each cross (Fig. 2 (b)), allowing us to track the position and orientation of each cross. An LED array provides homogeneous light from below. A circular polarizer sheet on the top of the LED array polarizes the light. A second crossed polarizer is rotated into or out of position in front of a camera, providing normal and polarized images (Fig. 2 (c) and(d)). The camera (Canon R � 70D) has a resolution of 3648 × 5472 pixels, corresponding to ∼ 31, 000 pix 2 per cross.

Protocol

We compressed the system quasi-statically in steps of 1mm. After each step, we obtained three images, one with only UV light, one polarized and a normal image with light from the LED array. The initial size of the shear box was 47.2 × 39.3cm 2 . 270 crosses were placed in the shear box. Before compressing the system, we ensured that there was no initial stress in the system, which meant that no fringes were visible in the polarized image of the initial state.

Results and discussion

We focus on the jamming transition where we contrast system properties that were achieved by compressing either uniaxially or isotropically, i.e. biaxially. Below, we distinguish these two protocols by unicompression or bicompression. For these two protocols, we determined the contact number, the force chain network, and the system pressure, and used these properties to characterize the jamming transition.

The system pressure (G 2 ) evolution as a function of the packing fraction for multiple compression cycles is shown in Fig. 3. The blue line corresponds to unicompression and the red one to bicompression. The green region indicates the transition area of the granular crosses from an unjammed to a jammed state. Manifestly, the jamming point is not a unique function of packing fraction, and it also depends on the compression protocol. The difference between unicompression and bicompression can be seen from the inset images in Fig. 3. The force chain network generated by unicompressing is more anisotropic than the bicompression case, reflecting the more anisotropic character of unicompression. The difficulty of rearrangement in the cross packing also caused a relative fast G 2 development by unicompressing with a slope of about 7.8. For the bicompression system, the slope of G 2 is about 7.2, and the system evolution is smoother.

The crosses move little during cyclic compression, G 2 decays over multiple cycles, as seen in 2D disc system [13] and in simulations of frictionless spheres [? ]. However, G 2 appears to reach a plateau, within moderate fluctuations, for large enough N c for bicompression, although not necessarily for unicompressions. Fig. 4 (b) and(c) show the force chain networks at the packing fraction, φ = 48%. The photo-elastic response at cycle 43 is slightly weaker than after the first cycle.

We compared the position and orientation of each cross in the unicompression system, as shown in Fig. 5, to see that the position of each particle changes very little, even after 43 compression cycles; consequently, the coarse-grain density also changes little over the 43 cycles. However, the particle orientations change substantially. It is reasonable to conclude that the pressure decay at the maximum packing fraction is strongly correlated with changes in the particle orientations.

Conclusion

We have experimentally investigated the jamming transition of a packing formed from 2D photoelastic crosses subject to multiple cycles of compression that is either uniaxial or biaxial. The jamming point of this system is between 47.5% and 48.0%, and the specific jamming density is sensitive to details of the compression protocol, including the cycle number and whether the compression is unixial of biaxial. The rate of increase of G 2 with φ is higher for unicompression than for bicompression. Several measures are sensitive to the cycle number, including the jamming point and the pressure at a given packing fractiontwo properties that are correlated. There is little change in the positions of individual particles, but there is significant rotation of crosses during the course of multiple cycles, and this rotation provides a mechanism for stress relaxation.

Figure 1 .

 1 Figure 1. (a): Sketch of cross particle geometry. (b): Photoelastic response generated by compressing a cross particle along each arm. (c): Schematic of uniaxial compression, using the biaxial apparatus. (d): Sketch of biaxial compression. The gravity direction is given by g .

Figure 2 .

 2 photoelastic particles

Figure 3 .

 3 Figure 3. System pressure, expressed by G 2 in arbitrary units, vs. compressive strain for multiple cycles. Insets: Image indicated by red arrow is a photo-elastic image at the jamming point for unicompressing while the one indicated by blue arrow is for bicompression. The force chain networks of the two images indicate that a more isotropic network occurs for the bicompression case. The green area shows the jamming transition zone.

Figure 4 .

 4 Figure 4. (a): Maximum G 2 at the highest packing fraction of each compression circle. Note that the maximum compressions and vertical scales are different for the two compression protocols. The maximum G 2 decays with N c , the number of compression cycles for both experimental methods. Alternatively, this can be seen as an increase in φ J as a function of N c . (b): Photo-elastic image when the system firstly reached the maximum packing fraction, φ max = 48%. (c): Photo-elastic image at the maximum packing fraction after 43 shearing cycles. The white bar in (b) indicates 10cm. (b) and (c) share the same scale.

Figure 5 .

 5 Figure 5. Positions and orientations of the crosses for unicompression initially (red crosses) and after 43 compression cycles (green crosses).

Acknowledgements and References

This work is supported by NSF Grant No.DMR1206351, NASA Grant No. NNX15AD38G, the W.M. Keck Foundation and a RT-MRSEC fellowship. HZ thanks NSFC Grant No. 41672256 and the NSF(Jiangsu, China) Grant No. BK20140845 for financial support.


Aghil Abed Zadeh

Jonathan Barés

Robert P Behringer

Avalanches in a granular stick-slip experiment: detection using wavelets

Avalanches have been experimentally investigated in a wide range of physical systems from granular physics to friction. Here, we measure and detect avalanches in a 2D granular stick-slip experiment. We discuss the conventional way of signal processing for avalanche extraction and how statistics depend on several parameters that are chosen in the analysis process. Then, we introduce another way of detecting avalanches using wavelet transformations that can be applied in many other systems. We show that by using this method and measuring Lipschitz exponents, we can intelligently detect noise in a signal, which leads to a better avalanche extraction and more reliable avalanche statistics.

Introduction

Many systems, under slow loading, exhibit intermittent and discrete events of broad sizes, known as crackling noise [1,2]. This behavior has been experimentally investigated in a wide range of physical systems [3][START_REF] Planet | the Subsurface: Mixing, Spreading and Reaction in Heterogeneous Media[END_REF][START_REF] Barés | [END_REF] from earthquakes in seismic activities to Barkhausen noise in ferromagnets [6,7]. Slowly sheared granular media are examples of these systems. They respond to shear stress with discrete slip avalanches, also called stick-slip events [8][9][10][11]. The avalanches statistics of these systems and their universal and non-universal quantities are key properties of these systems. In this paper, we introduce a method for avalanche detection based on wavelet transformation to extract more reliable statistics. This method can be used for many other systems for a better avalanche detection.

As shown in figure 1, to observe the physics of stickslip for slowly sheared granular materials, we have developed a novel apparatus, which is capable of detecting avalanches both at the global and local scales. In our experiments, a constant speed stage pulls a frictional slider of length 25 cm which rests on a vertical bed of bi-disperse cylindrical photo-elastic particles with diameteres of 0.4 cm and 0.5 cm (number ratio of 2.7) in a 2D granular system. The verticle walls are low frictional to reduce the efect of boundary friction. The slider is designed to have very high friction with the top layer of granular media. The stage is connected to the slider by a spring, the stiffness of which can be tuned. We measure the force on the spring by a force sensor attached to it, which provides a global measure of shear stress. As the system is intended to be at constant pressure, the slider is not fixed in any of 2D directions and may result in few degrees of inclination during stick-slip that we think is not significantly important. The photo-elastic response also provides us a measure of local stress; however, in this paper, we focus on the global response of the system. For a given pulling stage speed v, spring stiffness k, and granular bed height l, the evolution of the force F, applied to the slider is measured with a frequency of 100Hz in our experiment. To quantify the global avalanche statistics based on F(t), we use specific signal processing. In experimental data, noise is unavoidable, but the way noise is detected or reduced can be crucial to quantitative results. Conventionally, these noise spikes are removed using Fourier transform domain filters, including averaging over many data points. The latter is a low-pass fil- ter in Fourier transform domain and increases the signalto-noise ratio. However, the trade-off is loosing the edge information. Moreover, if the duration of an event in a signal is of the order of the duration of a noise spike, the event may break into several events and make them hardly distinguishable; or some noise spikes may be detected as events. These flaws can be reduced by applying a wavelet transform domain filter that selectively distinguishes noise spikes from events.

In the first part of this paper, we summarize a conventional way of avalanche detection in a typical signal of our experiment. We discuss how several chosen parameters in the signal analysis process can affect the statistics. Then, we introduce a wavelet transform domain filter and how to selectively detect noise spikes locally using Lipschitz exponents [12], without loosing other signal information. This method may provide a tool for a better avalanche extraction in many other experimental systems.

Conventional method

The slow shearing of a granular materials leads to intermittent slow loading and fast unloading of the medium that is called an avalanche. Figure 2a shows a typical signal F(t) from the force gauge. The fast drops in F (like at time t = 13s) correspond to slip avalanches of the slider and similar behavior in observables have been detected in many physical systems [START_REF] Barés | [END_REF]8,13,14]. Hence, to detect an avalanche, we need to find where the slope of the force is negative. However, if we use point to point derivatives, as shown in figure 2b, the presence of noise may make it impossible to find the avalanches. To overcome this problem conventionally, as depicted in figure 2c, the signal is smoothed and then a threshold τ is chosen for finding where it crosses the smoothed data. This defines the beginnings and endings of avalanches. To smooth the signal from our data, we applied a convolution as presented in equation 1, in which h is the derivative of the force signal F and g is a normalized discrete Gaussian function. This method requires two variables, the standard deviation of the smoothing function, σ, and the threshold, τ. The statistics of avalanches vary by changing these two variables.

(h * g)[m] = n t=1 f [t]g[m -t]
(1)

g[t] = 1 σ √ 2π exp (- t 2 2σ 2 ) ( 2
)
Figure 3 shows the probability density function of slip sizes for different values of the threshold τ and standard deviation σ ( as in fig. 2c). The slip size is the displacement of the slider during an avalanche, which is calculated using the force signal ∆x = ∆F/k, in which k is the stiffness of the spring. As shown in the figure, the probability density function (PDF) is affected by τ and σ. Indeed, both the exponent and the non-universal quantities such as lower cut-offs and prefactors depend on these two parameters. As large avalanches are easier to detect and less affected by low-pass filters, the statistics of them do not change dramatically by changing these two variables. However, small avalanches statistics are more affected by these parameters as they contain more high frequency information. This variation can also be seen in other avalanche properties such as the average shape of avalanches (not shown here).

Wavelet transformation method

Wavelet transforms have been widely used in signal processing in computer science for detection of different types of singularities in a signal [12]. A wavelet transformation of a signal is simply its convolution with a function g as in equation 1. The convoluting wavelet that we use here is proportional to the first derivative of a Gaussian function, as in equation 3, which is a widely used wavelet for edge detection problems [12].

g[t] = t σ 2 √ 2π exp (- t 2 2σ 2 ) (3)
Each local maximum in the signal wavelet transformation cooresponds to a "singularity". The dependence of the local maxima on σ determines the type of the singularity [12]. Typical noise is associated with a back and forth peak which resembles a Dirac delta function, whereas, avalanches are sharp unique variations (slow loading and One of local maxima is marked by red circles. We track the position of such local maxima to detect whether it corresponds to an avalanche or not. The avalanche detection is based on the evolution of the maxima position vs σ. rapid unloading), which resemble step function. This difference can be captured by well chosen wavelets. Figure 4 depicts the wavelet transformation of the force signal F (of fig 2a) as a function of σ. We track each local maximum, one of which is marked red in figure, vs. σ to detect the singularity type.

The evolution of local maxima vs. σ for the figure 2a is shown in figure 5. Each line in figure 5 belongs to a local maximum. The slope of each line is called a Lipschitz exponent of the singularity cooresponding with the maximum of the signal derivative. This exponent is non-negative for an edge function and negative for a Dirac delta function [12]. Hence, it is also non-negative for a singularity corresponding to avalanche events in our case, and negative for those of noise. Moreover, to detect the avalanches more precisely,we need to have non-negative Lipschitz exponent and a clear exponential decaying tail. By using the Lipschitz exponent and detecting the decaying tail, we can extract the maxima that correspond to the avalanches, which are marked blue in figure 5.

This wavelet transformation method can improve on avalanche detection over the conventional methods, without losing high frequency information. It also provides a better detection of small avalanches. Figure 6 shows the detected avalanches in a force signal F. The avalanches are marked as blue and each avalanche is a connected set of blue points in the figure. Using this method, we can then calculate the probability density function of avalanche sizes and other avalanche statistics. Figure 7 shows the PDF of avalanche slip size for instance. This figure is calculated using several thousand detected avalanches for a spring stiffness of k = 70N/m and stage speed of v = 1mm/s. The distribution obeys a power-law for about two decades with a clear lower cut-off and an exponential upper cut-off at slip size 2 × 10 -5 m and 2 × 10 -3 m respectively. Because of better detection of avalanches, the wavelet transformations method provides a smaller lower cut-off and larger number of power-law decades, compared to the conventional method (fig 3). This local noise detection and removal provides a more reliable exponent calculation in the PDF of avalanche sizes over the conventional method, which usually uses global smoothing of the data.

Conclusion

Our experimental apparatus enables us to study the avalanches in a granular stick-slip experiment. The slow shear of granular system is followed by intermittent fast unloading of shear stress that is measured by a force signal. This force signal is conventionally smoothed (low-pass filter) to increase the signal-to-noise ratio for avalanche detection. However, this conventional method changes the statistics of avalanches as it removes all low frequency information. We introduced a smarter way of avalanche extraction from measured quantities using wavelet transformation. Similar methods have been tested in many other problems and here it provides a better avalanche detection and consequently more reliable statistics, by detecting the noise locally in the signal. The method introduced here can be applied in many other systems for a more precise avalanche extraction.

Figure 1 .

 1 Figure1. Side view of apparatus. A bi-disperse photo-elastic granular media is sheared by a slider that is moved by the force gauge connected to a spring. The slider exhibits a stick and slip dynamics. The force on the slider is recorded via a force gauge and the internal structure and stress is imaged via high speed cameras (not treated in this paper).

Figure 2 .

 2 Figure 2. a) A typical force gauge signal F(t) from the experiment with spring stiffness of k = 70N/m and stage speed of v = 1mm/s acquired at 100Hz. b) Point to point derivative of force signal. The signal-to-noise ratio decreases dramatically when differentiated. c) The smoothed derivative of force signal. Smoothing is applied by a convolution of part b with a Gaussian function of standard deviation σ = 0.05s. The threshold τ (red line) determines the beginnings and endings of avalanches.

Figure 3 .

 3 Figure 3. Probability density functions (PDF) of slip sizes during avalanches for: a) Several different values of threshold τ (for fixed σ = 0.05) and b) Several different values of standard deviation σ (for fixed τ = 0.1) . These two parameters affect the statistical measurements of the force signal.

Figure 4 .

 4 Figure 4. Wavelet transformation of the force signal F of figure 2a with the wavelet given in equation 3 as a function of σ.One of local maxima is marked by red circles. We track the position of such local maxima to detect whether it corresponds to an avalanche or not. The avalanche detection is based on the evolution of the maxima position vs σ.

Figure 5 .

 5 Figure 5. Position of the local maxima (as in fig.4) values as a function of σ for a typical force signal F. Maxima that are detected to belong to avalanches are marked blue. See text for more details.

Figure 6 .

 6 Figure 6. The detected avalanches in a force signal sample using the wavelet transformation method. Avalanches are marked blue. This method has the advantage of better detection of small avalanches and ignoring the noise spikes occurring during avalanches. These noise spikes can break the avalanches into several ones in the conventional method.

Figure 7 .

 7 Figure 7. Probability density functions (PDF) of slip sizes during avalanches, detected by the wavelet transformation method. The PDF obeys a power-law for more number of decades compared to conventional method (fig3), as a result of avalanche detection using wavelets.

e-mail: aghil.abed.zadeh@duke.edu e-mail: jonathan.bares@umontpellier.fr e-mail: bob@phy.duke.edu

Acknowledgements and References

We acknowledge support from NSF-DMS-1248071, NASA Grant NNX15AD38G, and William M. Keck Foundation.