D Dureisseix

J Colmars

A Baldit

F Morestin

H Maigre

Msc

Follow-up of a panel restoration procedure through image correlation and finite element modeling

Keywords: residual stresses, wood panel, painting, conservation, cultural heritage, stereocorrelation, identification, inverse problem

Residual stress estimation is an important question for structural integrity. Since residual stresses are self-balanced stress fields, a classical way to obtain information on them is to remove a part of the structure, and observe the structure displacement field arising from the stress redistribution. The hole-drilling method is such an approach. In some cases, as for the present one concerning a painted panel of cultural heritage, the hole-drilling method is not suited (a structure with a complex geometry, few tests allowed) but one can take advantage of structural modifications if they are monitored (here, a restoration act). We therefore describe in this article a model updating approach, focusing on the residual stress estimation and not on the material parameter identification.

This study couples an optical non-invasive shape measurement (digital image correlation, using a projected speckle pattern on the painted panel, with luminance compensation) and a numerical approach (3D finite elements) for the model updating. The 3D stereo-correlation is used to measure a partial displacement field between three different states of the structure (at three different times of the restoration act). The numerical part concerns stress evaluation, once the model and the experiments are compared using a geometric mapping and a spatial projection of discrete fields. Using modeling and identification, the simulation is used to obtain the residual stresses in the panel, before and after the restoration.

Conservation of paintings of cultural heritage on wood panels could gain from using tools of wood science and structural mechanics, to guide decision for curators and restorers. Indeed, the consequences of a restoration act on the future integrity of an artwork can be assessed with virtual (numerical) simulation once a predictive model has been designed. Such numerical simulations in the same context have recently been used for painted panels, for instance in [START_REF] Dureisseix | Mechanical modeling of the activity of the flexible frame[END_REF][START_REF] Chassagne | Mechanical consequences of hygroscopic variations on wooden panel paintings: development of a predictive model[END_REF][START_REF] Marcon | Hygromécanique des panneaux en bois et conservation du patrimoine culturel[END_REF], for musical instruments, [START_REF] Saft | Supporting the restoration of historical pianofortes by numerical simulation[END_REF][START_REF] Leconte | Attempt to find the ancient sound: from the modelisation to the reconstruction of the Erard Piano in Musée de la musique[END_REF], for other wooden structures, [START_REF] Braovac | The viking ship finds: Preservation challenges and the search for potential solutions[END_REF][START_REF] Chassagne | The contribution of numerical simulation for the diagnosis of the conservation of art objects: Application to Antonio Santucci's armillary sphere of the 16th century[END_REF], for ancient buildings, [START_REF] Mele | Modelling and analysis of a basilica under earthquake loading[END_REF][START_REF] Rafiee | Modelling and analysis of the Nîmes arena and the Arles aqueduct subjected to a seismic loading, using the Non-Smooth Contact Dynamics method[END_REF][START_REF] Valluzzi | Structural investigations and analyses for the conservation of the 'Arsenale' of Venice[END_REF], and even for natural parks [START_REF] Cavagnero | Numerical model application for the restoration of the Racconigi Royal Park (CN, Italy)[END_REF].

Since each artwork is a particular case, it requires an identification step to nurture the model. Objects of cultural heritage are often unique and precious artworks, and few mechanical tests can be conducted to identify the present state of the structure, that cannot be estimated with the evolution it was subjected to, due to the lack of past measurements. The present study mainly deals with a model design, based on finite elements, to couple simulations and experiments during a particular restoration act. The concerned artwork is 'Baptême du Christ', from an anonymous artist, stored in 'Palais du Roure', Avignon, France, Figure 1-left.

An early restoration act, performed in the 70s but typical of the 19th and beginning of the 20th century, was a french parquetage (or cradle) on the rear side of the painted panel. This cradle consists of eight vertical beech beams, glued on the panel rear side (485 mm × 405 mm × 12 mm), and crossed with eight horizontal beech beams, Figure 1-right. This kind of restoration aims to rigidify the wooden support to avoid excessive movements (mainly bending due to dissymmetry in moisture exchanges on both sides of the panel) that may endangers the pictural layer. The drawback of such a rigidification is the increase in internal stresses with humidity variations of the environment, see [START_REF] Rothe | Critical history of panel painting restoration in italy[END_REF]. Moreover this panel exhibits two cracks; they will not be taken into account herein.

The new restoration act dates back to 2007 and has been followed in this Figure 2: Former horizontal beams (left), cut to be extracted (central part missing, used for material identification) and replacement of horizontal beams study: the replacement of the horizontal beams of the cradle by new ones made in spruce. The initial horizontal beams were partly glued and partly clamped due to a permanent long-term deformation of the panel that installed itself after the initial restoration (mainly due to the mechanosorption effect). For the former beams to be replaced, they had to be cut, Figure 2. This illustrates the presence of internal (or residual) stresses. The restoration has been monitored with image analysis: the shape of a part of the front painted side has been measured (i) before the removing of the former horizontal beams, (ii) after this removing and (iii) after mounting the new horizontal beams. These experimental data has to be used in conjunction with a structural analysis, in order to estimate the residual stresses in the panel, which is one of the goals of this study. Section 2 describes the experimental technique, while Section 3 focuses on the finite element model and its comparison to the measurements. Finally, Section 4 uses these tools for the residual stresses estimation.

Shape measurement with 3D stereo-correlation

This non-invasive optical technique allows measuring a 3D shape of a part of the surface of a structure. For artworks of cultural heritage, this technique is useful since no contact with the artwork is needed; nevertheless, no continuous measures can be obtained along time, only several ones at particular instants, due to the delay needed to install and calibrate the acquisition chain which is composed (Figure 3) at least of:

• two stereoscopic cameras, • a video beam projector,

• a processing unit (computer) to capture and correlate images from the cameras.

With two camera images taken from different points of view, the 3D position of visible and identifiable points (or patches) can be obtained by image correlation. The basic images of the artwork cannot be used due to too much contrast on the pictural layer, between large surfaces of too small contrast (aplats or flat tints). Therefore, a more suited pattern (classically a speckle pattern) should be substituted to the original image of the panel. Since no physical speckle pattern can be marked on the painting, a virtual speckle pattern image is projected onto the painted surface, once the initial painting is virtually rubbed out. To do so, an initial image is taken, is numerically treated to produce a 'negative' image that is back-projected to the panel to compensate the initial picture luminance (this is the so-called 'extinction' of the painting). The virtual speckle pattern is then added to the projection to appear on the panel, Figure 4. The correlation of the images taken from this virtual speckle pattern by the two cameras allows to derive the 3D position of patches of pixels, as in [START_REF] Maigre | Image correlation applied to wood painting[END_REF]. This leads to approx. 80 000 3D point locations on a large part of the pictural layer (not up to the border, nevertheless).

Three different measures have been taken:

• one before the restoration act; this artwork state will then be denoted with Ω (1) ;

• one with the horizontal beams removed, which panel state denoted with Ω (2) ; and

• one with the new horizontal beams mounted, denoted with Ω (3) .

Apart from an image distortion correction (with the help of a calibration before any data capture), the position of the centers of pixel patches (which are the measured points) are finally converted into physical length in the coordinate system related to one of the cameras. Since, for each of the measured states, the panel may be not repositioned at the same location, each measure is assumed to posses its own coordinate system.

Ideal geometry and wood behavior

For practical reasons, a finite element model has been design on an ideal geometry (i.e. a perfectly flat and rectangular panel), Figure 5. This particular state of the artwork, though not physically interesting, will be the reference state for computations, and will be denoted with Ω (0) . More precisely, the notation Ω (j) will denote the state (j) on the structure composed of the panel and the vertical beams, while the set of horizontal beams will be denoted with Ω ′(j) when needed. These two structures can be seen as two subdomains, and the interface between them will be denoted with Γ, Figure 5.

Assuming that all the previously mentioned states Ω (j) are not too far from the reference state Ω (0) , the small displacement and small strain assumptions hold, and all the finite element displacement fields will be defined on this reference state, as for the material characteristic coefficients. In this study, we expect the restoration act to be sufficiently short to neglect the relative humidity changes of the environment which the wood is sensitive to. Since only elasticity is modeled herein, other mechanical behavior such as viscoelasticity is neglected for sake of simplicity, though its characteristic time may be of the same order of magnitude that the restoration duration. Nevertheless, this assumption is conservative, since the residual stresses tend to relax in presence of viscoelasticity.

For the panel, once the sawn on the initial trunk is known, the elastic characteristics are selected as for a standard coniferous wood: spruce from [START_REF] Guitard | Mécanique du matériau bois et composites[END_REF], Table 1. The elastic behavior is orthotropic, and heterogeneous (the local anisotropic basis changes with the considered point). The position of the tree center is determined with observation of the growth rings on the RT (radial-tangential) section: an off-plane of 125 mm is obtained, Figure 6.

The beams are assumed to be homogeneous, sawn along the L (longitudinal) direction. Displacement fields with respect to the reference state are denoted . Since state Ω (2) will play a particular role in the following, the displacements with respect to Ω (2) are also of interest: they are denoted with U (j) ; indeed:

U (j) = Ū (j) - Ū (2)
.

Matching experimental points and finite element nodes

Each measured set of points is related to its own coordinate system, which are in turn different from the coordinate system used for the ideal geometry Ω (0) . Therefore each measured set (j) has to be matched to the state Ω (0) in order to define the partial finite element field of measured values Ū (j) m . Note that these measurements are partial information since the displacement is measured only To get these partial finite element fields of measured values, three reference points are selected in order to be easily located on pixelized images, both on the states Ω (j) with coordinates X (j) and on a frontal image identified as Ω (0) with coordinates X (0) , Figure 7.

The transformation mapping X (j) to X (0) is expected to be composed of a translation T (j) and a finite rotation R (j) . The translation links the centroids of the set of reference points. The finite rotation may be defined for instance with a polar decomposition of the mapping; in this case, we consider that the three reference points define a plane into which each point M (j) has barycentric coordinates λ: M (j) = X (j) λ. With independent points, λ = A -1 X (j) T M (j) with A = X (j) T X (j) . The transformation is expected to map M (j) to M (0) = X (0) λ = F (j) M (j) with the gradient of the transformation:

