MULTIBODY DYNAMICS 2007, ECCOMAS Thematic Conference C.L. Bottasso, P. Masarati, L. Trainelli (eds.) Milano, Italy, 25–28 June 2007

NUMERICAL STRATEGIES AND SOFTWARE ARCHITECTURE DEDICATED TO THE MODELLING OF DYNAMICAL SYSTEMS IN INTERACTION. APPLICATION TO MULTIBODY DYNAMICS Fr´ed´eric Dubois⋆ and Mathieu Renouf† ⋆ LMGC/Universit´ e

Montpellier II - CNRS - UMR5508 CC048 - place Eug`ene Bataillon, F34096 Montpellier cedex - France e-mail: [email protected], webpage: http://www.lmgc.univ-montp2.fr/∼dubois † LaMCoS/INSA

Lyon - CNRS - UMR5259 18-20, rue des sciences, F69621 VILLEURBANNE cedex - France e-mail: [email protected] webpage: http://mathieu.renouf.site.voila.fr

Keywords: DEM, FEM, Multicontact, Dynamics, Multiphysics Abstract. Many mechanical, tribological and geophysical problems bring difficulties such as large deformations, fracture, multiple localisations or wear. To fill the lack of knowledge and to bring information in fields where experimental approaches are limited (for example, contact interface measures), numerical tools have been developed and adapted. But if for some applications efficient numerical methods are now well identified, it still lacks tools able to cope with a class of multiproblems: multibodies, multicontacts, multiscales and multiphysics. Considering this kind of problems, needs to widen the purpose of classical modelling approaches, mixing continuous medium and discrete medium descriptions. The relevance of continuous models at the structure level and the discrete models at the microstructure level shows the necessity to converge to hybrid models able to switch between continuous and discrete descriptions if necessary. Moreover, it needs to improve theoretical and numerical related tools. The aim of the proposed paper is to present the state of a common work in view to build an unified modelling framework able to cope with extended discrete element method: hybrid DEM/FEM models, multiphysics applications, constrained/unconstrained modelling of contact, implicit/explicit time integrators, various bulk behaviours and interaction laws, etc. Various applications are chosen to illustrate the interest of such unified modelling framework in mechanics as well as in geophysics and tribology.

1

Fr´ed´eric Dubois and Mathieu Renouf

1 INTRODUCTION In the present paper, focus is made on the way of modelling divided materials or structures assimilated to a large collection of complex mechanical components in interaction. Among mechanical systems, a large number presents naturally a divided feature: granular media, masonries, metal at the grain level, etc. Some of them, due to external loading, may evolve from a continuous state to a, globally or locally, divided one: fracture, wear, cracks, etc. The reverse evolution, e.g. when a divided material evolve to a continuous one, may also happens: ceramics, Tribological Surface Transformation. Finally, one can introduce also a fictitious divided feature in the mechanical system such as in domain decomposition approaches or multigrid methods. Among multibody systems, two classes could be defined. The first one concerns systems for which the number of bodies is greater than the number of contacts, this includes dynamics of machines [26] or animated computer graphics [12]. The second one concerns systems for which the number of bodies is less than the number of contacts. In this case, the name of multicontact may be more relevant than the name multibody. In this paper, the focus is made on this secondclass only. Concerning a relevant framework able to model the first class of problems, one can consider for example the SICONOS platform (http://siconos.inrialpes.fr). The formulation of a problem involving numerous complex mechanical components in interaction raises different issues, which can be divided in two subsets: • the system description based on modelling tools. It includes: – a suitable description of each component of the system, including a modelling framework (the Lagrange equations in mechanics, the heat transfer equation , etc.), the choice of the behaviour (rigid, deformable with or without multiphysics coupling ), the space discretization, etc. – a suitable description of potential interaction zones (i.e. contactors). For example, it concerns the geometric description of an interaction zone of a component (i.e. boundary) and the mappings between the degrees of freedom of the components and the degrees of freedom of their interaction zones. – handling of interactions: detection, computation of some information (local frame, relative motions, topology of interaction network, etc.). – a suitable description of the interactions behaviours between the different components. • the numerical strategies to solve the problem based on simulation tools. It includes: – time evolution as Event-Driven or arbitrary Time-Stepping; – time integration; – resolution of the interaction problem (constrained or unconstrained global problem and smooth or nonsmooth contact laws). Among the different ingredients of the strategies, two of them appear as the most important: the time integrator (implicit or explicit) and the way to deal globally with interactions (constrained or unconstrained). Note that these two points are improperly presented 2

Fr´ed´eric Dubois and Mathieu Renouf

as linked, in the way that explicit (resp. implicit) time integrators are often used with unconstrained (resp. constrained) contact treatment. Mixing is possible as presented later. What ever happens, an explicit time integrator depends on dynamics through stability conditions. The multidisciplinary aspect of modelling such complex systems requires the integration, within a single environment, of various techniques as diverse as contact modelling, multiphysics modelling, computational geometry, intensive computation, etc. It starts with the developments of Discrete Element Methods, as the pioneer work of Cundall who developed the Distinct Element Method [10] to study blocky rock systems by considering each block as a rigid element. Later, the method was naturally extended to the simulation of granular media [11]. The Molecular Dynamics (MD) approach [3] is another contribution to DEM, different of the Cundall approach. MD, which consists initially in simulating the dynamics of atom and molecule in order to deduce macroscopic properties of the matter, was extended to the simulations of non-atomistic particles such as the study of granular flow [34]. DEMs have benefited of numerous improvements. Among them, the Contact Dynamics method initially developed by Moreau [23] can be emphasis. It appears as a different method, essentially because of a specific nonsmooth dynamical framework based both on measure differential inclusions formulation and on a relevant numerical strategy: implicit time-integration scheme, implicit treatment of the contact problem. It allows the numerical treatment of nonsmooth events such as impacts, Coulomb friction and unilateral constraints. Further works lead to the extension of the CD method to multicontact simulations of collections of deformable bodies [15]. Less than antagonist, these different methods are complementary, and are often related to a specific context. In sake of generic and scalable concerning modelling (space discretization, bulk behaviours, interaction laws, etc.) and simulation (time evolution, time integration, solving methods, etc.), relevant software architecture is necessary. A framework has been built on a dedicated class design; it will be presented in this paper. Actually, the LMGC90 platform (http://www.lmgc.univmontp2.fr/∼dubois/LMGC90) is based on this framework. The paper is organized as follows. Section 2 describes the mechanical framework of the extended DEM used in the paper. Section 3 is dedicated to the platform architecture. Section 4 presents some multiphysics extensions. Then several numerical studies are presented in Section 5. Finally, Section 6 provides concluding remarks. 2 MECHANICAL FRAMEWORK OF DISCRETE ELEMENT METHODS The purpose of this section is to review various aspects of DEMs relative to modelization and simulation. It will prepare the next Section devoted to the platform architecture. Classically, Discrete Element Method are based on some hypotheses as: the contact area is small behind the size of the component, the deformation concerns only the contact point neighbourhood or the interactions are binary (no long distance effects). These allow to consider that locus of interaction may be supposed as punctual, that components of the system may be considered as rigid or that interaction model depends only on the properties of the two concerned components through an interaction law giving the contact force as a function of the motion of the two components. 3

Fr´ed´eric Dubois and Mathieu Renouf

The pioneer DEMs ([11],[3]) are relying on these hypothesis. Furthermore it supposes that interactions, even connected by a component, are independent. It allows computing contact forces explicitly as function of contact deformations. These entire hypotheses are questionable, and various approaches are now available to extend the use of DEMs. Contact Dynamics is one of them, but there is huge amount of noticeable contributions particularly concerning the space discretisation of the system Munjiza [25] or Li et al. [20].

2.1 Modelling background 2.1.1 Dynamics Assuming a smooth evolution of a multicontact system allows describing the motion of each mechanical component by : ˙ = Fext (t, q, q) ˙ + R, M¨ q + Fint (t, q, q)

(1)

˙ the internal force and the nonlinear inerwhere M represents the inertia matrix, Fint (t, q, q) ext ˙ the external forces while R represents the contact forces. The vector tia terms, F (t, q, q) ¨ denote respectively the q represents the vector of generalised degrees of freedom and q˙ and q generalised velocities and the accelerations. Considering a linear deformable mechanical component, a classical Finite Element space discretization may be performed. In this case the internal forces can be decomposed as: ˙ = Cq˙ + Kq, Fint (t, q, q)

(2)

where C and K represents respectively the viscosity and the stiffness matrices. Considering a non-linear deformable mechanical component is detailed in Jean [15] Considering a rigid component, a more suitable dynamics equation may be used introducing the two vectors u and ω, related respectively to the translation and the rotation velocities of the mass-centre. The Eq.(1) is replaced by the well-known Newton-Euler system of equations:  ˙ +R Mu˙ = Fext (t, q, q) . (3) Iω˙ = −ω ∧ (Iω) + MFext + MR where M and I represents respectively the mass and the inertia matrices. MFext and MR represent respectively the momentum due to external and contact forces. Note that for 2D components or 3D components with geometric isotropy, the vector ω ∧ (Iω) is equal to zero. Assuming a smooth dynamic motion for each component of the system is relevant if one takes into account the deformation of contact neighbourhood and if the time discretization is sufficiently fine to describe collisions (i.e. is a fraction of collision time). Otherwise it is necessary to adapt the formalism in order to cope with discontinuities arising from contact-impact ([22],[24]). Since in multi-contact systems, multiple collisions and other non-smooth phenomena are expected, the velocity may be discontinuous and the acceleration can not be defined as the usual second time derivative of q. Assuming q˙ to have bounded variation entails the exis˙ denoted dq. ˙ Consequently the classical equation of motion tence of a differential measure of q, is reformulated in terms of a measure differential equation: ˙ ˙ Mdq˙ + Fint (t, q, q)dt = Fext (t, q, q)dt + dhR, 4

(4)

Fr´ed´eric Dubois and Mathieu Renouf

where, as previously, M represents the inertia matrix and Fint , Fext the internal and external forces. dt is the Lebesgue measure on the real space R and dhR is a differential measure of the cumulated impulsions hR experienced by the system from the initial instant. If a collision occurs in the system at some instant tc , the measure dhR possesses an atom of the form pδ(tc ), where δ(tc ) denotes the Dirac measure at point tc while p consists of the generalised components of the contact percussions arising in the system at this instant. In the smooth case the classical form is recovered. 2.1.2 Local-Global mapping The Eqs.(1,3 or 4) involve only global variables (component level). To deal with multicontact problem it is necessary to introduce local variables (contact level): The local relative velocities vα (α, index of the contact) concatenated in a large vector v, and the local contact forces rα in a vector r. The corresponding global variables are the generalised velocities of the system ˙ and the contact reactions R applied to all the bodies. (collection of components) q, (a)

y 11111111 00000000 00000000 11111111 00000000 11111111 00000000 i 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 11111111 00000000 00000000 11111111 q y 11111111 00000000 11111111 00000000 11111111

t ij nij

111111 000000 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000 111111 000000j 111111

qx

H Ry Rz

1111 0000 0000 1111 0000 1111 0000 1111 0000 1111 0000 1111 0000 1111 0000 1111

qz o

(b)

rn Rx

1111111 0000000 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111 0000000 1111111

rt

x

Figure 1: Representation of: (a) the global and local frames and (b) the linear mapping

Local variables are related to global ones via the linear mapping H(q)  R = H(q)r . v = H∗ (q)q˙

(5)

where H∗ is the transpose of H. The mappings are built using kinematic considerations such as the local frame defined at each contact point. Note that the mapping H is equal to the mapping G defined in [7]. The global to local mappings can be written straightforwardly in terms of measures. 2.1.3 Interaction laws One may consider various interaction laws (see for example [6], [3], [11] and [1]). Among them, one can consider unilateral contact law, Coulomb’s friction law, cohesive law, etc. This is not the topic of the present paper to detail these laws. However, it can be underlined that classical DEMs replace the geometric condition of impenetrability by some close-range repulsion laws; the abrupt feature of dry friction are similarly smoothed. This is the price paid to obtain dynamics tractable by standard integration techniques. The drawback is that, for the sake of precision, very steep repulsion laws have to be 5

Fr´ed´eric Dubois and Mathieu Renouf

used, generating stiff differential equations, which require very short step-length for their integration, possibly resorting to the introduction of artificial damping and artificial inertia to secure numerical stability. The Contact Dynamics approach can face unilaterality and dry friction without smoothing approximation. However, some classical laws need to be reconsider as velocity complementarity condition instead of the Signorini - gap complementarity condition and impact laws need to be introduced when bodies are modelled as perfectly rigid. See the work of J.J. Moreau for details [22, 23, 24]. 2.2 Numerical Strategies 2.2.1 Overview Pioneer DEMs ([3, 11, 14]) are based on very restrictive hypothesis, and they refer to explicit time integration scheme and explicit treatment of contact using smooth interaction laws (see further for details). This feature implies a permanent correction of an unbalance situation of the system, where kinematic constrains are hard to verify, even if decreasing time-steps to zero tends to uncoupled interactions over a time step. This kind of methods refers to unconstrained problems. explicit treatment of contacts

implicit treatment of contacts kinematic constraint

no h << t c classical dynamics

smooth contact law

Cundall Iordanoff

yes h >> t c

h << t c

h >> t c

classical dynamics

nonsmooth dynamics

nonsmooth contact law

smooth contact law

Carpentel et al Kishino

nonsmooth contact law

smooth contact law

Moreau Jean

Figure 2: Classification

Using an implicit treatment of interaction makes it possible to verify both the set of kinematic constraints and the space discretized mechanical system of equations. Afterwards it is possible to use an explicit or an implicit time integrator and a smooth or a nonsmooth description of interaction laws. To illustrate the different possibilities one can consider: • the approach developed by Carpenter et al. [7] based on an explicit time integration and an implicit treatment of contact. Contact conditions are applied either by Lagrange multiplier or by penalty method. • the so-called Granular Element Method, proposed by Kishino [18, 19] based on implicit time integration and an implicit treatment of smoothed contact. The solving method is similar to the penalisation techniques used in Finite Element Method. 6

Fr´ed´eric Dubois and Mathieu Renouf

• in a different way, Moreau developed the Contact Dynamics method [23] based on the convex analysis framework and dedicated to the simulation of rigid body systems. The time integration used is implicit as well as the contact resolution. Further works lead to the extension of the method to multi-contact simulations of collections of deformable bodies [15] and the method becomes the so-called Non Smooth Contact Dynamics method (NSCD). 2.2.2 Unconstrained dynamics When an unconstrained treatment (i.e. explicit) of contact is considered, time-steps are small enough to consider smooth dynamics. Typically the duration of a contact runs on several time steps. By the way, different explicit time integrator can be used to describe the motion of particle: time-centred scheme [10], Gear integrator [3], velocity Verlet scheme [14]. Smooth contact law describes the contact interactions as the Hertz contact model, the JKR model or more specific laws related to the medium properties. It is important to notice that such approaches are very simple to code but may be complicated to drive [4]. 2.2.3 Constrained dynamics Explicit Integration The most intuitive way to avoid the discontinuities due to contact-impact is to decrease the time step and to introduce smoothed interaction laws. Considering the time step, it leads, particularly when considering deformable components, to use explicit time integrator. Carpenter et al. [7] propose a method which respects the mechanical constraints (impenetrability) even using an explicit time integrator. The starting point of the method is the constrained equation of motion discretized in space (FEM, rigid,etc.). Considering the time increment [ti , ti+1 ] it leads to  ˙ i , qi ) = Fext M¨ qi + Fint + Hi+1 ri i (q i , (6) H∗i+1 qi+1 ≥ 0 ˙ i , qi ) represents the internal forces. As in the previous section, quantities indexed where Fint i (q by i (resp. i + 1) refer to time ti (resp. ti+1 ). In this approach, a second order scheme is taken into account, using the first and second time derivatives of the configuration parameter that appears to be the primary unknown of the problem as contact forces. Thus acceleration and velocity are related to configuration parameter using, for example, a β method,  1 2β   q˙ i = {q˙ i−1 + h(1 − β)¨ qi−1 + (qi+1 − qi )} 1 + 2β h , (7) 2   q ¨ i = 2 (qi+1 − qi − hq˙ i ) h where β (∈ [0.5, 1]) represents a numerical damping parameter. To take into account the displacement constraints in system (6), the displacement qi+1 is decomposed as ree qi+1 = qfi+1 + qci+1 ree where qfi+1 represents the prediction of the displacement without contact forces (similar to the free velocity of system (10)) and qci+1 the correction due to the contact forces only. In the case ree of the central different method (β = 0.5), qfi+1 is defined by ree ˙ i , qi )} + 2qi − qi−1 qfi+1 = h2 M−1 {Fext − Fint i i (q

7

(8)

Fr´ed´eric Dubois and Mathieu Renouf

To solve the frictional contact problem and to determine the solution of the equation of motion, the couple of unknown (ri , qci+1 ) solution of the following system must be found  2 ∗ ree h Hi+1 M−1 Hi+1 ri = H∗i+1 qfi+1 , (9) qci+1 = h2 M−1 Hi+1 ri Implicit Constrained nonsmooth mechanics A θ scheme is used for the time discretization of Equation (4). Scheme stability condition implies that θ remains between 1/2 and 1. The θ-method is an implicit scheme, equivalent to the backward Euler’s scheme when θ = 1. Proceeding to the time discretization, the contact problem is solved over the interval ]ti , ti+1 ]. Omitting internal forces this leads to  q˙ i+1 = q˙ fi ree + M−1 hRi+1 (10) qi+1 = qi + hθq˙ i+1 + h(1 − θ)q˙ i with ext q˙ fi ree = q˙ i + M−1 h(θFext i+1 + (1 − θ)Fi )

where q˙ f ree denotes the free velocity (velocity computed without contact forces). Quantities indexed by i (resp. i + 1) refer to time ti (resp. ti+1 ). For rigid body system, the inertia matrix M is diagonal and easily invertible, internal forces vanish and the external forces are given by a function of time. The extension to linear deformable bodies is straightforward. For nonlinear deformable bodies, a Newton scheme allows us to obtain an equivalent set of discretized equations [15]. The second Equation of system (10) is used to update the system. When θ = 1 the contact detection is performed using the configuration parameter at time ti . On the contrary, when θ 6= 1, the contact detection is performed using a prediction of the configuration parameter (at time ti + h(1 − θ)). In this case the θ-method corresponds to the well-known leapfrog scheme. 2.2.4 Contact problem resolution For unconstrained mechanical problems, no contact solver are necessary. As all contact are independents and smoothed, the computation of contact forces is explicit. For constrained mechanical problems, a contact solver must be selected. There exist a large amount of methods to solve multicontact problem (see [27] or refer to the NSSPACK library of the Siconos Project): Block Nonsmooth Gauss-Seidel [15], Conjugate Projected Gradient [28], Lemke [9], Generalised Newton method [2], etc. For the present applications, the solver used is the Block NonLinear Gauss-Seidel, robust and generic in regards of local interaction laws. The headlines of the solver are exposed below. Using the Equations (5) in the first Equation of system (10), the global discretization of the equation of motion and the contact laws can be summarised in the following system:  Whri+1 − vi+1 = −vf ree (11) ContactLaw[vi+1 , ri+1] where W (= H∗ M−1 H) is the Delassus operator, which models the local behaviour of the solids at the contact points. The right-hand-side of the first Equation in (11) represents the free relative 8

Fr´ed´eric Dubois and Mathieu Renouf

velocity only accounting for the internal and external forces F(t). The second Equation in (11) requires that the contact law must be satisfied by each component of the couple (vi+1 , ri+1 ). The global splitting scheme is written down as follows for the first Equation of system (10): P P Wαα rk+1 − vαk+1 = −vα,f ree − β<α Wαβ rk+1 − β>α Wαβ rkβ α β (12) = bkα where the index k refers to the splitting method iterations and nc the number of contacts. The time index is omitted to make pleasant reading. This solver has proved to be very robust and efficient on a large collection of heterogeneous problems [16, 33] and benefits of a parallel version [29] to ensure reduced simulation time. 3 PLATFORM ARCHITECTURE 3.1 Foreword The LMGC90 platform is a workbench devoted to the building of simulation tools relying on the DEM or hybrid FEM/DEM approach for the simulation of interacting physical systems (components): multibody contact dynamics, thermal convection, etc. In a sake of evolutivity the software has been designed using an object oriented approach. The global architecture, summarized by Fig. 3, aims to be very intuitive. It mimics the carving of the various stage of the problem (see the previous Section). In the following a general description will be given. Diving deeper in the architecture would be complex due to technical aspects, but interested people will find information on LMGC90 web site (http://www.lmgc.univmontp2.fr/∼dubois/LMGC90). Its global architecture takes benefits of published work on DEM/FEM codes; for examples [21, 36]. Parts of the global architecture of the platform are also quiet similar to existing tools, for examples SOFA (http://sofa-framework.org) and SICONOS (http://siconos.inrialpes.fr). Practically the code proposes a global pattern to model and solve the problem, with various functionnalities, which the user may enrich with its own routines through plug points. Basically the user needs to define the modelling ingredients (the bulk behaviour of the components of the system, the kind of interactions between these components, etc.) and the numerical strategy to simulate the evolution of the system. 3.2 Modelling framework The modelling framework provide low level functionalities necessary to describe the system and provide necessary functionalities to the strategy drivers. The architecture relies on the following class decomposition • Body. This abstract class encapsulate the dynamical system modelling. For mechanical components it concerns the bulk behaviour (rigid, deformable,etc.). One object of this class is represented by geometric information (nodes, lines, surfaces, volumes, meshes, etc.). Relying on these geometric informations, the bulk models are built. The algebraic description of the system relies on a System Of Equations class (SOE). Various mappings are also provided. This description can be linked to external dedicated software as PELICANS (https://gforge.irsn.fr/gf/project/pelicans) or Code Aster (http://www.codeaster.org). 9

Fr´ed´eric Dubois and Mathieu Renouf

Figure 3: Overview of LMGC90 class diagram

• Contactor. This abstract class encapsulates the interaction zone modelling. Contactor are relying on the body geometry. It manages essentially the mapping between body degrees of freedom and interaction degrees of freedom. It provides also some geometric information as local frame. • Interaction Handler This abstract class encapsulates the managment of interactions: detection (new/delete interaction) initialization and storage, etc. • Interaction This abstract class encapsulates the computation of interactions. It collects, in an “anonymous” manner, the interaction data. It provides a mapping through SOE. Various interactions are available. It manages the computation of interactions. For explicit interaction, it computes the interaction forces and adds it tho the SOE right hand side vector. For implicit solvers it build the reduced system and calls a dedicated solver class. Their are various classes managing specific operations : contact solvers, I/O, etc.

3.3 Numerical strategies framework The Numerical Strategies framework drives the simulation. Mainly one need to specify a time evolution strategy. For example when choosing a time stepping strategy one may use a linear contact dynamics which implies that the bulk problem will be supposed linear, the integrator will be a θ method, a Newmark scheme, etc. New time evolution strategies are easy to derive due to the high-level modularity of the platform. However, it specifies some low level functionalities: SOE storage, contact solver, etc.

10

Fr´ed´eric Dubois and Mathieu Renouf

4 MULTIPHYSICS FORMULATION The last, and still unstable, evolutions of LMGC90 platform concern the multiphysics models. 4.1 Electrical When the electrical feature of granular assemblies under dynamical solicitations is investigated (powder, steel bead assemblies...), it is difficult to determine which effects are purely electrical than the ones that are strongly dependent of mechanical effects. Although numerous experimental works have allowed stepping forward in the understanding of the static electrical behaviour of such media, it becomes difficult to determine which effects control the multi-physical behaviour of the medium, especially under dynamical solicitations. Thus an electrical-mechanical model embedded in discrete element schemes and able to cope with the oxidation phenomenon has been developed and implemented in the platform. The formulation of the electrical problem relies upon an analogy between an electrical and a contact network, each particle of the granular assembly is considered as a node of the electrical network while each contact is considered as a branch of it. A nonlinear electrical problem is built using classical electrical laws (Kirshoff and Ohm laws), the definition of local conductances and some hypotheses to take into account the oxidation phenomenon (see [31] for more details). Considering V as the electrical potential expressed at the centre of particles and I0 a vector related to the current intensity introduced in the system, the problem is written as: GV = −I0 ,

(13)

where the matrix G is equal to the product NCNT , with C the conductance matrix and N and NT two linear mapping with a role equivalent to the one of H mapping defined by System (5). When oxidation is considered, the problem consists in solving a succession of linear problem (13) with a matrix G which take into account oxidation transformation. Actually, the electrical problem at time ti+1 is built using the mechanical solution at time ti+1 and the electrical solution at time ti . The structure of problem (13) does not allowed to solve simultaneously the mechanical problem and the electrical problem. The first one is formulated at the contact level while the second one is formulated at the body level. Moreover, on one time step, there is no fixed point iteration between the two problems: The solution of the electrical problem is not used to update the parameter of the mechanical problem (adhesion for example). Indeed, as time steps are small (˜10−6 ), a strong coupling is not necessary and all linear and nonlinear characteristics of the model have been validated by experimental comparisons. 4.2 And what else ? Thermal effects are also considered. A classical thermo-mechanical model can be used to model deformable bodies [35]. Its extension to rigid bodies as been performed. An application is given in the next Section. In order to model fluid-structure interactions a meshless approach, quiet similar to discrete element, will be used. 5 APPLICATIONS The panel of applications chosen in the present paper aims to illustrate how the number of complex systems involved in this extended framework is large: geophysics, Tribology, coupled11

Fr´ed´eric Dubois and Mathieu Renouf

behaviour problem are especially concerned. 5.1 Geophysics Many geophysical problems bring difficulties such as large deformations and multiple localisations. If for some applications efficient numerical methods are well identified, there exist some problems where no well-suited methods have been explicitly determined. It is the case of two phenomena: Propagation fault and forced-fold evolution. The study of folds relies upon post-mortem observations on geological structures. With numerical tools, it becomes possible to reproduce them and to do correlation between their formations and the mechanical properties of the geological structure. The reader can refer to [30] for a complete study including numerical observations and mechanical analysis. (a)

(b)

Figure 4: Forced-fold evolution: a) The initial structure and b) the structure at the end of the simulation

The typical sample used for forced-fold evolution simulations is illustrated in its initial configuration by Fig. 4a. It is composed of 43 000 hard cylinders. Their radius range from 0.42 to 0.56 m. Their density is equal to 2 800 kg.m−3 . The contact law is a frictional contact law (Signorini-Coulomb law with µ = 0.4) with inelastic quasi-shocks [15]. All coloured layers have the same mechanical properties. The simulation consists in applying a constant velocity on the left wall during the whole process. The simulation is decomposed in 20 000 time steps of 2.5 10−2 s. For the NLGS algorithm, a minimal number of iterations equal to 100 is imposed to ensure the quality of the solution. The elapsed time of the simulation on a single process (3.6 GHz frequency) is equal to 96 h approximatively. The Figure 4b represents the deformation of the structure at the end of the simulation. Internal faults and avalanches are generated by the structure evolution. Faults generate vertical motions with few incidences on the evolution of the front while avalanches affect the position of the front in the direction of wall motion. 5.2 Wear modelling During wheel-rail interaction (strong pressure and high local shear velocity), as many contact problems, the contact between wheel and rail, called first bodies, generates a natural third-body composed of particles detached from the two first bodies. Thus the created body plays the important role of contact interface between the two first bodies and it is the clue to understand the thin link between friction and wear [5]. The sample illustrated by Fig. 5a represents a sample used to study lateral third-body flow during a wheel-rail interaction under constant loading and cycling displacements [32]. The portion of rail is filled with 73 000 disks which radius ranges from 5.10−2 to 10−1 mm. The 12

Fr´ed´eric Dubois and Mathieu Renouf

F (a)

V

1111111111111111111111111111111 0000000000000000000000000000000 0000000000000000000000000000000 1111111111111111111111111111111 0000000000000000000000000000000 1111111111111111111111111111111

(b)

particle without cohesive interactions particle without damaged interactions

Figure 5: Fretting-like simulation - a) initial sample composed of 73 000 disks and b) zoom on the interface to observe internal particle flows

rigid body, which represents the wheel, is composed of 450 disks. Their radii have the same property than the radii of disks located on the bulk. A cohesive zone model law governs particle interactions. By this way, the whole medium has a continuous characteristic and is able to be strongly wear away locally. The wheel exerts on the rail a constant force and is submit to lateral displacements. These contact conditions reproduce the fretting-like motion which occur during a real wheel-rail interaction. The sample illustrated by Fig. 5b is a zoom of the contact interface after approximatively one hundred cycles. Red particles have undamaged interactions while blue particles have any undamaged interactions: All interactions have been broken and the particles can flow easily inside the interface. The generated interface to reach this configuration, the simulation run during 30 days on a single process (3.6GHz frequency). 5.3 Fracture In experiments of Kalthoff and Winkler [17], a plate (L = 100mm) with two edge notches is impacted by a projectile as shown in figure 6a. A symmetry condition is taken into account to model half of the structure. The projectile impact effect is applied using a step velocity on the cracked edge. The impact velocity is chosen as V = 16.5m/s. The material properties are: ρ = 8000kg/m, E = 190GP a, ν = 0.30. A 100 x 100 cross-triangular mesh is used for the space discretization. The Non Smooth Fracture Dynamics (NSFD) approach is used for the simulation, which is based on the three features: a multibody description with mixed boundary conditions between each mechanical component, an enhanced cohesive zone model and a specific non smooth dynamical framework (see [1] for details). One can observe, at this low strain rate, a brittle fracture mode with a propagation angle of ∼ 70◦ (figure 6b). It is coherent with the experimental results. 5.4 Metallic powder coupled behaviour Among most interesting multiphysical problems, the electrical transport in metallic powder is one of them. This interest dates more than one century with the works of Branly (1873) and is always of the moment: More recent researches have been carry on to fill up the lack of 13

Fr´ed´eric Dubois and Mathieu Renouf L

(a)

(b)

L/2 V

2L

L/2

Figure 6: a) Experimental set-up for crack growth in an edge-cracked plate and b) crack pattern with a 100 x 100 cross-triangular mesh at t = 29.µs

knowledge on the electrical behaviour of metallic powder and granular packing [13]. Moreover the solution of problem such as the electrical-mechanical behaviour of the wheel-rail contact (mechanical shunt), the development of synthetic diamond crystals or the cold compaction of composite powders depends on progresses in this research field to make a step forward. Thus, as presented in the section 4.1, an approach coupling electrical and mechanical behaviour has been developed [31] and validate using comparison with experiments. (a)

(b)

F + −

Figure 7: Cohesive sample submit to an electrical current: a) simulation procedure and b) visualisation of the electrical potential expressed on the node of each particle; the oxide layer present in the sample created an electrical heterogeneity in the sample. Walls have been removed to observe particles freely.

The Figure 7 represents a sample used to validate the mechanical-electrical algorithm: Walls have been removed to observe particles freely. 15 000 spheres are compacted between two parallel plates in a cylinder of diameter equal to 0.09m. Particle diameters are equal to 2 10−3 m (±0.2 10−3 m). The sample is crossed by an alternative electrical current which varying from 0 to 0.2 A. When the current increase, breakdown of oxide layer occur in the sample and its resistivity decreases. 5.5 Masonry Masonries are heterogeneous material composed of blocks, with various shapes and sizes, separated by joints. Numerous models have been proposed in the literature in order to study the 14

Fr´ed´eric Dubois and Mathieu Renouf

mechanical behaviour of masonry. When study concerns the built heritage (stone bridge, cathedral, etc.), robust DEM are necessary to represent properly the whole structure [8]. When the study concerns complex applications, it could be necessary to develop new model to obtain a solution. 5.5.1 Monument

(a)

(b)

Figure 8: The Pont du Gard structure: Discrete Element assemblies composed of 35 000 rigid regular polyedra.

The global structure of the Pont du Gard is composed of 35 000 rigid regular polyedra which interact by a Signorini-Coulomb law. The friction coefficient between each element is equal to 0.5. The simulation running on a PC takes one day. The simulation have been performed to study the internal load equilibrium. 5.5.2 Industrial structures This application concerns the flue walls of industrial baking furnace. These walls are made of firebricks linked together by mortar in horizontal joints. During service, the flue walls are submitted to thermal cyclic loading (see Fig. 9). Complex deflection of the structure can be observed due to the thermal field expansion effect, which are restrained, in their free thermal expansion by the action of head walls. In this study, we introduce the thermal expansion effects on the 3D rigid bricks through an explicit computation of the “contactor” shape, we take into account the mortar damage using a cohesive zone model [16] in order to obtain the irreversible deformation of the wall due to cyclic thermal loading. 5.6 FEM-DEM coupling In tribology, as in many contact applications, it becomes necessary to take into account the behaviour of the structure as well as the one of the interface. The rheology of the interface affects the behaviour of bodies in contact but bodies in contact affect also the rheology of the interface: They are both the limit conditions of the other. As a consequence, to face this problem, it is necessary to simulate the discontinuous interface and the continuous bodies. The Figure 10 illustrates the FEM-DEM coupling for the compaction of a heterogeneous medium.

15

Fr´ed´eric Dubois and Mathieu Renouf

Figure 9: Industrial masonry structure submitted to thermal cycles. Left: snapshot of the thermal field applied on the structure. Right: cumulated horizontal displacements after some thermal cycles.

Figure 10: Zoom of the FEM-DEM simulation used to study the influence of interface heterogeneities on the behaviour of bodies in contact. Visualisation of equivalent constraint of particles and stress field of elastic bodies

The cohesive medium is composed of 3 000 rigid disks. The two plates are two elastic bodies. The size of an element of the mesh is equal to the smaller particle diameter. The law used between disks and mesh is the same that the one between disks: A nonsmooth Lehnard jones law, translation of the Signorini graph. 6 CONCLUSION A general framework, relevant for a class of multiproblems (multibodies, multicontacts, multiphysics and multiscales), has been proposed. It is shown how existing strategies may be described in this framework. These strategies are mainly different due to constrained/unconstrained modelling of contact and implicit/explicit time integrators. Moreover the interest of such unified modelling framework has been emphasized with multiphysics modelling as electrical-mechanical coupling. Various applications demonstrate the interest of such unified modelling framework in mechanics as well as in geophysics and tribology.

16

Fr´ed´eric Dubois and Mathieu Renouf

7 ACKNOWLEDGEMENTS All Simulations have been performed with the OpenSource platform LMGC90 (http://www.lmgc.univmontp2.fr/∼dubois/LMGC90). This work is supported by the Centre Informatique National de l’enseignement Sup´erieur (CINES) under projects mgc2547 and lmc2644. REFERENCES [1] V. Acary and Y. Monerie. Nonsmooth fracture dynamics using a cohesive zone approach. Technical report, INRIA - RR6032, 2006. [2] P. Alart and A. Curnier. A mixed formulation for frictional contact problems prone to Newton like solution method. Comput. Methods Appl. Mech. Engrg., 92(3):353–375, 1991. [3] M.P. Allen and D.J. Tildesley. Computer simulation of liquids. Oxford University Press, 1987. [4] T. Belytschko, W.K. Liu, and Moran B. Nonlinear Finite Elements for Continua and Structures. John Wiley & Sons, 2000. [5] Y. Berthier. Wear - Materials, Mechanisms and Practice, chapter Third-body reality Consequences and use of the third-body concept to solve friction and wear problems, pages 291–316. John Wiley & Sons, 2006. [6] B. Cambou and M. Jean. Microm´ecanique des mat´eriaux granulaires. Herm`es ScienceParis, 2001. [7] N. J. Carpenter, R. L. Taylor, and M. G. Katona. Lagrange constraints for transient finite element surface contact. Int. J. Numer. Methods Engrg., 32:103–128, 1991. [8] B. Chetouane, F. Dubois, M. Vinches, and C. Bohatier. Nscd discrete element method for modelling masonry structures. Int. J. Numer. Meth. Engrg, 64(1):65–94, 2005. [9] R. W. Cottle, J. Pang, and R. E. Stone. The linear complementarity problem. Academic Press, Inc., Boston, MA, 1992. [10] P. A. Cundall. A computer model for simulating progressive large scale movements of blocky rock systems. In Proceedings of the symposium of the international society of rock mechanics, volume 1, pages 132–150, 1971. [11] P. A. Cundall and O. D. L. Strack. A discrete numerical model for granular assemblies. Geotechnique, 29(1):47–65, 1979. [12] C. Duriez, F. Dubois, A. Kheddar, and C. Andriot. Realistic haptic rendering of interacting deformable objects in virtual environments. IEEE Transactions on Visualization and Computer Graphics, 12:36–47, 2006. [13] E. Falcon, B. Castaing, and C. Laroche. ”Turbulent” electrical trasnport in copper powders. Europhys. Lett., 65(2):186–192, 2004.

