A NonSmooth Contact Dynamics-based multi-domain solver

This paper presents a micromechanical modeling strategy for complex multibody interactions and the associated numerical framework. The strategy rests on a periodic multibody method in the framework of the NonSmooth Contact Dynamics approach of Moreau (1988) extended to classical domain decomposition problems. Many complex interactions can be taken into account: interactions between discrete elements, between discrete or rigid bodies, (quasistatic) contact or impact, friction or adhesion, decohesion (cracking), etc. The associated numerical platform, Xper, is composed of three independent libraries with Object Oriented Programming. The ability of this computational approach is illustrated by two examples of fracture in heterogeneous materials.


Introduction
Dynamic crack propagation in heterogeneous materials is a complex problem and has become a challenge in many engineering domains. Since the heterogeneities (microcavities, hard inclusions, brittle precipitates, ...) have major effects on the overall dynamic fracture, a relevant way to investigate the problem is to consider models at the scale of the heterogeneities. At this scale, the fracture can be described using micromechanical concepts, and recent advances in computer technologies make possible the simulation of the overall nonlinear response of a given microstructure.
In this context, a new computational micromechanical approach is developed to analyse the effects of the microstructure heterogeneity on the overall material behavior submitted to static or transient loadings. This approach is based both on the concept of Frictional Cohesive Zone Model and on a multibody method in the context of the NonSmooth Contact Dynamics (NSCD). In particular, the NSCD approach aims to solve dynamic frictional contact problems without regularization nor penalization techniques (Moreau, 1988;Jean, 1999). Since periodic formulations are well adapted to micromechanical studies, a two field Finite Element is written and the framework is extended to this formulation (Perales et al., 2006;Perales et al., 2008).
The scope of the NSCD framework can simply be extended to classical domain decomposition problems. The domain decomposition methods (DDM) solve a standard boundary value problem by splitting it in smaller problems on subdomains and managing "continuity" conditions between subdomains. Many techniques exist to enforce continuity of the solution in quasi-static (Dodds et al., 1980;Magoules et al., 2006;Glowinski et al., 1990) or dynamic (Herry et al., 2002;Gravouil et al., 2001). In this paper, the NSCD framework is written as a dual Schur approach where the continuity condition is replaced by any interaction law. The applications of this method concern parallelism, multiphysics or multibody interactions.
Since the implementation of the micromechanical framework would involve a high programming cost, the development strategy is to reuse and to extend existing specialized libraries. The software developed here, called Xper, is based on the coupling of three libraries (using Fortran90/C++ mixed programming) and thus takes advantage of the ability of each of them. Each library has a clear meaning from the mechanical point of view : -LMGC90 is dedicated to the surfacic interaction part : model, solving method. It relies on the NonSmooth Contact Dynamics approach (LMGC90, 2009), -PELICANS is dedicated to the bulk part : periodic Finite Element modeling (PELICANS, 2009), -MatLib is dedicated to complex constitutive models (Stainier et al., 2003). It is embeded in PELICANS.
The ability of the proposed strategy and of the platform is illustrated on the dynamic fracture of metal matrix composite with brittle inclusions under transient loading and on intergranular fracture in a periodic medium.
The paper is organized as follows. In section 2, the main elementary classes of multibody interactions that can summarize most of the complex interaction problems are recalled. In section 3, the interactions modeling (multibody, NSCD, periodic and DDM) is described. The aforementioned case study are presented in section 4.

The modeling strategy
The main purpose of this paper is to model and simulate divided and fractured media as a multibody system with interactions. From a mechanical point of view, these interactions correspond to any compression-sliding-traction situation. In particular, we have here in mind to model contact, impacts, dry or lubrificated friction, adhesion, decohesion, multiple cracking, failure or any combinations as frictional contact with cohesion (Raous et al., 1999;Perales et al., 2006). Moreover these interactions can take place between rigid or deformable media. From a geometrical point of view, any complex interaction belongs to one of the three following situations : (1) the body-to-body interaction, (2) the cluster-to-cluster interaction, (3) the body-to-cluster interaction. In order to classify these three situations, the two first are detailed below.
-Body-to-body interaction. For example, this situation occurs when any surface of a Finite Element mesh is susceptible to surface interaction with its surrounding (see Figure 1). This class of interactions includes the failure of bulk materials or of heterogeneous materials, or transgranular failure. In this kind of applications, the framework is known as the cohesive/volumetric finite element approach (Xu et al., 1994;Camacho et al., 1996;Raous et al., 1999;Jean et al., 2001;Perales et al., 2006). This concept is extended to any interaction of rigid bodies. Therefore, the body-to-body interaction stands for mesh-to-mesh interaction, rigid-to-rigid interaction or mesh-torigid interaction, 'Independent' body Surface behavior Figure 1. body-to-body interaction : each finite element is a body.