F (j) = X (0) A -1 X (j) T .
A polar decomposition of F (j) produces the rotation R j . Unfortunately, this rotation may be composed with a planar symmetry, which is cumbersome to eliminate, see [START_REF] Nesme | Efficient, physically plausible finite elements[END_REF][START_REF] Schmedding | Inversion handling for stable deformable modelling[END_REF]. We therefore preferred a simple Gram-Schmidt orthogonalization E (j) of the basis generated by X (j) (and E (0) accordingly). In such a case, the previous approach leads to:

F (j) = E (0) E (j) -1 which is exactly the rotation R (j) .
Then, the coordinates of all measured points can be transformed to lie in the coordinate system of Ω (0) . Finally, the normal coordinate to the plane of the panel is interpolated at each possible finite element node of the painted side to get a field of normal displacement z T Ū (j) m , Figure 8, where z T denotes the normal to the panel.

Model updating

Once the partial shape of the different states Ω (j) of the structure composed by the panel and the vertical beams, z T Ū (j)

m or z T U (j) m = z T Ū (j) m -z T Ū (2)
m , have been obtained, these different states have to be reconstructed from the ideal geometry: this is the model updating stage. Since the small displacement assumption holds, the shape updating reduces to find a finite element displacement on the undistorted ideal geometry: Ū (j) (M ), M ∈ Ω (0) . The column vector of the corresponding nodal displacements is denoted with ū(j) . These degrees of freedom (dof) can be split in several sets. First, the dofs at the interface between the panel (with vertical beams) and the horizontal beams are denoted with a subscript Γ. The remaining dofs are denoted with a subscript i; they are themselves split into measured dofs with a subscript c (these dofs are only the out-of-plane component of the displacement at measured nodes), and the non-measured dofs with a subscript r. The measured partial shapes at finite element nodes are stored in column vectors ū(j) m or u

m . Finally, we can define C (j) as the Boolean mapping matrix on c dofs (measured dofs), such that ū(j) c = C (j) ūj . The superscript (j) is omitted in the following.

For each state Ω, a part of its boundary is the interface with the horizontal beams Γ; the remaining part, ∂ 2 Ω, is supposed to be traction-free. We will neglect the body forces due to gravity in the following; indeed, all the measurements are performed on the panel in up-right position, therefore the vertical compressive stress on a cross-section S is bounded with σ = mg/S; with a mass of the artwork m ≈ 3 kg, one gets σ ≈ 6 kPa which will be negligible with respect to the residual stress to be estimated. On the interface Γ, the structure is submitted to the action of the set of horizontal beams; these may consist of a displacement U Γ and a force density F Γ defined on this interface. The techniques that will be used in the following are derived from model identification techniques as in [START_REF] Calloch | Identification de modèles de comportement de matériaux solides : utilisation d'essais et de calculs[END_REF][START_REF] Ladevèze | Updating of finite element models using vibration tests[END_REF][START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF]. A couple of displacements (U , U Γ ) will be said kinematically admissible (KA) simply if U = U Γ on Γ. A couple of stress and interface forces (σ, F Γ ) will be said statically admissible (SA) if: div σ = 0 in Ω, σn = 0 on ∂ 2 Ω, and σn = F Γ on Γ. n is the unitary normal vector to the boundary, pointing out of Ω.

A mechanical state of the studied structure consists of two couples, (U , U Γ ) KA and (σ, F Γ ) SA, satisfying the constitutive relation, here the elastic behavior: σ = Dε(U ) in Ω. In other words, the following constitutive relation error should be null:

Ω e 2 (σ, U )dΩ = 0 (1)
with

e 2 (σ, U ) = 1 2 [σ -Dε(U )] : D -1 [σ -Dε(U )] (2)
When dealing with model updating, the mechanical state has moreover to cope with the measured quantities, here: the measured displacement field z T U m on ∂ m Ω. Both the model and the measured cannot be perfect, therefore one seeks for a compromise between all the constraints to be satisfied. The model updating proposed approach consists in balancing the verification of the constitutive relation and the measures, by searching the couples (U , U Γ ) and (σ,

F Γ ) as: argmin (U , U Γ ) KA (σ, F Γ ) SA Ω e 2 (σ, U )dΩ + α ∂mΩ e 2 m (U )dS (3)
in which the term e m is an error with respect to the measures, for instance:

e 2 m (U ) = 1 2 [z T (U -U m )]ω[z T (U -U m )] (4)
ω is a scalar field of weighting terms, taking its values in [0, 1], to take into account confidence levels in the experimental values. α is a singe scalar coefficient (homogeneous to a stiffness) to balance the two terms in the functional, or the cost function, to be minimized [START_REF] Braovac | The viking ship finds: Preservation challenges and the search for potential solutions[END_REF]. These parameters will be precised in the following.

Note that this problem is not in a closed form until additional informations are stated on the interface fields (U Γ , F Γ ).

Dealing with finite element numerical approximations is easy for the kinematically admissible couples: the displacements (U , U Γ ) are replaced by their finite element approximations counterparts (u, u Γ ). For the stress field, an additional approximation is to assume a special form, deriving from a displacement field V : σ = Dε(V ), which is in turn classically discretized by finite elements into v. D(M ) is the Hooke operator field on the panel and the vertical beams, modeling the elastic behavior of the wood material.

The static admissibility is therefore replaced with its finite element approximation: denoting with f Γ the generalized nodal forces associated to F Γ , (v, f Γ ) will be said statically admissible (SA), if:

Kv = K ii K iΓ K Γi K ΓΓ v i v Γ = 0 f Γ or Kv =   K rr K rc K rΓ K cr K cc K cΓ K Γr K Γc K ΓΓ     v r v c v Γ   =   0 0 f Γ   ( 5
)
where K is the finite element stiffness matrix, arising from Hooke operator D.

Note that this implies:

v = v i v Γ = A Γ v Γ with A Γ = -K -1 ii K iΓ 1 (6)
A T Γ KA Γ = K ⋆ ΓΓ (where a superscript T denotes the transposition) is the so-called Schur complement of K on Γ dofs, and finally the discretized static admissibility is merely a relationship between displacements and forces at the interface:

K ⋆ ΓΓ v Γ = f Γ .
The model updating now consists in finding:

argmin u (v, f Γ ) SA f (u, v) (7)
with

f (u, v) = 1 2 (u -v) T K(u -v) + 1 2 α(Cu -u m ) T ω(Cu -u m ) (8)
u is a kinematically admissible (KA) field, expected to be close to the measurements u m . This is a balance between a smooth field (according to equilibrium equations, thanks to the first term), and a field matching measurements (according to the second term). ω is a diagonal weighting matrix (with entries ω in [0,1]). We may interpret either the field u or v as a smoothing and a prolongation of u m on the whole structure.

The coefficient α can be estimated by using a collocated prolongation of the measurement: û = A c u m , where A c is defined similarly as A Γ , but with the set of measured dofs on Ω replacing Γ:

α = ûT K û u T m Ku m = u T m A T c KA c u m u T m Ku m = u T m K ⋆ cc u m u T m Ku m (9)
Eventually, it can be adjusted by iterating the identification procedure a few times.

The choice of the diagonal weighting ω is part of the modeling. It may take into account the reliability in the measurements (a value 1 is a maximal confidence, a value 0 is a minimal one). With no available quality estimator field of the measurements, we only recall that the correlations are less accurate on the boundary of the measured area. We therefore choose to decrease the value of ω, from 1 in almost all the measured area, except on two layers of elements on the boundary of the measured area, driving it linearly to 0 at the boundary within these two layers. As both α and ω depend on the set of measured dofs, they are different for each measured state.

An equivalent formulation of the updating problem ( 7) is:

argmin u K ⋆ ΓΓ v Γ = f Γ g(u, v Γ ) (10) with g(u, v Γ ) = 1 2 (A Γ v Γ -u) T K(A Γ v Γ -u) + 1 2 α(Cu -u m ) T ω(Cu -u m ) (11)
This problem can be interpreted as: control a structure (panel and vertical beams) with dofs on Γ to make the displacement on the measured set close to the measured values u m . As this will be outlined in the following, this problem may be ill-posed, and may require regularization.

Each time an updating is produced, the distance to the measures can be evaluated with:

η m = (Cu -u m ) T ω(Cu -u m ) u T m u m (12)
3.4 Relaxed state Ω (2) This state, after horizontal beam removal, is similar to the relaxed state in the hole-drilling method used to determine residual stresses. This method has been used for a long time, and more recently, it has been used in conjunction with field measurements in [START_REF] Schajer | Full-field calculation of hole drilling residual stresses from electronic speckle pattern interferometry data[END_REF][START_REF] Nelson | Residual stress determination using hole drilling and 3D image correlation[END_REF][START_REF] Baldi | Full field methods and residual stress analysis in orthotropic material. I Linear approach[END_REF] and with finite elements for analyzing the measures, with or without inverse identification, see [START_REF] Zhang | A boundary element solution of an inverse elasticity problem and applications to determining residual stress and contact stress[END_REF][START_REF] Shaw | A finite-element technique to analyze the data measured by the hole-drilling method[END_REF][START_REF] Schajer | Advances in hole-drilling residual stress measurements[END_REF][START_REF] Di Scalea | Whole-field strain measurement in a pin-loaded plate by electronic speckle pattern interferometry and the finite element method[END_REF][START_REF] Cárdenas-García | Solution of the moiré hole drilling method using a finite-element-method-based approach[END_REF]. Since residual stresses are self-balanced on the considered structure, their energy is null on any displacement change on the same structure. The hole-drilling method therefore relies on a geometric modification. Here, the removal of the horizontal beams may be considered as such a method. Nevertheless, no residual stresses on the configuration Ω (2) can be reached without a deeper material removal, which is not allowed. So, we can only access to additional residual stresses, or equivalently, we may consider residual stresses on state Ω (2) as null:

σ (2) = 0.
The measured maximal average curvature of the painted side is about ρ

m = 0.187 m -1 . (2)
The previous model updating approach is used with f Γ = 0, and so: v Γ = 0 and v = 0. In such a case, the displacement ū(2) mapping Ω (0) to Ω (2) is selected as u of the previous approach, for which Equation (10) leads to:

K(2) u = C (2) α (2) ω (2) ū(2) m with K(2) = K + C (2) (α (2) ω (2) )C (2) T (13)
Since ū(2) = u is the displacement with respect to Ω (0) , note that the measured displacement with respect to Ω (0) is used in the right-hand-side. As soon as

α (2) ω (2) is non null, K (2)
is obviously symmetric, positive, definite (SPD).

This identification is performed with an error estimation η

(2) m = 16.7 %. Figure 9 plots the weighting field ω (2) and the contribution to the error with respect to the measures [e

(2) m ] 2 .
3.5 Initial state Ω (1) This state, prior to the restoration, corresponds to the (assumed perfect) gluing between the panel and the beams. Model updating from Ω (0) to Ω (1) requires the field ū(1) . Alternatively, once ū(2) has been settled, the more meaningful

m = ū(1) m - ū(2)
m on the measured set of state Ω (1) . The maximal average curvature of the painted side is about ρ

(1) m = 0.107 m -1 , therefore the effect of the horizontal beams is to reduce the free curvature of state Ω (2) .

Since perfect gluing is assumed, the force f Γ is a priori non null, and is not known. In the previous updating approach, we therefore consider v Γ as an unknown. The minimization of the functional g in Equation ( 10) with respect to both u and v Γ , after algebraic manipulations using the property

A T Γ K = [0 K ⋆ ΓΓ ]
, leads to: v Γ = u Γ , and: K(1) u = C (1) α (1) ω (1) u (1) m [START_REF] Leconte | Attempt to find the ancient sound: from the modelisation to the reconstruction of the Erard Piano in Musée de la musique[END_REF] with

K(1) = K +   0 0 0 0 α (1) ω (1) 0 0 0 -K ⋆ ΓΓ   (r,c,Γ)
As soon as α (1) ω (1) is non null, the left-hand-side can be shown to be symmetric, positive but semi-definite only, due to the minus sign before the Schur complement on Γ. Its kernel is the set of control dofs u Γ , whose recovering u = A Γ u Γ is null on the measured set, i.e. C (1) A Γ u Γ = 0; therefore, if C (1) A Γ is not injective, the problem is ill-posed. An interpretation is the following: if there are too many control dofs on Γ, or too small measured dofs, the control is not unique. The model should therefore be sufficiently refined on the measured region. An alternative solution, used here, is to regularize the functional to be minimized: a regularizing term, avoiding uncontrolled solutions on Γ, but leaving rigid body motions free (these are filtered with measures) can be added to the functional. The simpler choice is:

h(u, v Γ ) = g(u, v Γ ) + 1 2 u T Γ K ⋆ ΓΓ u Γ (15)
which leads to the problem, similar to Equation (13):

K(1) u = C (1) α (1) ω (1) u (1) m with K(1) = K + C (1) (α (1) ω (1) )C (1) T (16)

Figure 10: Weighting function ω (1) (left) and contribution to the error with respect to the measures [e

m ] 2 (right), for the initial state and u (1) = A Γ u Γ . This identification is performed with an error estimation η