17

Fr´ed´eric Dubois and Mathieu Renouf

[14] I. Iordanoff, B. S`eve, and Y. Berthier. Solid third body analysis using a discrete approach: Influence of adhesion and particle size on the macroscopic behavior of the contact. ASME J. Tribol., 124:530–538, 2002. [15] M. Jean. The non smooth contact dynamics method. Compt. Methods Appl. Math. Engrg., 177:235–257, 1999. [16] M. Jean, V. acary, and Y. Monerie. Non-smooth contact dynamics approach of cohesive materials. Phil. Trans. R. Soc. Lond. A, 359:2497–2518, 2001. [17] JF Kalthoff and S. Winkler. Failure mode transition at high rates of shear loading. In International conference on Impact Loading and Dynamical Behavior of Materials, volume 1, pages 185–195, 1997. [18] Y. Kishino. Disk model analysis of granular media. In Micromecanics of Granular Materials, pages 143–152, 1988. (eds. M. Satake and J.T. Jenkins),Elsevier. [19] Y. Kishino. Granular flow simulation by granular element method. In Slow Dynamics in Complex Systems, November 2003. Sendai, Japan. [20] S. Li and W.K. Liu. Meshfree Particle Methods. Springer, 2007. [21] F.T. McKenna. Object-Oriented Finite Element Programming: Frameworks for Analysis, Algorithms and Parallel Computing. PhD thesis, University of California, Berkeley, 1997. [22] J. J. Moreau. On unilateral constraints, friction and plasticity. In G. Capriz and eds. G. Stampacchia, editors, New Variational Techniques in Mathematical Physics (C.I.M.E. II ciclo 1973, Edizioni Cremonese, Roma, 1974), pages 173–322, 1974. [23] J. J. Moreau. Unilateral contact and dry friction in finite freedom dynamics. In J.J. Moreau and eds. P.-D. Panagiotopoulos, editors, Non Smooth Mechanics and Applications, CISM Courses and Lectures, volume 302 (Springer-Verlag, Wien, New York), pages 1–82, 1988. [24] J. J. Moreau. Numerical aspect of sweeping process. Comptut. Methods Appl. Math. Engrg., 177:329–349, 1999. [25] A. Munjiza. The Combined Finite-Discrete Element Method. John Wiley & Sons, 2004. [26] F. Pfeiffer and C Glocker. Multibody dynamics with unilateral contacts. Non-linear Dynamics. John Wiley and Sons, 1996. [27] M. Renouf and V. Acary. Comparisons and coupling of algorithms for collisions, contact and friction in rigid multibody simulations. In C.A. Mota Soares et.al. (eds.), editor, III European Conference on Computational Mechanics Solids Proceedings, Lisbon, Portugal, 58 June 2006. [28] M. Renouf and P. Alart. Solveurs parall`eles pour la simulation de syst`emes multi-contacts. Rev. Eur. Elem. Finis, 2004. [29] M. Renouf, F. Dubois, and P. Alart. A parallel version of the Non Smooth Contact Dynamics algorithm applied to the simulation of granular media. J. Comput. Appl. Math., 168:375–38, 2004. 18

Fr´ed´eric Dubois and Mathieu Renouf

[30] M. Renouf, F. Dubois, and P. Alart. Numerical investigations of fault propagation and forced-fold using a non smooth discrete element method. Rev. Euro. Method. Num., 15:549–570, 2006. [31] M. Renouf and N. Fillot. Coupling electrical and mechanical effects in discrete element simulations. Int. J. Numer. Meth. Engng. submited january 2007. [32] M. Renouf, A. Saulot, and Y. Berthier. Third body flow during a wheel-rail interaction. In C. A. Mota Soares et.al., editor, III European Conference on Computational Mechanics Solids, 58 June 2006. [33] G. Saussine, C. Cholet, F. Dubois, C. Bohatier, P.E. Gautier, and J.J. Moreau. Modelling ballast behaviour under dynamic loading. part 1: A 2d polygonal discrete element method approach. Comput. Methods Appl. Mech. Engrg., 195(19-22):2841–2859, April 2006. [34] L. E. Sibert, D. Ertas, G. S. Grest, T. C. Halsey, and D. Levine. Analogies between granular jamming and the liquid-glass transition. Phys. Rev. E, 65(5):1307–1310, 2002. [35] L. Stainier, F. Dubois, and R. Peyroux. Une bibliothque portable de modles constitutifs pour la mcanique non-linaire des solides: concepts et implmentation. In Actes du sixi`eme colloque national en calcul des structures, volume 1, pages 111–118. CSMA-AFM-LMS, 2003. [36] C.-T. Yang and S.-H. Hsieh. An object-oriented framework for versatile discrete objects simulation using design patterns. Computational Mechanics, 36(2):85–99, 2005.

19