Meshed body
Surface behavior For the particular case of periodic media, the clusters can be subject to periodic conditions (see Figure 3). This class includes the failure of periodic media or of Representative Volume Elements (Perales et al., 2008;Pelissou et al., 2009).
The overall behavior is then obtained by coupling the standard volumetric behavior inside the meshes (or rigid bodies) and the surface properties between the bodies, taking complex interactions into account. In this study, this coupling is based on the NonSmooth Contact Dynamics (NSCD) framework (Moreau, 1988;Jean, 1999).
In the framework of the NSCD method, a two level resolution is carried out : the standard or periodic volumetric problem is solved at global level -by PELICANS (PELICANS, 2009) and MatLib (Stainier et al., 2003) -and the non smooth contact problem is treated at a local level -by LMGC90 (LMGC90, 2009) -.
The corresponding two field modeling framework is detailed next.

Global level and periodic problem
A specific bulk model is derived to model a dynamical system under periodic conditions. Since the classical dynamic problem can be seen as a subproblem of the periodic one, we focus on the formulation of the periodic problem. In what follows, the main points of the formulation are described (see (Perales et al., 2008) for more details).
Consider a periodic multibody medium Ω 0 = e Ω e 0 . At any boundary of a body Ω e 0 , mixed boundary conditions (given by the interactions depending on the displacement jump) are introduced. In this framework, the deformation gradient field F = ∇u + I and the first Piola-Kirchhoff stress field Π are assumed to be periodic with the same period as the medium. The average fields over the periodic medium are denoted byF,∇u andΠ. The fields F, ∇u and Π fluctuate around their average values. The local deformation gradient field can be thus split into an overall field (the field if the medium were homogeneous) and a fluctuation denoted ∇u # , which takes the presence of heterogeneities into account. Since fracture is expected, the heterogeneities are not only due to the inclusions but also the cracks in the structure. The global displacement field u admits the following decomposition : u = (F − I) · X + u # where X is the initial position vector and I is the second-order identity tensor. The Finite Element formulation becomes a two field Finite Element formulation. The local periodic dynamic problem can be written : find the periodic displacement field u # , the deformation gradient field F and the stress field Π verifying : -relations for each body Ω e 0 : -average relations in the medium Ω 0 : where T is a boundary force given by the considered interaction law, N is the unit outward normal vector of the body, ρ is the density, the jump symbol [f ] = f + − f − is defined as the difference of a field f on the two facing surfaces (here the superscripts + and -denote the two opposite surfaces),F imp ij are the components of the prescribed macroscopic transformation gradient andΠ imp kl are the components of the prescribed macroscopic stress, with ij = kl and {i, j, k, l} ∈ {1, 2, 3}. Note that the classical dynamic problem can be obtained by replacing the periodic displacement field u # by the standard displacement field u in the first equation of [1].
The framework is dedicated to the study of periodic problems embedding non regular interactions. Considering the standard NSCD algorithm, the non smooth contact problem is treated at the local level (Jean, 1999;Jean et al., 2001). At this level, two main points have to be underlined in the writing of the periodic problem : -dynamics are dedicated to the treatment of the non regular conditions and has to be only introduced at the local level, -the presence of the heterogeneities are taken into account only by the fluctuation field u # .
After choosing the admissible spaces, the broken Sobolev space m ∀Ω e 0 , v periodic and the space of linear transformation V = L(R m ) (m is the space dimension) for the velocity field and its periodic part respectively, and considering the kinematically admissible virtual velocity field v * = (v * ) # +Ḟ * · X, the weak unit cell value problem is then obtained : [3]