(1) m = 16.6 %. Figure 10 plots the weighting field ω (1) and the contribution to the error with respect to the measures [e 3.6 Final state Ω (3) The measured average maximal curvature is ρ

(3) m = 0.075 m -1 . Therefore, though the new beams are less stiff that the older ones, their effect is to render the panel plane. This can be understood by the fact that the relaxed state of the horizontal former beams was itself permanently bent. Therefore, though the stiffness of the structure is lower with the new beams, this restoration act does not automatically reduce short-term residual stresses.

Up to this point, the final state Ω (3) , after replacement of horizontal beams, does not require additional measurements, since the new beams, made with spruce, are supposed to be known with a perfect initial state without residual stress. The beam replacement can therefore be simulated. In a first step, the matching of displacement in the assembly reads: W (3) = Ū (3) at least along the normal vector n to Γ (for a frictionless contact model) or for all components (for a perfect gluing), where W (3) is the displacement field on the new horizontal beams. The panel and vertical beams strain is measured with reference to Ω (2) , and are therefore ε(U (3) [START_REF] Baldi | Full field methods and residual stress analysis in orthotropic material. I Linear approach[END_REF] ). This leads to a classical elastic problem on Ω (0) with a prestress σ(2) = Dε( Ū (2) ) on the panel and vertical beams, and ε( Ū (3) ) as unknown. This problem is semi-definite positive and the displacement solution is obtained up to an undetermined global rigid body motion (rbm). With this ideal geometry of the new horizontal beams, and a perfect gluing condition, the average simulated curvature is 0.011 m -1 and the error estimation η

) = ε( Ū (3) ) -ε( Ū (
(3) m = 62 % (the unknown rbm of the problem for state Ω (3) is selected in order to minimize η

(3) m ).
To get a more realistic model, one has to take into account the fact that, in order to allow the mounting of the new beams, their height had to be reduced when compared to the old ones. With an estimated height reduction of 1.25 mm, the simulation of the state Ω (3) is performed with an elastic 2-body problem with Figure 12: Numerical results (coarse discretization): residual stresses / MPa on the panel (displacement amplified 10 times). From left to right: Ω (1) ; Ω (2) ; Ω (3) , perfect assembly; Ω (3) , contact with gap restoration is 18.6 MPa, and is 1.7 MPa after restoration (for the assembly with 1.25 mm gap and the finer discretization).

Conclusions and outlooks

In this article, we exemplify that a restoration act on a painted panel of cultural heritage can be considered as a mechanical test, and that with field measurements and modeling, information can be drawn on the structure. Experimental data and numerical modeling are used together to analyze this test, with a model updating approach. This allows to estimate the residual stresses in the structure before and after the restoration.

Though it is the case here, the replacement of old horizontal beams with less stiff new ones does not necessarily decrease the residual stress, since the old beams may exhibit permanent curvature. An other solution would have been a dedicated shape of the replacement parts on an artwork, as in [START_REF] Ciatti | Dipinti su tavola, la tecnica e la conservazione dei supporti[END_REF][START_REF] Uzielli | Conservazione dei dipinti su tavola[END_REF]. Indeed such a shape could be designed for the new beams: for instance, taking only into account a curvature ρ The somehow large discrepancies between the model and the measurements for the restored panel may arise from different sources:

• the model does not take into account 2 cracks on the panel;

• the mounting of the new beams can have damaged the panel;

• a too rough estimation of elastic parameters (the identification performed here concerns the shape updating, or the residual stress identification, and not the material characteristics);

• the neglected viscoelastic behavior.

The final replacement of the horizontal beams of the cradle are used as a verification step for the model. An other interesting tool to check the validity of the model would be a local contribution to the error in order to locate the areas where the model is incompatible with the measures. This could be the role of the term e 2 =1 2 (uv) T K(uv) in the cost function to be minimized, if the matching to the measures on the restored artwork is also performed.

Once established, such a model can serve as a predictory tool for assessing risks on the artwork, useful for restorers and conservators, by allowing virtual simulation of the influence of several conservation conditions, and possibly several restoration acts.

Figure 1 :

 1 Figure 1: Baptême du Christ, Palais du Roure, Avignon, France (left) and its cradle (right)

Figure 3 :

 3 Figure 3: Principles of the measurement technique (left) and stereo cameras (right)

Figure 4 :

 4 Figure 4: Left camera, after beam removal: initial image capture (left) and projected speckle pattern on the luminance-compensated image (right)

Figure 5 :

 5 Figure 5: Finite element model of the ideal geometry for the panel and vertical beams Ω (0) (left two images), interface Γ (center), horizontal beams (right) Ω ′(0)

Figure 6 :

 6 Figure 6: Determination of local anisotropic basis with growth ring observation

Figure 7 :

 7 Figure 7: Reference points (left) and measured field projected onto the finite element mesh (right) for the initial state, with its own reference points

Figure 8 :

 8 Figure 8: Measured fields projected onto the finite element mesh: initial state (left), after beam removal (middle), and after beam replacement (right). Displacement amplified 8 times.

Figure 9 :

 9 Figure 9: Weighting function ω (2) (left) and contribution to the error with respect to the measures [e (2) m ] 2 (right), after beam removal

( 3 )

 3 m for sawing them, could decrease again the residual stresses.

Table 1 :

 1 Material elastic parameters after[START_REF] Guitard | Mécanique du matériau bois et composites[END_REF] (for standard relative humidity H and specific density ρ)

	species	spruce beech
	H / %	13.1	9.4
	ρ / g/cm 3	0.31	0.63
	E R / MPa	816	2 040
	E T / MPa	304	867
	E L / MPa	8 020 14 100
	G RT / MPa	48	500
	G T L / MPa	461	980
	G RL / MPa	558	1 850
	ν RT	0.67	0.73
	ν LT	0.33	0.46
	ν LR	0.34	0.36

http://www.woodculther.org

Acknowledgements

The authors wish to thank Emmanuel Maurin (LRMH, Laboratoire de Recherche des Momunents Historiques) for fruitful discussion on the restoration technique, Daniel Jaunard and Patrick Mandron (restorers), and Marie-Claude Leonelli (DRAC, Direction Régionale des Affaires Culturelles) for getting access to the panel painting. This work has been supported by European COST Action IE0601 'Wood Science for Conservation of Cultural Heritage' 1 .

m = 26 %. Figure 11 plots the weighting field ω (3) and the contributions to the error with respect to the measures [e

m ] 2 , depending on the assembly model.

Estimation of the residual stresses

With the previous developments, and still assuming small perturbations in all of the following, the computed residual stresses on the panel are easily obtained:

The von Mises component of these residual stresses are depicted in Figure 12. Table 2 recalls the obtained numerical values, for the present discretization with n e = 7788 elements, and with a finer discretization of n e = 35094 elements (two thirds of them are 8-node cubes, the remaining ones are 6-node prisms). The residual stress estimation is still satisfactory, except for the perfect gluing state Ω (3) whose model is erroneous.

Though it is not a systematic result, the present restoration acts decreases the residual stress in the panel: the maximum von Mises stress σ V M before


V Visseq

Alexandre Martin

Damien Iceta

Emilien Azema

D Dureisseix
email: david.dureisseix@insa-lyon.fr

P Alart

Dense granular dynamics analysis by a domain decomposition approach

Keywords: domain decomposition method, nonsmooth contact dynamics, discrete elements, frictional contact, granular packing, LMGC90 PACS: 45.70.Cc, 83.10.Rs, 83.80.Fg MSC: 65Y05, 70E55

In this article a domain decomposition approach is combined with the nonsmooth contact dynamics approach for analysing the global behaviour and the micromechanical structure of large-scale dense granular systems. Previously introduced and theoretically investigated, this method is herein investigated precisely on two aspects: the properties of the interface operators, especially when applied to the corners of the subdomains, and the influence of the substructuring on the solution of a mechanical test. Such topics are specific to the dense granular systems characterized by the discreteness and the nonsmoothness of their behavior.

Introduction

In connection with a domain decomposition strategy, the granular dynamics reveals two main features: discreteness and nonsmoothness.

The non-overlapping decomposition of a granular domain is all the more delicate since such a medium is a non-structured discrete system. Contrary to trusses or tensegrity structures studied in [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF] for which an elementary pattern may be defined during the whole process, a granular system involves a permanent evolution of the connectivity of the particles, specially when granular flows occur. Consequently a box-like partitioning insuring a locality of data, useful for a parallel implementation, provides two possible approaches. A primal strategy leads to a 'non-perfect' interface between the subdomains made of nonsmooth interactions. Because such a method is a simple algebraic partition of the equations and is easy to implement, it has been applied to an industrial problem, the simulation of railway ballasts [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF]. However when some large rigid bodies constitute the boundaries of several subdomains of the system (as the sleepers on a railway track), the size of the interface increases drastically. The dual strategy is less intuitive because it requires to split the grains at the interface. Contrary to the primal approach the interface behavior is now 'perfect', in the sense that only linear equations are describing it (local equilibrium and velocity continuity). Indeed we have to glue the interface grains by adding a new, but linear, equation which modifies strongly the global nonlinear (nonsmooth) solver and which complicates the implementation. However this dual approach has two advantages: (i) the occurrence of large rigid bodies do not perturb the size of the interface; (ii) the perfect boundary of each subdomain should allow to introduce an automatic homogenization process to switch possibly from a discrete model toward a continuous model. This second approach is detailed in the following.

Once the sub-structuring has been performed, a nonsmooth solver has to be combined with the domain decomposition strategy. The nonsmooth relations are derived from the NonSmooth Contact Dynamics (NSCD) approach which is well suited to the simulation of granular systems. NSCD or Contact Dynamics in short, has been developed by J. J. Moreau and M. Jean over the last two decades [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Numerical aspects of sweeping process[END_REF]. It is suited to many applications but has proven to be particularly useful when collections of rigid or deformable bodies are packed together in a dense assembly and subjected to dynamic loading deformations. Numerical simulations have to be performed using a fully implicit resolution of the contact impulses. This allows us to deal properly with nonlocal momentum transfers involved in multiple collisions, contrary to classical molecular dynamics schemes that consider the system evolution as a succession of binary collisions. The approach proposed by Cundall [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF], inspired from the molecular dynamics, as the event-driven methods, requires very short time steps for accounting for the successive collisions and such a strategy leads to a prohibitive computational cost, especially for the dense granulates with a large number of simultaneous collisions. The NSCD approach refers to a time-stepping method that requires at each time step the solution of nonsmooth equations by an iterative solver. The computational cost may be quite high, but the gain is substantial. Simulations of very large granular systems can range from 10 m of a ballast railway submitted to cyclic dynamic loading, to the behavior of the Nîmes arena and Arles aqueduct (France) subjected to seismic loading, which are examples of two challenges in computational mechanics.

The Domain Decomposition Methods (DDM) in the context of multiprocessor computations are well established from theoretical and practical standpoints when dealing with a linear system derived from a discretization of a continuous problem [START_REF] Tallec | Domain-decomposition methods in computational mechanics[END_REF]. For nonlinear continuous problems the DDM seems to be efficient when it is used only to solve an intermediate linear problem embedded in an iterative process as a Newton type method [START_REF] De Roeck | A domain decomposed solver for nonlinear elasticity[END_REF]. Unfortunately the simulation of the granular systems with nonsmooth interactions between grains does not use such a nonlinear solver and the combination with a DDM has to be rethought. Indeed the NonLinear Gauss-Seidel (NLGS) algorithm may be considered as the generic nonsmooth solver because it allows to embed a large range of interaction laws such as adhesion, cohesion, capillarity... without modifying deeply the algorithm. In line with the NSCD approach, the velocity-impulse formulation is extended herein to a multidomain reformulation preserving the generic algorithm. More precisely the multidomain reformulation is based on a FETI-type approach where the subdomains are 'glued' by Lagrange multipliers which are inter-domain forces. This choice is made in accordance with the NSCD approach where the impulses are privileged variables. The so-called Non-Smooth Contact Domain Decomposition (NSCDD) solving method consists of a two-stage algorithm. One of these stages recovers the generic NLGS method applied subdomain per subdomain in conjunction with the NSCD formulation; for details about convergence, refer to [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF]. This generic algorithm is presented in Section 2 and a theoretical study of the convergence is developed in [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF].

The DDM introduces different types of interface according to their dimension. For a three-dimensional continuous problem, we distinguish facets, edges and corners. Specific strategies are developed for dealing with the corners resulting from the domain decomposition of structures discretized by a finite element method [START_REF] Farhat | FETI-DP: A dual-primal unified FETI method -Part I: A faster alternative to the two-level FETI method[END_REF][START_REF] Lesoinne | A FETI-DP corner selection algorithm for three-dimensional problems[END_REF]. For discrete systems the distinction is less clear and we have developed in [START_REF] Alart | A scalable multiscale LATIN method adapted to nonsmooth discrete media[END_REF] the concept of 'weak' interfaces in the context of static problems solved by a LATIN type method. We investigate in Section 3 the features of the interface problem solved at the second stage of the generic algorithm when some grains, located at the corners, are connected to more than two subdomains.

Finally the evaluation of the efficiency of a new multidomain solver in comparison with a previous monodomain one is a difficult topic because a dense granular system is an evolutive nonsmooth problem leading to a large multiplicity of solutions [START_REF] Moreau | Indetermination in the numerical simulation of granular systems[END_REF]. Consequently we have not relevant error estimates as underlined in [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. Only the quality of the computation may be appreciated using a set of qualitative indicators as presented in Section 4 for studying the global behavior and the micromechanical structure of a granular sample submitted to a biaxial test. A time consuming analysis is performed on a Sequential Multidomain implementation in order to estimate the CPU time gain that we can expect from a multiprocessing implementation on a distributed memory architecture. The present approach has been implemented into the LMGC90 platform [START_REF] Dubois | LMGC90. In 10e colloque national en calcul des structures -CSMA2011[END_REF].

2 Dual domain decomposition method for granular systems

Reference problem

When a time-stepping scheme is used, we denote known quantities at the beginning of the time slab [t i , t i+1 ] with a superscript (-); the quantities at the end of the time slab (without a superscript) have to be determined.

Figure 1: Coordinate basis; R g : global coordinate basis (for interface quantities), R i : local (related to the grain i) coordinate basis (for grain dynamics), R α : local (related to an interaction α between two grains) coordinate basis (for interactions). Details of the local contact frame (n, t) at a contact α between two touching particles i and j.

Grain nonsmooth dynamics. Considering a rigid model for the grains, the dynamics of the granular medium is written as the following vector equation

[34] M V -R = F d (1)
where the prescribed right-hand side is

F d = R d +M V -.
V is the velocity of the grain (it contains the translational degrees of freedom, and the rotational ones); R is the resultant impulse on the grain due to interactions with other grains. The matrix M contains both the mass (for the translational degrees of freedom) and the inertia (for the rotational degrees of freedom). A choice leading to get a constant, and diagonal, matrix M consists in expressing the global coordinates of rotation vectors in the inertia eigenbasis of each grain, Figure 1. With such a choice, R d = R ext + R rot , where R ext are the prescribed external forces and R rot are the fictitious forces defined as

R rot = 0 ω × Iω ,
where I is the moment of inertia and ω is the angular velocity. These fictitious forces are non linear with respect to the degrees of freedom. In the case of dense granular media, angular velocities are small enough to express R rot in an explicit way, by choosing the value obtained at t i , as reported in [START_REF] Radjaï | Discrete-Element Modeling of Granular Materials[END_REF]. This renders them explicitly known, and so, they are assembled into the right-hand side R d , that contains also the prescribed impulse fields. The assembly of these equations (independent for each grain) for all the involved grains is formally written in the same way [START_REF] Acary | Contribution à la modélisation mécanique et numérique des édifices maçonnés[END_REF].

Contact interaction.

Here we focus on simple unilateral contact which is naturally expressed as a complementary condition linking contact force to a gap. For dynamics, Moreau proved via a viability lemma [START_REF] Moreau | Numerical aspects of sweeping process[END_REF], that we can use a velocity-impulse complementary law. The constitutive relation is summarized in the following formal equation:

R(v, r) = 0 (2)
v is the velocity jump at the contact point between the two interacting grains, r is the impulse at the same contact point. R is usually a non linear and multivalued relationship between the previous two dual quantities. For instance, for frictionless contact, it relates the normal components with the KKT complementary conditions v n ≥ 0, r n ≥ 0 and v n f n = 0. Other models, such as extensible cables and frictional contact can be found in [START_REF] Nineb | Domain decomposition approach for nonsmooth discrete problems, example of a tensegrity structure[END_REF][START_REF] Saussine | Modelling ballast behaviour under dynamic loading, part 1: A 2D polygonal discrete element method approach[END_REF]. The assembly of the interaction-related quantities for all interactions is also written formally in the same way [START_REF] Acary | Numerical simulation of monuments by the contact dynamics method[END_REF]. Both v and r are expressed in the local coordinate basis to the contacts between the interacting grains, Figure 1. Therefore, they are linked to the global kinematic and static quantities with a compatibility condition:

v = H T V and R = Hr (3)
H is the signed mapping between the global quantities related to the grains in their local basis R i with the local relative quantities related to the interactions in the local basis R α .

Reduced dynamics. Taking the dynamics (1) and the compatibility conditions (3) into account, the reduced dynamics involving material variables can be obtained:

W r -v = -v d ( 4
)
where W is the Delassus operator: W = H T M -1 H, and

v d = H T M -1 F d .
To close the problem, one adds the constitutive relation [START_REF] Acary | Numerical simulation of monuments by the contact dynamics method[END_REF], and the reference problem reads:

W r -v = -v d R(v, r) = 0 (5)
This problem is classically solved within the NSCD (NonSmooth Contact Dynamics) method with a non-linear Gauss-Seidel (NLGS) solver [START_REF] Moreau | Numerical aspects of sweeping process[END_REF][START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF].

Extension to deformable grains. Even if it will be not tested in the numerical results presented in the following, the case of modeling the grain behavior as an elastic deformable solid, with a finite element discretization, can be derived easily. This leads to even more large problems for which a domain decomposition method has also potentials for larger gains. Such a modeling is suited in particular for granular materials where deformation and eventually fracture of grains is under concern. Now, the kinematic global unknowns V are the whole set of translational degrees of freedom of the nodes, K is the classical finite element stiffness matrix of the grain and M is the mass matrix of the grain. Some care must nevertheless be taken with this mass matrix to get a discretized well-posed problem, see for instance [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF].

A co-rotational formulation [START_REF] Acary | Contribution à la modélisation mécanique et numérique des édifices maçonnés[END_REF][START_REF] Acary | Numerical simulation of monuments by the contact dynamics method[END_REF] has several advantages: if the rotations are finite, but the deformation is small, expressing the degrees of freedom in the inertia eigenbasis of the grain allows to consider constant operators M and K. In such a case, as previously, the Coriolis and centrifugal effects are explicitly written, and are part of the given right-hand side of generalized forces F d .

A two-scale description consists in setting V = R s V s + V e where V s is the previously described small-sized vector of the global rigid body movement of the grain, at its center of mass. R s is the extension of this movement to all the nodal degrees of freedom of the discretized body (in its inertial eigenbasis). V e is an additional movement mainly containing the elastic deformation of the grain (in its inertial eigenbasis as well), to be precised in the following. The non smooth dynamics of the grain therefore reads:

M V + R int -R = F d (6)
with the internal impulse

R int = ti+1 ti KU dt ( 7
)
where U is the nodal vector of the displacement in the inertial eigenbasis). Using the test functions V = R s V s + V e , the dynamics leads to

M V e + R int + M R s V s -R = F d (8) M V s + R T s R int + R T s M V e -R = F d (9)
where

M = R T s M R s , R = R T s R, F d = R T s F d . With a constant stiffness matrix, one coupling term is R T s R int = R T s K ti+1 ti U dt = 0
since KR s = 0, and, to ensure the uniqueness of the two-scale description, we choose as an orthogonality condition between the two kinematic spaces: R T s M V e = 0 that cancels the second coupling term in [START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF], which is therefore identical to the rigid model ( 1). Once V s is obtained, the "deformable" component V e then arises from [START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF]:

M V e + R int -R = F d -M R s V s
The last step is to link the displacement update U to the velocity V . A possibility is to obtain it from two sources: U = U s + U e , where U s is a rigid body rotational finite movement (useful to update the finite rotation of the inertial eigenbasis), and U e corresponds to the complementary part. U s can be obtained with the Hughes-Winget scheme [START_REF] Hughes | Finite rotation effects in numerical integration of rate constitutive equations arising in large deformation analysis[END_REF] or the Rodrigues formula, while U e can be obtained with a θ-method as a time integration scheme. Neglecting the residual KU - s in algebraic developments, one gets the internal impulses as

R int = hKU - e + hθ[h(1 -θ)KV - e + hθKV e ]
where the time step is h = t i+1 -t i . The corresponding dynamics therefore reads:

M V e -R = F d -M R s V s ( 10
)
where F d is a given right-hand side, and M = M + h 2 θ 2 K. One can check that, with the coupling term M R s V s , the dynamics [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] gives a solution that satisfy to the orthogonality condition. Indeed, by pre-multiplying (8) with R T s , one can easily check that it leads to R T s M V e = 0. Finally, the reduced dynamics can also be drawn for this model case, as previously with a new definition of matrix H which is now correlated to the whole set of degrees of freedom, but with the same expression as in (3):

v = H T (R s V s + V e ) and R = Hr, R = R T s Hr (11) which leads to W r -v = -ṽ d R(v, r) = 0 (12)
with, in fine, W = H T M -1 H and ṽd = H T M -1 F d . Therefore, the problem characteristics are very close to the ones obtained from the case of rigid grains and algorithmic choices in a domain decomposition approach should be valid for both modelings.

Domain partitioning

The domain has to be split into subdomains in order to use parallel computing. This decomposition is performed as frequently as needed to take into account the migration of grains from one subdomain to another. Such a strategy may be implemented with minimal computational efforts using sophisticated routines out of the purpose of this paper. Since the nonsmoothness may occur in interactions between grains, we choose to distribute interactions among subdomains as in [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF] (we proceed by distributing the middle points between the centers of mass of interacting grains, according to their coordinates, using an arbitrary regular underlying grid, Figure 2 and Figure 3). Indeed, with such a choice, the "boundary" grains are duplicated in the two subdomains. If a grain indexed with i is connected with m i subdomains, m i is called its multiplicity number.

For consistency for the rigid model of the grains, the masses and moments of inertia are distributed among the neighboring subdomains according to their multiplicity number. More precisely the distribution of masses and inertia is an algebraic partitioning and not a geometrical partitioning. It is more meaningful to speak about a duplication of the boundary grains than a splitting of them.

For the elastic deformable model of the grains, this splitting can be performed with a classical mesh decomposer. The interface between two subdomains is defined as the set of these grains, that joins the subdomains. The nonsmoothness is therefore localized only within the subdomains. This modeling choice is identical to [START_REF] Barboteu | A domain decomposition strategy for nonclassical frictional multi-contact problems[END_REF] and somehow the dual of that proposed in [START_REF] Ladevèze | A multiscale computational approach for contact problems[END_REF], where nonlinearities (contact on crack lips) are isolated at the interfaces. Note that a direct algebraic partitioning of the reference problem can also be chosen, leading to a dual partitioning and a different algorithm [START_REF] Hoang | A domain decomposition method for granular dynamics using discrete elements and application to railway ballast[END_REF]. Some advantages and disadvantages have been mentioned in the introduction but such a topic has to be investigated more deeply in forthcoming works.

FETI-like formulation and NSCDD algorithm

In each subdomain E, the problem is identical to the global one (with the subscript E), provided that a term arising from the inter-grain interface is added.

It is therefore described in the following only for the rigid model of the grains. It can be built from the interconnecting condition (on the velocities of boundary grains) that has been added to "glue" neighboring subdomains, where A ΓE is a signed boolean matrix with a finite rotation, to map the grain velocities V E to the global coordinate basis into which the interconnectivity is expressed: Γ denotes the global interface between all the neighboring grains. Formally the previous summation is performed on all the subdomains (number equal to n sd ), even if, for a given interface grain the only neighboring subdomains have to be considered. Then we obtain a FETI-like formulation [START_REF] Dureisseix | A numerically scalable domain decomposition method for the solution of frictionless contact problems[END_REF] for the reference problem using a multiplier field F Γ and the notations ÂT

n sd E=1 A ΓE V E = 0 (13)
ΓE = H T E M -1 E A T ΓE : W E r E -v E -ÂT ΓE F Γ = -v d E R(v E , r E ) = 0 E = 1, . . . , n sd n sd E=1 A ΓE V E = 0 (14)
The reduced problem on (r

E ,v E ,F Γ ), with the notations f = E A ΓE M -1 E F d E and X = E A ΓE M -1 E A T ΓE , reads: W E r E -v E -ÂT ΓE F Γ = -v d E R(v E , r E ) = 0 E = 1, . . . , n sd XF Γ - n sd E=1 ÂΓE r E = f ( 15
)
As for many domain decomposition approaches, the goal is to be able to localize the same typical problem that is under consideration on each subdomain independently, while designing a suited coupling recovery algorithm between subdomains, i.e. on the interface.

Here, the algorithmic formulation described in Algorithm 1 has been implemented into the LMGC90 platform [START_REF] Dubois | LMGC90. In 10e colloque national en calcul des structures -CSMA2011[END_REF] for time-evolution problems (N is the number of time steps). At each new time step of the incremental solving procedure, the mapping H and the contact graph have to be reevaluated within a contact detection phase. Eventually, the domain could also be repartitioned according to the new contact graph. For each time step, the iterative resolution proceeds with several stages. First, the interface impulses obtained at the previous iteration are disassembled into F E = -A T ΓE F Γ that is considered as given additional external impulses in the subdomain E, and added to the prescribed values F d E . At each iteration, the solver is itself a predictor-corrector scheme, for which a "free" grain velocity is first computed as

V d E = M -1 E (F d E + F E )
. At the interaction level, one then computes vd

E = H T E V d E .
The correction phase is composed with an incomplete solving procedure of the nonsmooth dynamics problem on each subdomain, with n GS prescribed iterations of non-linear Gauss-Seidel algorithm. This provides the local impulses rE (satisfying the reduced dynamics, even if the solve is incomplete). The resultant impulse per grain is RE = H E rE , and the correction reads: VE = V d E + M -1 E RE . Up to this point, it is interesting to note that the interface problem in [START_REF] De Roeck | A domain decomposed solver for nonlinear elasticity[END_REF] can be stated in a correction form, using ∆F Γ = FΓ -F Γ : noting that rE satisfies to the grain dynamics, this interface problem can be restated as

X∆F Γ = f + n sd E=1 ÂΓE rE -XF Γ = n sd E=1 A ΓE VE ( 16
)
the last term being the residual on the interface, i.e. the velocity jump. At each time step, inner iterations are stopped when the classical NLGS convergence criterion (on the subdomains, [START_REF] Cambou | Micromechanics of Granular Materials[END_REF]) and the gluing criterion over the interface [START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF]

W E rE -vE = -v d E R(v E , rE ) = 0 Update (r E , v E ) ← (r E ,

Interface topics

For discrete systems, several specificities of the interface treatment are detailed in this section, in particular the structure of the interface operator X.

Structure of the interface

For discrete systems, the global interface Γ is constituted of grains supporting contacts in more than one subdomain. The number of subdomains a grain i is connected to is called its multiplicity m i . As for classical domain decomposition methods [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF] one gets m E = 1 + A T ΓE A ΓE as the diagonal matrix whose entries for each grain i kinematic dof is m i . Depending on their multiplicity, the grains are called "internal grains" if m i = 1 (otherwise, they will be called "interface grains"), "face grains" if m i = 2 and "corner grains" if m i > 2. Contrary to the continuous media case where face, edge and corner nodes can be distinguished in 3D, the discrete systems we are considering here do not differentiate edge and corner topology.

Corner grains are split in m i parts and links are stated as gluing conditions between these parts (the impulses in such gluing conditions are stored in F Γ ). Sufficient gluing conditions should be stated for each interface grains (face or corner) to ensure to recover the reference problem solution. Several options are:

• Discard the treatment of corner grains. This option can be used for several DDM for the continuous media case when interface fields are defined at the finite element level on edges of elements rather than at nodes [START_REF] Dureisseix | Une approche multi-échelles pour des calculs de structures sur ordinateurs à architecture parallèle[END_REF][START_REF] Bernardi | A new nonconforming approach to domain decomposition: The mortar element method[END_REF] since the measure of corner nodes is zero. For a discrete model as considered herein, this is not an available option since then, continuity at these grains won't be taken into account.

• Consider as many gluing conditions as the mutiplicity of the considered corner grain: n i l = m i . In this case, there is a small overconstrained gluing condition (only m i -1 links are sufficient to glue m i parts).

• Consider an even larger number of gluing condition, similarly to redundant corner conditions in FETI methods [START_REF] Rixen | A simple and efficient extension of a class of substructure based preconditioners to heterogeneous structural mechanics problems[END_REF]. The maximal number of conditions that can be established between m i parts is 1 2 m i (m i -1). In order to avoid singularity of the interface operator X, and to allow several solving procedures for the interface problem, we choose to prescribe the necessary and sufficient number of gluing conditions on corner grains, i.e. m i -1 conditions only.

i l = m i -1 = 1 2 m i (m i -1) = 1; (b) m i = 3 and n i l = m i = 1 2 m i (m i -1) = 3; (c) m i =

Analysis of the interface operator X

Internal grains (m i = 1) indeed have no contribution to X. When only face grains exist (m i = 2), X has been proved to be diagonal per block, i.e. each grain is decoupled with the other ones [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF] and each block is at most of size (b, b) where b is the number of kinematic or static components for each grain: b = 3(D -1) for rigid grains where D = 2, 3 is the dimension of the considered physical space. In a dynamic framework "sthenic" may be preferred to "static" in reference to a description of dual variables to the kinematics. This is a specific issue of the dynamics of rigid grains: the interface problem does not condense information from the inner part of subdomains. Moreover, each block may itself be diagonal for special cases (circular or spherical grains). In these cases, solving the interface problem is trivial. Corner grains do indeed modify the structure of this operator. It is still block-diagonal, but a full block occurs on links related to each corner grain i with a size (n i l b, n i l b). Consider the contribution of corner grain i from the various subdomains that share this grain; this corresponds to a block in matrix X denoted as:

X i = n sd E=1 A i ΓE (M i E ) -1 (A i ΓE ) T
There are no gluing impulse F Γ between different corner grains, so blocks X i are decoupled to each other. M i E is the mass matrix for the part of grain i located in subdomain E. Since they are always symmetric, positive definite, and for sake of simplicity, they are omitted in the following developments (they are all considered as identity matrices, which does not change the structure of the interface operator), moreover, all the following matrix entries will correspond to a block of size (b, b) corresponding to a whole set of kinematic or static (or sthenic) components. With an arbitrary sign convention in A i ΓE , and for n i l = m i -1 gluing conditions, which is an open cyclic connectivity graph, the block X i is a permutation of the following pattern:

X i =       2 -1 -1 . . . . . . . . . 2 -1 -1 2       (17)
which is clearly invertible. Adding an additional gluing condition will close a sub-cycle in the connectivity, and therefore produces a sub-block in X i with the pattern

        2 -1 1 -1 . . . . . . . . . 2 -1 -1 2 1 1 1 2         (18)
which turns out to be singular (last row is the sum of all the previous ones). Since this is a minor of X i , this last matrix will be singular as soon as n i l > m i -1. The convergence study of the NLGS algorithm is quite complicated and will not be discussed herein. Nevertheless, simulations using corner grains will exemplify the impact of those grains on the physical properties of the granular system are discussed in the following.

A mechanical study as a validity test

The issue of this section is to test the robustness of the DDM approach with respect with various "well-known" aspects of the mechanical behavior of model granular media. To do so, at a first hand, we describe the numerical samples, and compare the macroscopic response of sheared granular packings for different decompositions. The microstructure (i.e. the spatial organization of the particles and their contacts) is analyzed at a second hand as a function of the number of subdomains.

Simulation of a biaxial test

A dense packing composed of 12000 disks is first set up by means of a layerby-layer deposition model based on simple geometrical rules [START_REF] Voivret | Space-filling properties of polydisperse granular media[END_REF]. The particles are deposited sequentially on a substrate. Each new particle is placed at the lowest possible position on the free surface as a function of its diameter. This procedure leads to a random close packing in which each particle is supported by two underlying particles and supports one or two other particles. To avoid long range ordering a small polydispersity in size is used. Following this geometrical process, the packing is compacted by isotropic compression inside a rectangular frame of dimensions l 0 ×h 0 in which the left and bottom walls are fixed, and the right and top walls are subjected to a compressive stress σ 0 . The gravity g and friction coefficients µ between particles and with the walls are set to zero during the compression in order to avoid force gradients and obtain isotropic dense packings, see Figure 5(a). The isotropic samples are then subjected to a vertical compression by downward displacement of the top wall with a constant velocity v y for a constant confining stress σ 0 acting on the lateral wall, see Figure 5(b). The friction coefficient µ between particles is set to 0.35 and to zero with the walls. Since we are interested in quasistatic behavior, the shear rate should be such that the kinetic energy supplied by shearing is negligible compared to the static pressure. This can be formulated in terms of an inertia parameter I defined as [START_REF] Gdr-Midi | On dense granular flows[END_REF]:

I =        ε m p in 2D, ε m pd in 3D, ( 19
)
where ε = ẏ/y is the strain rate, m is the mean particle mass, d the mean diameter and p is the mean pressure (defined as the force per unit width for the 2D case). The quasistatic limit is characterized by the condition I 1; in the proposed simulations, I is below 10 -4 . This simulation is repeated 5 times, with various numbers of subdomains ranging from 0 (corresponding to the reference simulation) up to 4 (they are tagged in the following with S0 up to S4). Figure 6(a) depicts the five decompositions we choose at the initial state. Figure 6(b) is a zoom of the case S4 to illustrate the particle arrangement.

Macroscopic behavior

In this section, we consider the stress-strain and volume-change behavior according to the domain decomposition. We therefore need to evaluate the stress tensor and solid fraction during deformation from the simulation data at microscopic scale.

Definition of some macroscopic parameters

In granular media, the expression of the stress tensor σ in the volume V is an arithmetic average involving the branch vectors α (joining the centers of the two neighbouring particles) and the contact force vectors f α at contact α. It is given with [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF][START_REF] Staron | Friction versus texture at the approach of a granular avalanche[END_REF]:

σ = 1 V α∈V f α ( α ) T (20)
Under biaxial conditions with vertical compression, we have σ 1 ≥ σ 2 , where the σ k are the stress principal values. The mean stress is p = (σ 1 + σ 2 )/2, and the stress deviator is q = (σ 1 -σ 2 )/2. According to the Mohr-Coulomb model, the shear strength of dry granular materials can be linked to the internal angle of friction ϕ as follows [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]:

sin ϕ = q p = σ 1 -σ 2 σ 1 + σ 2 (21)
The vertical macroscopic strain ε 1 is the cumulative value defined as:

ε 1 = h h0 dh h = ln 1 + ∆h h 0 ( 22
)
where h 0 is the initial height and ∆h = h 0 -h is the total downward displacement.

Shear strength and solid fraction

Figure 7 depicts the internal angle of friction sin ϕ as a function of the shear strain ε 1 for all the decompositions. The jump observed at ε 1 = 0 reflects both the rigidity of the particles and the large initial solid fraction of the samples. In all cases, the shear stress passes by a peak (q/p) peak ∼ 0.38 before relaxing to a stress plateau (q/p) * ∼ 0.28 corresponding to the so-called "residual state" in soil mechanics [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. We see that, up to the fluctuations, all curves join nicely on the same curve. Figure 8 shows the variation of the solid fraction ν = V p /V as a function of ε 1 for all the decompositions, where V p is the volume (area in 2D) occupied by the particles. The solid fraction decreases first from ν 0 0.84 to 0.825. It is remarkable that, during this phase the solid fraction is rigorously identical for all the samples. At larger strains, the solid fraction decreases much more slowly and, up to the fluctuations, saturates on the same curve. Indeed, at larger strains, dilation is localized within shear bands appearing and vanishing throughout the system underlying the saturation of ν. This is well illustrated in Figure 9 where two maps of the particle velocities are shown at ε 1 = 0.2 for S0 and S3 cases. We see clearly that the topology of the shear band are slightly different even if, on average, the solid fraction is identical. In fact, localization phenomena leads to multiple possible physical solutions, and it has been already been exemplified that the formation of the shear band depends on the details of the numerical parameter of the simulations (time step, solver, number of iterations...) [START_REF] Renouf | Conjugate gradient type algorithms for frictional multi-contact problems: Applications to granular materials[END_REF].

In this section we have shown that the macroscopic response of a sheared granular material is independent of the chosen number of subdomains. Nevertheless, a granular material is a typical example of multi-scale material: the macroscopic behavior results from the average properties of a collection of interacting particles through the contact network. This is clearly illustrated in the case of elongated particles where the residual shear strength increases linearly with elongation [START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF], whereas for angular or non-convex particle shapes, the residual shear strength increases first and saturate as the level of angularity or non-convexity of the particles is increased [START_REF] Azéma | The effect of the angularity of particles on the mechanical response of sheared granular packings[END_REF][START_REF] Saint-Cyr | Rheology of granular materials composed of nonconvex particles[END_REF]. An other surprising effect is that the residual shear strength is independent of the polydispersity [START_REF] Voivret | Multiscale force networks in highly polydisperse granular media[END_REF]. This wide variety of behaviors finds its origins at the scale of the particle and contact properties. Thus, we also need to test the robustness of the domain decomposition solver in terms of the granular microstructure.

Micromechanical analysis

The granular microstructure (granular texture), i.e. the spatial organization of the particles and their contacts, is basically controlled by steric exclusions between the particles and force balance conditions [START_REF] Troadec | Model for granular texture with steric exclusions[END_REF]. The strong inhomogeneity of contact forces is a well known feature of granular media. Figure 10 shows the contact forces for S2 and S4 cases at the residual state. For the same level of strain, the force-carrying backbones are different, even if the global inhomogeneity seems to be preserved. Of course, this is due to the fact that the resolution of the contact forces is performed domain by domain, and to the plurality of local solutions for frictional granular media. Nevertheless, this anisotropic structure, generally at the origin of the shear strength of granular media, can be described more rigorously in terms of various statistical descriptors pertaining to the forcebearing network of particles. In the following, we consider two aspects of the microstructure: (i) the angular average of the contact and force orientations, and (ii) the normal force distributions.

Contact and force anisotropy

A common approach is to consider the probability distribution of the contact normals n, which is usually nonuniform. As shown in Figure 1, for the 2D case the unit vector n is described by a single angle θ. The probability density function P (θ) of contact normals provides detailed statistical information about the texture. In the same way, expressing the force vector in the local contact frame (n, t), where t is an orthonormal unit vector oriented along the tangential direction, we can compute the angular distributions f n (θ) and f t (θ) of normal and tangential forces, respectively. The above three functions describe the general state of the packing. Under shearing, the packing organizes itself into a state where these functions are well approximated with their lowest-order Fourier expansion [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF][START_REF] Kruyt | Kinematic and static assumptions for homogenization in micromechanics of granular materials[END_REF][START_REF] Radjaï | Key features of granular plasticity[END_REF][START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF]:

   P θ (θ) = 1 2π {1 + a c cos 2(θ -θ c )} f n (θ) = f {1 + a f n cos 2(θ -θ f n )} f t (θ) = f a f t sin 2(θ -θ f t ) (23)
where a c , a f n , and a f t are the anisotropy parameters and θ c , θ f n , and θ f t represent the corresponding privileged directions. In a sheared state the privileged directions tend to follow the principal stress direction (i.e. . The multiplicity is: 1 for a light gray particle, 2 for a medium gray particle and 3 for a black particle.

θ c = θ f n = θ f t = θ σ ).
In practice the values of all anisotropy parameters can be computed from generalized fabric and force tensors presented in [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF][START_REF] Azéma | Stress-strain behavior and geometrical properties of packings of elongated particles[END_REF]. Figure 11 shows the variations of all these anisotropy parameters as functions of ε 1 for all the domain decompositions. Up to fluctuations, all the curves join also nicely on the same curve. We see that a c follows the same trend as the shear strength, increases first to a maximum value equal to 0.35 and then declines to a plateau at 0.26. In contrast, a f n decreases from 0.38 to a constant value 0.24, whereas a f t remains nearly constant with the strain.

These anisotropy parameters are very relevant to the analysis of the granular microstructure because they can bring to light the geometrical and mechanical origins of the shear strength. Indeed, it can be shown that the general expression of the stress tensor (21) leads to the following "stress-force-fabric" relationship (a term coined for the first time by Rothenburg and Bathurst in [START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF]):

sin ϕ 1 2 (a c + a f n + a f t ), (24)
where the cross products between the anisotropy parameters have been neglected. It is very important to test the validity of relation [START_REF] Jean | The non-smooth contact dynamics method[END_REF] in the context of numerical simulation by DDM because this equation reveals an explicit link between microscopic and macroscopic scales in granular media. We plot in Figure 12 the variation of the shear strength together with the harmonic fit of equation [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. We see that data are in quantitative good agreement with this harmonic approximation which is not affected by the number of subdomains.

Figure 12: Friction angle sin ϕ (symbol) together with the harmonic approximation (line) as functions of the vertical strain ε 1 for all the tested decompositions.

Force distribution

Force transmission has been investigated by experiments and numerical simulations for disks, elongated, polygonal and non-convex particles in 2D as well as for spherical, cylindrical and polyhedral particles in 3D [START_REF] Radjaï | Force distributions in dense two dimensional granular systems[END_REF][START_REF] Coppersmith | Model for force fluctuations in bead packs[END_REF][START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF][START_REF] Richefeu | Force distribution in cohesive and non cohesive granular media[END_REF][START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF]. The probability density function (pdf) of normal forces is characterized by two features that are specific to granular media: (i) the pdf is roughly a decreasing exponential function for forces above the mean value, (ii) in the range of weak forces below the mean value, the pdf does not decline to zero with the force. The relative scattering of data reported by different authors for weak forces shows the sensitivity of the pdf within this range to the microstructure details.

The probability density function (pdf) of normal forces normalized by the mean normal force f n is shown in Figure 13 in log-linear and log-log scales at large strains (the data are cumulated from several snapshots in the residual state) for all the simulations. As usually observed in the literature [START_REF] Radjaï | Force distributions in dense two dimensional granular systems[END_REF][START_REF] Coppersmith | Model for force fluctuations in bead packs[END_REF][START_REF] Mueth | Force distribution in a granular medium[END_REF][START_REF] Antony | Evolution of force distribution in three-dimensional granular media[END_REF][START_REF] Richefeu | Force distribution in cohesive and non cohesive granular media[END_REF][START_REF] Azéma | Force transmission in a packing of pentagonal particles[END_REF][START_REF] Azéma | Quasistatic rheology, force transmission and fabric properties of a packing of irregular polyhedral particles[END_REF], the number of forces above the mean value f n falls off exponentially whereas the number of forces below the mean value varies as a power-law:

P (f n ) ∝ e -α(1-fn/ fn ) for f n > f n (f n / f n ) β for f n < f n (25)
where α and β are the exponents. As we can observe, the pdf of the forces in each samples collapse to the same curve given precisely by [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF]. This shows that the inhomogenous character of the force distribution chain in granular media is not affected by the number of subdomains used for the solver. The genuine organization of contact forces in granular media was first analyzed by Radjaï et al. by means of contact dynamics simulations for packings of circular and spherical particles [START_REF] Radjaï | Bimodal character of stress transmission in granular packings[END_REF]. The most important result was that the contact network can be decomposed unambiguously into two subnetworks named "weak" and "strong" networks with complementary mechanical properties. More precisely, stronger forces chains are propped by large number of weak contacts, so that the shear stress is almost totally sustained by the strong contact network. This is well illustrated also in Figure 10 where strong forces are plotted in black whereas weak forces are plotted in red. We see that strong forces are mainly vertical (along the principal shear direction) whereas weak forces are, in average, perpendicular to the direction of shear. Data are also in good qualitative agreement with this feature, without much influence of the number of subdomains.

Time consuming analysis of the NSCDD Sequential Multidomain implementation

The issue of this section is to estimate the CPU time gain that we can expect from a multiprocessing implementation of the NSCDD algorithm on a distributed memory architecture. We present at first the implementation and the parameters chosen before to an analyse of the numerical tendencies of the proposed method.

Sequential Multidomain implementation and chosen parameters

A Sequential Multidomain implementation of NSCDD algorithm has been performed on the LMGC90 software to study the influence of a domain decomposition on the biaxial test presented above. To do so, the sequential LMGC90 database has been duplicated, according to the number of subdomains considered, to mimic the behavior of a multiprocessing environment. This approach allows to separate the topic of physical validity of the solution given by the proposed domain decomposition method from the technical aspects of MPI implementation.

For the simulations we selected n gs = 1 (cf. Algorithm 1), which means that one NLGS iteration in subdomains is always followed by an interface resolution, consistently with the study of the influence of n gs parameter on the convergence of the NSCDD algorithm reported in [START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF].

The cumulative elapsed time and timers of the main steps of the NSCDD algorithm for samples S0, S2 and S4 is given in Table 1, for the Biaxial test presented above, and performed over 250 10 3 time steps. The various stages are classified as:

• White: generic stages of NSCD algorithm which may be parallelized,

• Light-gray: stages introduced by the domain decomposition which may be parallelized (intermediary routines between generic stages of NSCD algorithm and specific stages of the NSCDD one).

• Dark-gray: stages of the NSCDD algorithm which must be done sequentially (if a slave/master communication scheme is presupposed) and implies exchanges between processors.

Analysis of the main referenced stages

"Domain partitioning and rough detection" (cf Table 1). As presented above, the proposed domain partitioning leans on contact distribution among the subdomains. More precisely, the considered contacts are those roughly selected according to a box method referenced as the "rough detection". The elapsed time for each samples is quite similar, a small increase is observed only for S4. In our simulations this phase is performed every 10 time steps leading to a very small contribution to the overall computation time. "Fine detection" and "Update positions". The increase in the elapsed time related to those stages arises from the number of particles in samples, which increases with the number of subdomains (n sd ) due to the duplication of interface grains, as illustrated in Table 2.

"NLGS preprocessing". This routine stores every extra-diagonal block matrix of the Delassus operator (W αβ , α = β). Its running time decreases according to the number of subdomains. Indeed the duplication of the interface grains in the neighboring subdomains reduces the number of the adjacent contacts (and so the number of W αβ , cf. Table 2). This phenomenon also explains the similar behavior of "NLGS iterations" stage.

"Compute vd E (F Γ )" and "Compute A ΓE VE ". As expected, the elapsed time increases regularly according to the size of the interface.

"Check Convergence". A similar elapsed time is observed for each samples. The numerical monitoring shows that the time consuming stages may be parallelized whereas the sequential stages requires at most 5% of the total CPU time in our study. Moreover, even for a Sequential Multidomain implementation, the total elapsed time may be reduced when using several subdomains in comparison with a single subdomain, in spite of the increase of the interface size (in terms of unknowns and equations). This is due to the simultaneous decrease of W (from 418 10 3 to 303 10 3 W αβ , cf. Table 2). However, the expected gain from MPI implementation may be quite different because of the potentially expensive exchanges between processors.

Conclusion

The present work gives a new illustration of the ability to use a domain decomposition method coupled with the nonsmooth contact dynamics approach for dealing with large-scale dense granulates. The proposed approach is as close as possible to the standard nonoverlapping DDM for large-scale linear problems, more precisely the FETI approach. As the interfaces are made of grains, the features of the interface matrix has been systematically studied, for example when a grain belongs to more than two subdomains. A mechanical analysis of a biaxial test exemplifies the relevancy of the results in spite of the chaotic behavior of such a system with a large multiplicity of solutions. The solutions may depend locally on the substructuring procedure but the global behavior of the granular medium is preserved. We rediscover the sensitivity, with respect to the discretization, of a ductile material involving localization effects such as shear bands, of plastic nature, that it is modeled by finite elements or by discrete elements. However the discrete elements are not determined by a discretization process but imposed by the microstructure. Moreover the forthcoming behavior is not strongly oriented by the early localization because of the appearance and the vanishing of multiple different shear bands in dry granulates.

The numerical efficiency, especially the scalability, is recovered if a single NLGS iteration is performed in each subdomain in the first stage of the algorithm [START_REF] Iceta | Mixed versus impulse-oriented domain decomposition method for granular dynamics[END_REF]. It is not necessary to iterate many times in the first stage because the second stage, characterized by the quasi-diagonal interface matrix, does not transfer long-distance correlation through subdomains. Likewise it is not possible to improve the convergence by adding a coarse problem as in classical computational structural mechanics based on the finite element method. Before extracting a macro-homogenized model from the interface problem it is neces-sary to enrich this interface problem. Such an approach has been developed in [START_REF] Alart | A nonlinear domain decomposition formulation with application to granular dynamics[END_REF] from a theoretical and semi-analytical point of view. As a conclusion of this study the convergence of the so obtained algorithm does not seem to be significantly accelerated whereas the computational cost of the second solution stage strongly increases. Consequently we propose now to develop, not a nonsmooth solver on a single time step but a multiscale time integration scheme over a time interval for granular systems. The principle would be to combine an explicit linear prediction of the interface forces and an implicit correction of the contact impulses inside the subdomains.

Finally the main drawback of the NSCD approach is the possible indetermination of the contact impulses generated by severe kinematic constraints, especially for dense granular systems. This is conveyed in the singularity of the Delassus matrix whose null space represents the self-equilibrated impulse networks. To overcome this indetermination, from an algorithmic viewpoint, represents an important challenge. To conceive a time integration scheme during a process requires to tackle this topic.

Figure 2 :

 2 Figure 2: Geometrical partitioning of the discrete domain and duplication of the interface grains.

Figure 3 :

 3 Figure 3: Illustrations of the proposed domain partitioning technique: four grains having a multiplicity of 2 (a) (cf. Figure 6(b)), four grains having a multiplicity of 2 and one having a multiplicity of 3 (cf. Figure 10(b)). Contacts are colored according to the subdomain they belong to. The grains having a multiplicity of 2 are hatched, the grain having a multiplicity of 3 is crossedhatched.

Figure 4 :

 4 Figure 4: Gluing conditions between interface grains. (a) m i = 2 and ni l = m i -1 = 1 2 m i (m i -1) = 1; (b) m i = 3 and n i l = m i = 1 2 m i (m i -1) = 3; (c) m i =3 and n i l = m i -1; (d) m i = 4 and n i l = m i ; (e) m i = 4 and n i l = 1 2 m i (m i -1) = 6; (f) m i = 4 and n i l = m i -1. (For convenience the grains have been split even if they are not geometrically cut.)

  Figure 4: Gluing conditions between interface grains. (a) m i = 2 and ni l = m i -1 = 1 2 m i (m i -1) = 1; (b) m i = 3 and n i l = m i = 1 2 m i (m i -1) = 3; (c) m i =3 and n i l = m i -1; (d) m i = 4 and n i l = m i ; (e) m i = 4 and n i l = 1 2 m i (m i -1) = 6; (f) m i = 4 and n i l = m i -1. (For convenience the grains have been split even if they are not geometrically cut.)

Figure 5 :

 5 Figure 5: Boundary conditions for (a) isotropic and (b) biaxial compactions.

Figure 6 :

 6 Figure 6: Examples of decomposition at the initial state (a) and zoom at the intersection of the four subdomains of case S4 (b). The multiplicity is (b): 1 for a gray particle and 2 for a black particle.

Figure 7 :

 7 Figure 7: Internal angle of friction sin ϕ as a function of the vertical deformation ε 1 .

Figure 8 :

 8 Figure 8: Solid fraction ν as a function of the vertical deformation ε 1 .

S0Figure 9 :

 9 Figure 9: Maps of the normalized angular velocities of the particle for S0 (a) and S3 (b) cases, at ε 1 = 0.2.

Figure 10 :

 10 Figure 10: Maps of force chains in a portion of S2 (a) and S4 (b) samples, at ε 1 = 0.2; (c) is a zoom on map at the subdomains intersection (b). Line thickness is proportional to the normal force. Strong forces (f n > f n ) are in black and weak forces (f n < f n ) in gray (see text). The multiplicity is: 1 for a light gray particle, 2 for a medium gray particle and 3 for a black particle.

Figure 11 :

 11 Figure 11: Contact anisotropy (a) and normal and tangential force anisotropies (b) as functions of the vertical deformation ε 1 for all the tested decompositions.

Figure 13 :

 13 Figure 13: Probability distribution function of normal forces f n , normalized by the average normal force f n in log-linear (a) and log-log (b) plots for all the tested decompositions.

3 Table 2 :

 32 Total number of NSCDD iterations, mean number of interface particles and mean number of extra-diagonal block matrices of the Delassus operator (W αβ ) over the 250 10 3 time steps of the processes related to samples S0, S1, S2, S3 and S4.

  are verified. Algorithm 1 NonSmooth Contact Domain Decomposition (NSCDD) for i = 1, . . . , N do contact detection (eventually parallelized) and possible new decomposition of the domain initialize unknowns at time t i : (r E , v E , F Γ ) while (convergence criterion not satisfied) do In parallel for E = 1, . . . , n sd Disassemble interface impulses F Γ into local impulses F E Compute the "free" velocity V d E and vd E ,v E ) with n GS non-linear Gauss-Seidel iterations on:

E

Compute (r

  vE ) Compute RE and correct the velocity on interface grains to get A ΓE VE In sequential (may be partially parallelized) Compute ∆F Γ as: X∆F Γ =

	impulses F Γ	n sd E=1 A ΓE VE and update interface
	end while
	Update grain positions in parallel
	end for

Table 1 :

 1 Percentage and absolute elapsed time in the main stages of the NSCDD algorithm related to samples S0, S2 and S4.

	Main			S0	S2		S4
	stages/Samples
			%	Elapsed	%	Elapsed	%	Elapsed
				time		time		time
				Preprocessing
	Domain	parti-	2 %	4 10 3 s	2 %	4 10 3 s	2 %	5 10 3 s
	tioning and rough
	detection
	Fine detection	5 %	14 10 3 s	6 %	14 10 3 s	8 %	19 10 3 s
	NLGS	prepro-	32 %	83 10 3 s	25 %	57 10 3 s	16 %	41 10 3 s
	cessing
				NSCDD iterations
	Compute vd E (F Γ ) NLGS iterations Compute A ΓE VE Interface problem	0 % 50 % 131 10 3 s 0 10 3 s 0 % 0 10 3 s 0 % 0 10 3 s	4 % 42 % 101 10 3 s 40 % 100 10 3 s 9 10 3 s 7 % 17 10 3 s 9 % 21 10 3 s 13 % 33 10 3 s 1 % 2 10 3 s 2 % 4 10 3 s
	Check	conver-	1 %	2 10 3 s	1 %	1 10 3 s	1 %	1 10 3 s
	gence
				Updates and outputs
	Update positions	9 %	23 10 3 s	9 %	21 10 3 s	10 %	26 10 3 s
	Write files		2 %	6 10 3 s	2 %	5 10 3 s	2 %	5 10 3 s
	Total		100 % 262 10 3 s 100 % 235 10 3 s 100 % 251 10 3 s

Acknowledgements. We warmly thank F. Radjaï for useful and stimulating discussions on the physics and micromechanics of granular materials and F. Dubois for help with numerical developments. This work was partly supported by OSEO, FEDER and the region of Languedoc-Roussillon (Degrip project).