Standard dynamic equation and the extension to periodic
Considering that some discontinuities may appear in velocity time evolution, the standard dynamic problems are written in a semi-discrete form for each body Ω e 0 : where M is the mass matrix, q,q and dq are respectively the discrete displacement, velocity and differential measure of velocity, dp represents the differential measure of interaction impulse and F (q,q, t) represents the internal and external forces without the contribution of interaction dp (Jean, 1999;Jean et al., 2001). The measure dp may contain both smooth and non-smooth contributions. In this framework, the derivatives are written in a distribution sense.
Practically the differential measure equation [4] needs to be integrated over an arbitrary time step (even reduced to a shock instant), which gives a momentum balance : The linearized formulation of [8] is obtained through a Newton-Raphson method (superscript k stands for iterations) (Jean, 1999;Jean et al., 2001) : [10] which leads to : where w k i+1 denotes the inverse of the iteration matrix, h is the length of time subinterval ]t i , t i + 1], subscript i the quantities at a time t i and i + 1 the quantities at time t i+1 , θ ∈ [0.5, 1] and p k+1 i+1 is the impulse value. In a periodic case (see section 3.1 for the formulation), the semi-discrete form [4]-[5] becomes for each body Ω e 0 (Perales et al., 2008) : where G(q # ,q # ,d,ḋ, t) and K(t) represents respectively the macroscopic stress and the macroscopic prescribed stress, andd andḋ are respectively the discrete average deformation gradient and the first time derivative.
The linearized formulation [10] is then rewritten as : The unknowns of the periodic problem areq # andḋ. A mapping P permits to recover the standard discrete velocityq fromq # andḋ : with P such that : where Rel ∈ R N dof u × R m is a discrete mapping such that the periodic velocity is is the discrete average mapping over Ω 0 , m is the space dimension and N dof u and N dof F are the number of degrees of freedom of the discrete velocityq and of the first derivative average deformation gradientḋ respectively.
At the iteration k + 1, the discrete velocity and the discrete free velocity write : [17] The local problem at the contact level is thus solved using the standard NSCD algorithm without any modification. The two field periodic formulation can thus be seen as a simple extension at the global level of the NSCD algorithm, that is to say the resolution of the periodic Finite Element problem.

Kinematic relations for frictional contact
Writing interaction behavior in a discrete form requires to define mappings H between bulk unknowns (q, p) and interaction unknowns (see Figure 4) : relative velocity (U) and impulse (R). Considering one interaction (index α) between two bodies (c and a) and using classical kinematic relations, the relative velocity writes : Due to duality consideration, one may write the contribution of the interaction α to the global impulse as : Using the linear mappings H, the dynamic system of equations is rewritten as : One can generalize the previous form in case of a multi-interaction situation, simply by including an additional contribution : To close the problem, one needs an interaction law relating relative velocity and impulse. Numerous choices are possible depending on the phenomenology of the interaction , for example, unilateral condition (velocity Signorini condition (Moreau, 1988)), friction (Coulomb's law) or cohesion (Jean et al., 2001;Perales et al., 2006) .
As detailed further, the general traction-sliding-compression interactions presented in the 'modeling strategy' (section 2) have to include sliding-compression relations which are causing impulse conditions. For sake of clarity, unilateral frictional contact laws of Signorini-Coulomb type are presented but other choices can be made (e.g. Tresca friction law). These relations write : where I K is the indicator function of the set K, μ is the Coulomb friction coefficient and U is decomposed such as U = U N n + U T with n the unit normal vector of the frictional contact zone.
The principle of the global to local mapping and the associated variables are summarized in a standard case and in a periodic case respectively in Figure 4 and in Figure 5. In particular, Figure 5 shows that the periodic mapping concerns only the global level of the algorithm and that the kinematic relations and local NSCD solver are not affected.q

NSCD as a Domain Decomposition Method
In classical domain decomposition methods (DDM) one assumes contiguous subdomains and tries to enforce continuity of primal (displacement) and/or dual (force) interface unknows. Considering quasi-static problems, there exist various techniques to ensure this continuity which differ in the treatment of the interface between subdomains. In the following, we consider techniques without overlapping between subdomains. Among them one can consider primal Schur-type approaches (where continuity of the displacement field is enforced at the interfaces) (Dodds et al., 1980), dual Schur-type approaches (where equilibrium of interface forces is enforced through Lagrange multiplier) (Magoules et al., 2006) and mixed approaches where one tries to achieve both conditions (Glowinski et al., 1990). Considering dynamical problems one can potentially use the three previous approaches. But as mentioned by Gravouil et al. (Herry et al., 2002;Gravouil et al., 2001), in dynamics, the choice of which kinematic quantity is continuous at the interface (displacement, velocity or acceleration) is difficult. In the case of contact problems various possibilities exist. Dureissex et al. (Dureisseix et al., 2001) propose an extension of a dual approach (FETI). Alart et al. (Alart et al., 2003) propose an approach where contact is embedded in a continuous domain. More recently, in Nineb (Nineb et al., 2007) and Icéta (Iceta et al., 2009), domain decomposition approaches are proposed for non smooth problems, where the continuity condition is enforced through the bulk model.
Our purpose in this part is to show that the NSCD method may be written as a dual-Schur approach where the basic continuity condition may be replaced by any interaction law. In that case the continuity condition between subdomains is not written on the bulk part but through an interaction law.
For the sake of simplicity we will only consider in the following two domains in interaction. In the context of the dual Schur formulation one may rewrite the standard discrete problem ([10] and [18]) involving two subdomains (a and c) : where the basic continuity condition of DDM written in velocity is replaced by a more general implicit interaction law (law(U k+1 , R k+1 ) = true). The free motion is computed at the global level and is independent of interaction conditions. It can be computed by any classical FEM library such as PELICANS. 2 The interaction motion can be computed at the local level solving two subproblems (omitting k+1 ) : The local solver needs the projection of the free velocity on the boundary of the subdomains and the projection of the inverse of the pseudo mass matrixM a orM c .

Multi domain solver
The classical NSCD approach, and its periodic extension, rely on a Non Linear Gauss Seidel (NLGS) algorithm to solve the problem. The spirit of the method is the following. Considering, one by one, the local systems to solve for each contact α : the contributions due to other contact (β = α) are frozen taking updated values if β < α or old values if β > α.

Global-to-local strategy and dedicated platform
The corresponding global-to-local strategy is summarized for one iteration of the Newton-Raphson algorithm in Figure 6.
The strategy rests on the local/global levels : -at the global level, the Finite Element method (including the periodic extension) is taken into account, -at the local level, the standard non smooth contact framework is used.
For lower programming cost, the software strategy development retained here was not to develop the entire software from scratch but to reuse and extend existing libraries. It permits to take advantages of each library update while developing the whole software.
The software platform, called Xper ('eXtended cohesive zone model and PERiodic homogenization'), respects the two levels strategy (Figure 6) : -at the global level, the extended Finite Element method is managed by the PELI-CANS library (PELICANS, 2009). This software, developed by the French 'Institut de Radioprotection et de Sûreté Nucléaire' (IRSN) is a toolbox for the implementation of various numerical methods dedicated to the solution of systems of partial differential equations (PDEs). To take into account complex non linear constitutive models, the library is coupled with the MatLib library, developed by Stainier (Stainier et al., 2003) al., 2003). LMGC90 is a platform for the modeling interaction problems including multi-physics.
For more details concerning the implementation of the software, see the appendix A.
In the following, both the abilities of the numerical strategy and the Xper software are illustrated on fracture of heterogeneous media.

Cohesive zone interaction law
Examples presented deal with fracture of heterogeneous materials. The fracture model is based on a cohesive-volumetric micromechanical approach involving Frictional Cohesive Zone Models (FCZM). These models rest on coupling tractionseparation interaction law to some frictional contact model. The cohesive-friction coupling is a key concept of fracture of heterogeneous media : whatever the overall loading of such a media (including pure mode I loadings), interfaces between phases with different Poisson ratio are locally subjected to combination of shear and tensile loading.
The interaction law relating the displacement jump U to the stress R, used to close the NSCD framework, is obtained introducing a cohesive stress, called R coh , in the Signorini-Coulomb problem [22]-[23] (Perales et al., 2006) : [34] where n is the unit normal vector of the cohesive zone, I K is the indicator function of the set K, μ is the Coulomb friction coefficient, C N and C T denote respectively the initial normal and tangential stiffness of the perfect interface (M P a/m).
The surface variable β, initially introduced by Fremond (Fremond, 1987), plays the role of a surface damage variable. The evolution law of this variable is governed by eqs. [36] and [37], where the function g describes the weakening process leading from perfect interface to crack (β = 1 : the interface is undamaged, 0 < β < 1 : the interface is partially damaged and β = 0 : the interface is fully damaged) : where δ 0 = R max 2 , 0 ≤ β 0 ≤ 1 is an initial surface damage, w is a reference fracture energy (J/m 2 ), R max is the maximum value of the cohesive stress (M P a), U max is the maximum value reached by U during the fracture process. In a 2D case, Figure 7 shows respectively (a) the normal behavior (with U T = 0) and (b) the tangential behavior (with R N constant) associated with the FCZM [33]- [37]. This Frictional Cohesive Zone Model takes into account the progressive damage between two bodies and the post-fracture frictional contact on the created crack lips. Different models can be used according to the material (brittle, ductile) by specifying the different evolution laws [36]-[37] (Michel et al., 1994;Tvergaard, 1990;Alfano et al., 2001;Perales et al., 2006) .

Examples
In what follows, the finite element discretization is based on linear displacement triangular elements that are arranged in a "crossed-triangle" quadrilateral pattern. The analysis considers 2D plane-strain conditions. The considered miscrostructure is composed of a metal matrix (Zircaloy-4) and inclusions (δ-hydrides). The Zircaloy-4 behavior is assumed to be elastoplastic (J2 plasticity, Young Modulus E = 99GP a, Poisson's ratio ν = 0.325, Yield in tension σ 0 = 450M P a, Hardening Modulus H Y = 850M P a) (Balourdet et al., 1999;Cazalis et al., 2005) and hydrides to be neo-Hookean (E = 135GP a, ν = 0.32) (Yamanaka et al., 1999;Yamanaka et al., 2001). We assume that the Zircaloy-4 and the hydrides inclusions have the same density ρ = 7800kg/m 3 . The FCZM coefficients of Zircaloy-4, zirconium hydrides and Zircaloy-hydride interface are respectively : C Zr N = 2 × 10 18 P a/m, w Zr = 0.5J/m 2 , R Zr max = 1GP a, C ZrH if the interface is considered as 'weak' or w Zr-ZrH = 1000w Zr , R Zr-ZrH max = 45R Zr max if the interface is considered as 'strong'. Moreover, we assume a low friction coefficient μ = 0.05 and same compliance for the normal and tangential behaviors C N = C T .
The PELICANS library manages the finite element part of the problem : the structure geometry, the inclusions distribution and morphology, the boundary conditions, the periodicity, the discretization and the global solver. The LMGC90 library takes into account the FCZM between each body. The material properties of the matrix and the inclusions are managed by the MatLib library.

Body-to-body interaction : influence of interface type on the fracture of a brittle heterogeneous material
This example deals with the fracture of a bimaterial. Each mesh is considered as an independent body (body-to-body interaction, Figure 1 case). In particular, the influence of the interface behavior on the fracture of metal matrix composites is investigated. The considered composite is representative of hydrided Zircaloy-based alloys at high burnup which compose cladding of nuclear fuel rods after many years in Pressurized Water Reactor.
The structure is composed by a metal matrix (Zircaloy-4) and rectangular aligned inclusions (δ-hydrides). The structure geometry is a square with length L = 50μm. Horizontal velocity boundary conditions on the left and right vertical faces are prescribed (respectively V = −1m/s and V = 1m/s). A precrack is introduced at the bottom perpendicular to the loading (Figure 8).
Two cases of bonding strength value between the two phases, strong and weak, are considered. Figure 9 shows the fracture features for the weak and strong matrix/inclusion cases. The crack path is significantly influenced by the interface matrix/inclusions bonding strength. For the weak interfaces, the cracks propagate inside the matrix and along the inclusions boundaries due to the formation of microcracks along weak interfaces. In case of strong interfaces, the cracks propagate through the inclusions due to the high cohesive strength with the matrix. These results are consistent with the following criterion (Raous et al., 2002;Raous et al., 2001;He et al., 1989;Siegmund et al., 1997;Martin et al., 1998;Xu et al., 1998) : the transition between deflection and penetration occurs for a ratio of interface fracture energy w int to inclusion fracture energy w i in the range [0.013, 0.25]. In other words, when w int < 0.013w i the matrix crack is delfected at the interface (case of weak interface) ; when w int > 0.25w i the matrix crack propagates through the inclusion (case of strong interface). Figure 10 shows the evolution of the stress during cracking process. During the loading, the stress increases linearly. When the precrack propagates into the matrix to the first inclusion, the stress becomes to decrease. In case of strong interfaces, it continues to decrease rapidly due to the propagation of the cracks through the inclusions while in case of weak interfaces, it increases due to the formation of microcracks and their coalescence (precrack is arrested at the matrix/inclusion interface).

Figure 8. Boundary condition (arrows), matrix (light gray), inclusions (dark gray) and precrack (white line)
The energy release rate decreases with increasing interfacial bonding strength. The energy release for strong interfaces is lower than for weak interfaces. These results show the importance of the interface bonding strength and, in conclusion, a strong interface is more deleterious than a weak interface.

Body-to-body with periodic boundary interaction : fracture of a ductile heterogeneous material
Since the previous example can be considered as a calculus of a structure behavior, we consider now the same type of heterogeneities but from the point of view of the This example deals with the fracture of a bimaterial with periodic conditions and without precrack (body-to-body interaction, Figure 1 case, with periodic conditions). The example is the same as the previous one (section 4.2.1), except that the boundary conditions are periodic and there is no precrack in the structure. In addition, the bonding stength value between the two phases is assumed to be strong.
A macroscopic strain gradient is prescribed along the horizontal direction.
Without any initial precrack, a multifracture initiates in the inclusions, especially at the locii of high concentration of inclusions ( Figure 11). Then, the cracks propagation occur through growth and coalescence and lead to the failure of the structure. Figure 12 shows the overall stress-strain curve. The overall behavior is ductile over an uniaxial stress about 530M P a (the yield stress of the matrix being 450M P a) : as expected, the presence of elastic inclusions with higher strength than the matrix increases the overall yield stress. This example deals with intergranular fracture in order to show the ability of the method in a cluster-to-cluster strategy. Each grain is a meshed domain related to each other with FCZM (cluster-to-cluster interaction, Figure 2 case). The considered structure is composed by metal grains (Zircaloy) and by two multi-grains inclusions (δhydride). The structure geometry is a periodic square with length L = 100μm (periodic boundary interaction, Figure 3 case), see Figure 13. A periodic Voronoï tessellation is used. A macroscopic strain gradient rate is prescribed along the horizontal direction.

Figure 11. Initiation of cracks (white line) in the inclusions
The bonding strength value between the hard grains and the soft grains is assumed to be weak. Figure 13 (right) shows the fractured periodic structure. Again as expected, the cracks initiate at the grain boundaries of the hard grains. Two main cracks propagate along the grain boundaries of the soft grains and they join together to form a single one which leads to the failure of the material. Figure 14 shows the overall stress-strain curve. The overall behavior is brittle due to the weak bonding strength of the grain boundaries.

Conclusion and perspective
This paper presented a numerical framework for the modeling of dynamic crack propagation in heterogeneous materials. The underlying model rests on a multibody approach and complex interactions behavior in a periodic extension of the standard NonSmooth Contact Dynamics framework. This framework has been rewritten as a dual Schur approach with complex interaction law between subdomains instead of "basic" continuity condition. It can be seen as an extension of the classical nonoverlapping domain decomposition methods. In particular, it allows to model complex in- teractions, discrete media, periodic media, parallelism or multiphysics. The associated developed software has the following capabilities : dynamic, finite strain, non-linear behaviors, periodic Finite Element, crack initiation and propagation, nonsmooth postfracture behavior. Two examples dealing with the influence of the matrix/inclusion interfaces on the crack propagation and the intergranular fracture in a heterogeneous material have been presented to illustrate the ability of Xper and of the modeling strategy.
This strategy and Xper code can be applied to nuclear safety, for example, to determinate the fuel behavior at high burnup in a nuclear reactor or the ageing effect on the behavior of nuclear power plant equipment.
A natural extension of the proposed modeling framework is the parallelism. Each body can be considered as an independent meshed subdomain. In that case, the subproblems are solved on parallel computers. An another application is the multiphysics in which each mathematical problem is posed on a different domain. Furthermore, one can introduce complex interactions law between the subdomains.

A. Overview of the software implementation
The architecture design of the software is based on object oriented techniques. Object oriented programming (OOP) provides a clear modular structure for programs and makes it easy to maintain and modify existing code. The fundamental concepts are (Martin, 1996;Martin, 2003;Meyer, 1997) : inheritance, encapsulation, abstraction and polymorphism. One of the principal advantages of OOP techniques over procedural programming techniques is that they enable programmers to create modules that do not need to be changed when a new type of object is added. A programmer can simply create a new object that inherits many of its features from existing objects. This makes object-oriented programs easier to modify.
For low programming cost, the strategy retained here is not to develop the software 'from scratch' but to reuse and extend existing libraries. It permits to take advantages of the libraries updates while developing the whole software. The framework can be splitted in two distinct levels (see section 3) : a global Finite Element resolution (including periodic extension) and a local non smooth contact resolution. The Xper architecture respects this local/global levels by the coupling of three libraries : LMGC90, PELICANS and MatLib. The libraries are described in the following.

A.1. Existing libraries
The LMGC90 library (LMGC90, 2009) is developed by Dubois and Jean . The software is an open platform under the CECILL License (CECILL, 2005) for modeling interaction problems between elements including multi-physics. It allows to model : granular material made of rigid or deformable bodies and with complex interactions (contact, friction, cohesion, wear, etc.), discrete media, masonry, fracture, etc. The modeling approach is based on an hybrid or extended Finite Element Method (FEM) -Discrete Element Method (DEM), using various numerical strategies such as Molecular Dynamics (MD) or NSCD. In particular, the NSCD algorithm allows to take into account cohesive zone models with contact and friction (see section 3 and 4.1).
The programming language used is Fortran90. Although Fortran90 is a procedural programming language, LMGC90 is implemented in the form of modules using Object Oriented Programming. In particular, the code is open for extension and closed for modification, it respects the 'Open-Closed' principle (Meyer, 1997). LMGC90 can be used as a library or as a standalone software through a macro language.

A.1.2. PELICANS
The PELICANS library (PELICANS, 2009), developed by the French 'Institut de Radioprotection et de Sûreté Nucléaire' (IRSN) under CECILL-C License (CECILL, 2005), is a toolbox for the implementation of various numerical methods dedicated to the solution of systems of partial differential equations (PDEs).
The PELICANS platform is written in standard C++ and is object-oriented. The constitutive classes are classified in two groups : 'plug-points' and 'service provider'. The 'plug-points' classes play the role of mother classes and permit user classes to plug into the platform. The 'service provider' classes make functionalities available to users. The platform respects the Component-Based Development principles, such as the Design by Contract, command-query separation, Liskov substitution principle, inheritance, naming and self-documentation issues (Martin, 1996;Martin, 2003;Meyer, 1997).

A.1.3. MatLib
The MatLib library is developed by Stainier (Stainier et al., 2003). MatLib is a material constitutive models library. It is based on a variational formalism of thermomechanical constitutive updates (Ortiz et al., 1999;Yang et al., 2006). This formalism rests on potential combination naturally involving an object oriented structure. The architecture permits the creation of new models from existing models and their implementation in the library. The models are accessible to users through a common interface.
The library is written in C++.

A.2. Software architecture
The Xper rests on the coupling of the three libraries (see Figure 15) using OO and Mixed Programming.
In the current version of the code, the LMGC90 library plays the role of the master program in the coupling. It manages, in particular, the time discretization and the Newton loops. In the future, the master program will be an interpreted program written with high-level and OO programming language like Python and independent of the three libraries.
The choice of the master program being made, each library has a clear mechanical meaning (Figure 15) : -PELICANS manages the periodic finite element at the global level, -MatLib provides the non linear constitutive models (global level), -LMGC90 manages the surface behavior between the bodies at the local level.
The volumetric behavior is obtained by the coupling of the finite element library developed from PELICANS and the constitutive models library MatLib. Since the two libraries are written in C++ language, the coupling is strong. For example, a MatLib object can be instantiated from a PELICANS class. A PELICANS-MatLib interface permits to manage the MatLib objects in the library developed from PELICANS, preserving the independence of each library.
The surface behavior is managed by LMGC90. The coupling between LMGC90 and PELICANS needs Fortran90/C++ Mixed Programming techniques. Two interfaces are developed for passing parameters from C++ to Fortran90 and back, one in PELICANS and the other in LMGC90.
This strategy permits to compile all the libraries independently, and then to link them.
Note that minor changes have been made in the LMGC90 and MatLib libraries. Note that for the preprocessing, the postprocessing and the linear algebra, the internal capabilities of the PELICANS and LMGC90 libraries are completed by the coupling to external software. These functionalities are thus available in the software.
To reduce the high cost of software development, existing codes have been reused instead of developing them 'from scratch'. According to this strategy, the implementation and the validation of the Xper project, including all the features presented, represents 1 man-year. This is a relatively low development costs compared to developing entirely the features of the libraries in a new platform.