S D Nakatsutsumi

S Baton

Hulin

X Vaisseau

A Morace

M Touati

M Nakatsutsumi

S Hulin

Ph Nicolaï

R Nuter

D Batani

F N Beg

J Breil

R Fedosejevs

J.-L Feugeas

P Forestier-Colleoni

C Fourment

S Fujioka

L Giuffrida

S Kerr

H S Mclean

H Sawada

V T Tikhonchuk

J J Santos

Collimated propagation of fast electron beams accelerated by high-contrast laser pulses in highly-resistive shocked-carbon

Keywords:

Collimated propagation of fast electron beams

accelerated by high-contrast laser pulses in highly-resistive shocked-carbon X.

The transport of intense currents of relativistic electron beams (REB) driven by intense lasers into dense matter is a key topic for applications such as secondary particle acceleration, x-ray and gamma-ray generation, isochoric heating to temperatures relevant in the context of planetology or astrophysics, and laser-driven thermonuclear fusion. The fast ignition scheme (FI) [1] proposes to ignite a pre-compressed capsule of deuterium-tritium fuel by the rapid energy-deposition from a REB. The key figure for FI-feasibility is the energy coupling of the ignitor laser pulse and the dense fuel-core: 20 kJ of 1-2 MeV electrons should be delivered into a 20 µm-radius hot spot [2]. The laser to REB energy-transfer can, in principle, be up to ∼ 40% [3][4][5][6], yet the source angular distribution of REB momenta is intrinsically broad [7][8][9]. The use of a hollow cone re-entrant into the capsule provides a clear pathway for the ignitor laser till the cone tip [10,11], moving the REB-source as close as 100 µm from the high-density core [12]. Yet, the challenge persists of confining the REB propagation within a small radius, while limiting both resistive and collisional energy losses over the compressed plasma [13][14][15] as well as REB-source perturbations either by cone-tip ablation by the laser intensity-pedestal [16][17][18] or by tip disruption due to the high-pressure of the asymmetrically imploded plasma [19].

No methods have been devised yet to control the REB divergence directly at the source, but a few strategies to magnetically collimate REB propagation have been proposed [20]. Magnetic-field (B-field) strengths in the range of kT, either self-generated or imposed, are needed to produce a sufficient effect on the MeV electrons. Selfgenerated B-fields arise from a resistive return current at high REB current densities, j h 10 11 A.cm -2 [21]. Target structuring with different materials allows to obtain radially-convergent resistivity gradients. Such fields produced an observable collimation effect on REBtransport in both planar [22] and cylindrical [23] experimental designs. Numerical simulations explored the effect in even more complex cone-target structures [24][25][26]. However it is unlikely that they will survive the implosion and ignitor-laser irradiation. Other strategies explored experimentally are the REB interaction with radially-converging resistivity-gradients by cylindrical target-compression [27] or a precursor less-intense laser pulse generating a seed B-field [28]. Yet, in all the cases an efficient collimation, with small electron losses, has yet to be demonstrated.

Other studies pointed out that materials of higher atomic number, such as Au, of higher heat capacity and ionization level, yield higher resistivity-gradients over the REB time-scale, favoring its collimation through selfgenerated B-fields [29,30]. Even though, such dense materials are prone to increase the collisional diffusion and energy losses of the beam electrons. A recent study pointed out that not so dense vitreous-carbon is highly resistive over a broader temperature range up to ∼ 100 eV, therefore facilitating rapid B-field growth [31] even under the conditions of very intense currents. Our experiment was focused on exploring REB transport in samples of such high-resistivity, while of relatively lowdensity, under conditions linked to the FI scenario: REB generated in cone targets and propagating into shockcompressed warm dense matter.

The experiment was carried out at the LULI2000 facility (Ecole Polytechnique, France). The targets consisted of a 10 µm-thick, 450 µm-long hollow Cu-cone buried in a 500 µm-thick vitreous-C layer of ρ 0 = 2 g.cm -3 -density [Fig. 1 a)]. A 10 µm Cu-tracer was located at a nominal distance L 0 = 50 µm away from the outer cone tip, followed by a 15 µm plastic (CH) layer acting as an ablator. We used a dual laser beam setup: a compression shock was driven at the CH-side of the target by a frequency doubled long pulse (LP) laser (λ LP = 0.53 µm), with a τ LP = 3 ns square temporal profile (∼ 100 ps rise-time). It was focused by means of a random phase plate onto a 340 µm flat-top spot containing ∼ 61% of the 496 ± 50 J on-target energy, yielding an intensity of 1.3 × 10 14 W.cm -2 . The λ SP = 1.053 µm, τ SP = 1.3 ps full-width-at-half-maximum (FWHM) short pulse (SP) laser was focused on the inner cone tip, of 50 µm diameter, in order to generate a REB propagating against the shock. A plasma mirror was used to direct the SP laser along the cone axis and to reduce its pedestal due to the amplified spontaneous emission (ASE) down to ∼ 10 10 W.cm -2 intensity. We obtained Gaussian-like focal spots of 10 µm FWHM containing ∼ 60% of the on-target 33 J SP-laser energy, yielding a peak intensity of 2 × 10 19 W.cm -2 . For each shot, the delay ∆t between the two laser pulses (SP peak intensity in respect to the LP half-rise time), was set according to the required compression state of the target with a precision ≤ 50 ps. The chosen ∆t were determined by a simulation of the shock formation and compression using the 2D radiative-hydrodynamic Lagrangian code chic [32], benchmarked by shock-breakout measurements on Al foils. The description of the shocked cone-targets' hydrodynamics was further benchmarked in a previous experiment by streaked optical pyrometry of the inner-cone tip and target side-on radiography using Cu K-shell ionization fluorescence as backlighter [33].

Fast electron generation and transport in the target were investigated from side-on imaging of the Cu-K α ionization fluorescence from both the cone and the rear-side tracer using a spherically-bent Bragg crystal [34], providing simultaneous access to the REB spatial distribution at the source and after crossing the C-layer. Fig. 1 b) shows sample data as a function of ∆t, from ∆t = 0 ns (REB injection into an unperturbed solid-cold vitreous-C target), up to ∆t = 4 ns (1 ns after shock-breakout time at the inner cone-tip surface; REB generation in an expanding plasma filling the cone volume). The use of the LP laser also for the ∆t = 0 -case allowed the formation of a long CH ablation-plasma that prevented fast electrons from recirculating after their first transit through the Cu layers [13,35]. The simulated density and temperature maps at the corresponding delays are presented on Fig. 1 c). At ∆t = 3 ns, just before shock breakout, the C-layer density is increased by a factor ∼ 2.7 over its full thickness while its temperature reaches ∼ 3 eV. Looking only at the cone emission for times before shock-breakout, it mainly originates from the tip, highlighting the efficiently localized energy coupling of the high-contrast SP laser.

The emission from the rear side tracer is also clearly visible for the two earlier times, showing that a significant fraction of the REB energy was transported into the target depth. For ∆t = 2.5 ns, the rear side tracer is being pushed towards the cone tip by the compression shock, as seen on the hydro simulation. The tracer fluorescence yield is higher with respect to the ∆t = 0 image, in relation to a reduction of the fast electrons energy-loss due to the resistive effect (collisional losses are kept the same for samples of identical areal mass) [13]. An important feature is that the tracer emission size is of the same order as that from the cone tip, consistent with its 50 µm-diameter, despite the several tens of microns propagation distances, as it is shown on Fig. 2 a) by the measured diameters of the tracer-fluorescence (blue full triangles). This data set demonstrates a collimated propagation of the REB over the carbon layer, either entirely cold-solid (∆t = 0 ns), or partially (∆t = 1 ns, ∆t = 2 ns) or entirely in warm-dense state (∆t = 2.5 ns, when the shock front is entering the cone tip). Such collimation was confirmed by doubling the cone tip-Cu tracer distance to L 0 = 100 µm, resulting in a Cu-K α diameter of ∼ 60 µm for the unperturbed target [blue full square on Fig. 2 a)]. Given the longer propagation distance, this is still fairly close to the cone-tip diameter. For ∆t ≥ 3 ns the emission from the Cu tracer can no longer be distinguished from the emission of the cone tip as the distance between them, as predicted by the hydro simulation and target radiography results [33], is smaller than the ∼ 10 µm-resolution of the Cu-K α imager. The tip fluorescence yield reaches its highest value at ∆t = 3 ns due to the peak Cu-density under the shock.

The fluorescence patterns are significantly different when injecting fast electrons after the shock breakout inside the cone (∆t > 3 ns). Fast electrons were then generated in an expanding Cu plasma progressively filling the volume of the cone. As evidenced from the image at ∆t = 4 ns (note the different spatial scale), this dramatically increased electron divergence, their dilution into the plasma-filled cone volume and the length of their interaction with the cone walls. Consistently, the Cu-K α yield decreases. A steep decrease in the global energy coupling to the target depth after shock breakout is corroborated from spatially integrated measurements of the Cu-K α yield as a function of the delay between the two laser pulses, using a potassium acid phthalate (KAP) conical crystal at the 5th diffraction order (not presented) [36].

To further understand the physics at play, in particular the mechanisms of REB collimation/dispersion, the experimental data were compared to a suite of numerical simulations. The REB source was computed with the 1D particle-in-cell (PIC) code ocean [37] modeling the interaction between the SP beam and the Cu generation layer. Collisions were not taken into account since the goal of these PIC simulations was only to estimate the energy spectra of the REB source. The SP temporal profile was a 1.5 ps FWHM Gaussian with a peak intensity of 10 19 W.cm -2 . We considered four different compression delays corresponding to those shown on Fig. 1, the ∆t = 0 and 2.5 ns cases corresponding to an uncompressed cone tip. In the cases before shock breakout time the pre-plasma gradient length was set to 0.2 µm, due to the high-contrast of the laser pulse. The REB source characteristics were extracted after crossing the same quantity of matter in each simulation (same 4.5 mg.cm -2 areal mass), corresponding to a 5 µm depth into the cone Cu layer for the unperturbed cone-tip cases. The REB-source energy distributions, averaged over the 5.3 ps-duration of the simulations, are presented in Fig. 3 a). Since the distributions are fairly similar before shock breakout time, only two cases are shown in the figure .  The REB distribution exhibits a 1.5 MeV temperature component when generated on a several microns length plasma ramp (after shock breakout time), which is not present for interactions at the early times. The REB mean energy rises from ∼ 20 keV to ∼ 170 keV before and after shock breakout time, respectively. The predicted total number of accelerated electrons is higher in the last case by virtue of a higher laser energy absorption over the long scale-length plasma, yet the energy coupling into the dense target regions is reduced due to the fast electrons dilution into the plasma ramp and cone walls, as inferred from the Cu-K α images described above.

The REB transport was simulated with a reduced ki- netic code solving the linearized Vlasov-Fokker-Planck equation with an M1 approximation [38]. The 2D hydro simulation results, described above, were used as initial density and temperature profiles at the corresponding REB injection times. The SP laser energy was 33 J and laser-to-fast electron conversion efficiency η L→e = 20% before shock breakout time (in agreement with Davies' scaling law [39]). The fast electron angular distribution was described by a mean propagation angle θ 0 ∼ 25 • tanh [r/r b ], with r the transverse coordinate and r b = 8.5 µm, and a mean dispersion angle ∆θ = 25 • , in agreement with Green's scaling law [8]. Transport simulations were ran up to 4 ns with space and energy grids of 1 µm and 10 keV, respectively. Carbon resistivity at solid-density was inferred from quantum molecular dynamics calculations [42] and extended to high temperatures 500 eV using the Spitzer model [43] [Fig. 3 b) (blue solid line)]. As the compression was relativelly small, the same expressions were used for the resistivity of shocked carbon. The Cu and Al resistivities were calculated with the Drude model, η = m e ν/e 2 n e , with n e the electron density and ν the collision frequency given by a harmonic mean between the Lee & More model [40] for hot plasma temperatures and the Eidmann-Chimier model [41] for lower solid-liquid temperatures. The Cu-K α emission was characterized using a model based on the Hombourger's cross sections [44] and the computa-tion model proposed in [45]. The careful description of the REB source, the target hydrodynamics and the transient resistivity evolution resulted in a good agreement between the size of the synthetic [blue empty symbols in Fig. 2 a)] and the measured Cu-K α spots at the different positions of the rear side Cu-tracer. This confirms the collimation of fast electrons in both compressed (warm and dense) [Fig. 3 c), L 0 = 50 µm-case at ∆t = 3 ns] and cold-solid carbon layers [Fig. 3 e), un-driven L 0 = 50 µmcase]. This behavior can be explained by looking at the corresponding maps of azimuthal self-generated B-field [Fig. 3 d) and f)]. The collimating field of 500 T seen in the high-resistivity C-layer mainly originates from the curl of the REB current density, largely compensating the decollimating effect due to resistivity gradients arising from the rapid REB energy deposition along the cone axis.

In order to confirm the crucial role played by the vitreous-C resistivity, we performed two additional simulation tests consisting i) in suppressing the bump seen at ∼ 6 -70 eV temperature range in the vitreous-C resistivity [black dotted line on Fig. 3 b), labeled CH-like], and ii) in using instead the resistivity of Al (a conductive material) [red solid line on Fig. 3 b)]. The suppression of the C resistivity bump alone (first test) already reduces by a factor of 2 the amplitude of the self-generated B-field, yielding a much less collimated REB propagation [Fig. 3 g) and h)]. In the second test, the B-field is even weaker and the REB is significantly divergent [Fig. 3 i) and j)]. These two tests yield synthetic Cu-K α spot sizes [respectively the black open circle and the red semi-square for L 0 = 50 µm, and the black and red half-filled squares for L 0 = 100 µm in Fig. 2 a)], larger than those obtained for vitreous-C in both experiments and simulations.

The enhancement of the REB energy density flux associated with the collimated transport in vitreous-C compared to other materials is however mitigated by an increase of the resistive energy losses due to the vitreous-C higher resistivity: These increase by ∼ 25% over 50 µm of C when compared to the less resistive Al case. It is nonetheless worth to point out that the additional number of electrons stopped in vitreous-C consists of the less energetic component of the REB spectrum, which would in any case be stopped before reaching the compressed DT-core in a FI fusion target.

In all cases is possible to identify an enhanced radial diffusion of the lower-energy electrons when reaching the Cu-tracer (denser than the C layers) as well as the development of surface B-fields at the C-Cu interfaces attributable to the corresponding resistivity discontinuities. Yet, these features do not interfere with the conclusions concerning the role of resistivity in the electron transport through the C-layers (shocked or undriven) and from the comparison of transport in C, CH or Al.

In conclusion, we characterized the spatial distribution of relativistic electron beams generated by coupling a high-intensity laser pulse into a vitreous carbon target in both cold-solid and warm-dense conditions. The use of a high-contrast laser pulse allowed its efficient coupling to the tip of a cone drilled into the target before shock breakout time. A strong dilution of the fast electron source was observed when firing the laser pulse after shock breakout, confirming the deleterious effect on the REB energy-density flux due to plasma-filling of the cone [16][17][18][19]. More importantly, a clear collimation of the fast electron beam over several tens of microns distances was experimentally achieved in both cold-solid and warm-dense vitreous C propagation layers, and reproduced by numerical transport simulations. This remarkable behavior is attributed to the high-resistivity of vitreous C, generating several hundreds of Teslas collimating B-fields, even under compressed matter conditions and over a wide range of temperatures up to 100 eV. An additional set of transport simulations using different materials confirmed the importance of the C high resistivity in this range of temperatures. As a counterpart, the high resistivity of vitreous-C was responsible for additional resistive energy losses when compared to the less resistive cases. Noteworthily, the collisional losses are kept low by virtue of the C low-density, an advantage compared to other studies proposing high-resistivity propagation media using high-Z materials of larger heat capacity and higher ionization level [29,30]. The ex-plored behavior of vitreous-C can be of considerable importance for an advanced design of a FI target by inserting a low-density and high-resistivity plug to a (thinner, ∼ µm) cone-tip of a denser material: The light and resistive plug, placed close to the REB source, would be thick enough to mitigate the REB natural divergence, but thin enough to moderate the subsequent rise of resistive energy losses.

We gratefully acknowledge the support of the LULI pico 2000 staff during the experimental run. This work was performed through the funding from the French National Agency for Research (ANR) and the Competitiveness Cluster Alpha -Route des Lasers, Project TERRE ANR-2011-BS04-014. The work was carried out in the framework of the Investments for the future Programme IdEx Bordeaux LAPHIA (ANR-10-IDEX-03-02) and the EUROfusion Consortium, and has received funding from the European Union's Horizon 2020 research and innovation program under Grant No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

ρρFIG. 1 .FIG. 2 .

 12 FIG. 1. (color online) a) Experimental setup: target design, geometry of laser irradiation and Cu-Kα sight axis (top view). b) Experimental side-on images of the Cu-Kα fluorescence, attributable to the time-integrated REB spatial distributions, as a function of ∆t. Note the different signal-intensity scales. Dashed lines indicate the initial position of the cone walls and the front surface of the Cu tracer. The dotted lines indicate the shock-front positions. c) Density (top) and temperature (bottom) maps as a function of the delay between the laser pulses ∆t, obtained with 2D rad-hydro simulations. The maps' aspect ratio has been modified for clarity.

FIG. 3 .

 3 FIG. 3. (color online) a) Fast electron source energy distribution extracted from PIC simulations before shock breakout time (black solid line) and 1 ns after shock breakout time (blue dashed line). b) Electric resistivity functions used in transport simulations: nominal resistivity of vitreous-C, from [42] (blue solid line), CH-like resistivity obtained by suppressing the bump in the 6 -70 eV range (black dotted line) and Al resistivity (red solid line). On the top, c), e), g) and i): fast electron density maps extracted from transport simulations at 2.7 ps, the time of the peak azimuthal B-field [maps on the bottom, d), f ), h) and j); note the different color scales]. Simulations with vitreous-C resistivity: c) and d) L0 = 50 µm-case just before shock breakout time (∆t = 3 ns, the C-layer is fully compressed), e) and f ) un-driven L0 = 50 µm-case. g) and h) Un-driven L0 = 50 µm-case using the CH-like resistivity. i) and j) Un-driven L0 = 50 µm-case using the Al resistivity. The z = 0-coordinate corresponds to the inner cone-tip surface.


Malik Muhammad Haris

Domenico Borzacchiello

José Vicente Aguado

Advanced parametric space-frequency separated representations in structural dynamics: A harmonic-modal hybrid approach

come   L'archive ouverte pluridisciplinaire


François Boulier

The Management of Parameters in the MAPLE DifferentialAlgebra Package

This technical report contains a description of the management of the parameters in the MAPLE DifferentialAlgebra package and, in particular, in the RosenfeldGroebner function.

Need of Algebraic Specifications

The Rosenfeld-Groebner algorithm is an elimination method in differential algebra [START_REF] Fels | Differential Algebra[END_REF][START_REF] Robert | Differential Algebra and Algebraic Groups[END_REF]. Many different versions of it were written [START_REF] Boulier | Étude et implantation de quelques algorithmes en algèbre différentielle[END_REF][START_REF] Boulier | Representation for the radical of a finitely generated differential ideal[END_REF][START_REF] Boulier | Réécriture algébrique dans les systèmes d'équations différentielles polynomiales en vue d'applications dans les Sciences du Vivant[END_REF][START_REF] Boulier | Computing representations for radicals of finitely generated differential ideals[END_REF] and implemented: a first version in the diffalg MAPLE package [START_REF] Boulier | diffalg. Package of MapleSoft MAPLE standard library from MAPLE V to MAPLE 13[END_REF], a second version wihing the BLAD libraries [START_REF] Boulier | The BLAD libraries[END_REF], available in MAPLE through the DifferentialAlgebra package [START_REF] Boulier | DifferentialAlgebra. Package of MapleSoft MAPLE standard library since MAPLE 14[END_REF]. In this paper, we describe so far never published features of the BLAD/DifferentialAlgebra version, which are related to the management of parameters.

This software was designed with the long term idea that it could be used within a computational process which would involve differential elimination as an intermediate step. Such a situation would possibly occur within a software dedicated to the investigation of the properties of models presented or internally encoded by systems of nonlinear differential equations.

In such a situation, it is indeed essential for the output of the software to be specified algebraically. The argument is easy to understand over the more familiar Gaussian elimination: consider a square linear system A x = b to be solved. If one uses the L U decomposition of A then one gets an algebraic relation A = L U between three matrices which permits to transform the problem "solving A x = b" into "solving L y = b (first) and U x = y (second)". If one uses the tacky Gaussian elimination, the former reduction problem gets difficult to explain, because of the lack of algebraic relation between the input matrix A and the elimination output.

The Different Types of Parameters

A polynomial differential system which encodes a mathematical model is likely to involve a large set of parameters (in models all constants are denoted by letters: masses are denoted m, lengthes . . . ) and, during the simplification process undertaken by the Rosenfeld-Groebner algorithm, the right way to manage these parameters may depend on the parameter. We may distinguish the following different types of situation:

A this type corresponds to parameters numerically known though denoted by a letter for legibility (the constant g of Newton F = m g law);

B this type corresponds to parameters which are not supposed to be perfectly defined (in biological modeling or in epidemiology, models do not follow from physical laws and are rough approximations of the real phenomenon). It thus hardly makes sense to assume that these parameters satisfy any specific relation;

C this type corresponds to parameters constrained by algebraic relations which are known in advance (the three edges of a rectangle triangle); D this type corresponds to parameters constrained by algebraic relations which are not known in advance.

One may then possibly want to discuss the result of the elimination process by considering the vanishing or non vanishing of various algebraic relations arising during computations.

Mathematical Handling of the Different Types

From a computational point of view, the management of parameters essentially amounts to control the way the Rosenfeld-Groebner algorithm splits cases. From an algebraic point of view, this case splitting mechanism is controlled by acting on the base field of the polynomial ring. This is an important point of the design of the DifferentialAlgebra package and a controversial one! Indeed, in algebraic text books, polynomials are always considered as elements of a single, well-defined polynomial ring. One is not supposed to change the ring to which a polynomial belongs. Switching a polynomial from one ring to another is actually possible but one then needs to define a ring homomorphism and consider the image of the polynomial through this homorphism. In constrast, in applied text books -or in computer algebra software such as MAPLE -equations arise first and the authors (who may be brilliant scientists) may not even bother to define the algebraic structure these equations belong to. The DifferentialAlgebra package approach is somewhere in between: the equations may be entered first; the polynomial ring they belong to is specified at simplification time as an argument to the RosenfeldGroebner function. The result of this design is 1) soundness: a polynomial or an ideal always belongs to a precise ring when a computation is undertaken; 2) flexibility: the ring need not be uniquely defined 1 .

Let us thus assume that we want to process a polynomial differential system with numerical coefficients in the field of the rational numbers and depending on a single differential indeterminate y and a single parameter θ. Assume that we have a single derivation δ acting on the differential indeterminate, yielding derivatives. Introduce a symbol x such that δ x = 1 so that the differential indeterminate can be viewed as an unknown function y(x) and the derivation as d/dx. The "single" assumption if for legibility only.

For types A and B, the differential polynomial ring of the computation is

K{y} where K = Fr(Q{θ, x}/{δθ, δx -1}) .
With words, the base field K is the field of fractions of the differential polynomial ring Q{θ, x} quotiented by the prime differential ideal {δθ, δx -1}. In the ring K{y}, every nonzero algebraic relation constraining θ is a nonzero element c of the field K. It is thus invertible and the equation c = 0 is equivalent to 1 = 0. The Rosenfeld-Groebner algorithm will thus discard any case leading to such a relation. For type C, let us denote t the ideal of Q[θ] of the known relations which constrain the parameter θ (in the case of a single parameter, t is a principal ideal but this needs not be true in the general case). Then, the differential polynomial ring of the computation is K{y} where K = Fr(Q{θ, x}/{t, δθ, δx -1}) .

Observe that, if t is any polynomial of t then δt is a linear combination of monomials which all admit δθ as a factor. Thus {t, δθ} ∩ Q[θ] is equal to the radical of t. Assume t is prime (this is the usual case). Then {t, δθ, δx -1} is so and K is a field (dropping this assumption, K would be isomorphic to a direct product of fields -one field per isolated associated prime ideal of t). In the ring K{y}, an algebraic relation constraining θ is either zero (if it belongs to {t, δθ}) or invertible (if it does not). The Rosenfeld-Groebner algorithm will thus discard any case leading to a relation which does not belong to {t, δθ}.

For type D, the differential polynomial ring of the computation is

R{y} where R = Fr(Q{x}/{δx -1}){θ}/{δθ} . 1 Let us stress that we do not think that the algebraic text books presentation implies that a polynomial belongs to a precised single ring. The authors (who may also be brilliant scientists) may just want to fix one of many different possible rings before investigating algebraic issues (the choice of the ring being out of the scope of the book). A somewhat similar situation arises in numerical computing: in numerical text books, algorithms are studied independently of any precise question and it is common to focus on relative errors exclusively; however, as soon as a precise question is investigated, absolute errors enter in consideration.

With words, the parameter θ is withdrawn from the base field. It is considered as a differential indeterminate (as y is) implicitly constrained by the relation δθ = 0. The Rosenfeld-Groebner algorithm will thus allow any case leading to any relation in Q[θ] \ Q.

Computational Handling of the Different Types

Thanks to the regular differential chains algorithmic theory, the Rosenfeld-Groebner algorithm does not rely on any Gröbner basis computation2 . See [START_REF] Boulier | An Equivalence Theorem For Regular Differential Chains[END_REF] for recent developments on this theory.

Regular differential chains are particular cases of triangular sets. The implementation idea is the following. Assume that the Rosenfeld-Groebner algorithm processes some differential polynomial system Σ in one of the rings introduced in section 3. This algorithm actually computes a tree of quadruples A, D, P, S where A is both the set of the (somehow) already processed equations and a regular chain. In the basic case (no parameters), the algorithm starts with the root A, D, P, S = ∅, ∅, Σ, ∅ and stops on leaves of the form A, D, P, S = A, ∅, ∅, ∅ where A is a regular differential chain.

In the case with parameters, is is assumed that the relations needed to define the fields K or the ring R introduced in section 3 are presented as a regular differential chain A (this actually is a restriction on the ideal t of the type C). Then the algorithm starts with the root A, D, P, S = A , ∅, Σ, ∅ and stops on leaves of the form A, D, P, S = A, ∅, ∅, ∅ where A is a regular differential chain containing A as a regular differential subchain (see [START_REF] Boulier | An Equivalence Theorem For Regular Differential Chains[END_REF]Proposition 38] for a related theorem). Here are a few remarks:

• all base field defining relations are thus handled by the general algorithms for regular differential chains: zero decision, regularity decision, normal forms methods;

• relations such as δθ = 0 or δx = 1 are detected and passed as arguments to any algorithm (such as the differentiation algorithm of differential polynomials) which would otherwise waste computation time;

• depending on the types of parameters, some cases may not be generated;

• unless required, the equations of A such as δθ = 0 are not displayed by a DifferentialAlgebra function such as Equations;

• in the partial differential case, this mechanism is extended to handle functions depending on a restricted set of independent variables. For instance, assume two derivations δ x and δ t encoding partial derivatives with respect to x and t are defined; then a function θ(x) is handled as a differential indeterminate θ subject to δ t θ = 0.

Academic Examples

In this paper, we only provide academic examples. A more interesting one for type C is described in [START_REF] Boulier | Chemical Reaction Systems, Computer Algebra and Systems Biology[END_REF] and shows how the Michaelis-Menten formula can be obtained by a differential elimination process. Relations among the parameters come from hypotheses on conservation laws and equations meant to rename constants.

Type D Example

Let us start by the simplest type D. The empty parentheses following θ in the definition of the differential polynomial ring indicate that this "differential indeterminate" does not depend on x i.e. is a constant. The differential indeterminate u can be viewed as an unknown function u(x). The symbol u[x] denotes the derivative δu: The possible vanishing of constant factors of p is indeed discussed (polynomials in regular differential chains should be viewed as left hand sides of equations):

> Equations (ideal); 2 [[theta], [theta -1], [u[x] -4 u], [u]]
In the two last regular differential chains, the equation δθ = 0 is not displayed. We can get it if desired:

> Equations (ideal, fullset=true); 2 [[theta], [theta -1], [u[x] -4 u, theta[x]], [u, theta[x]]]

Types A and B Example

Let us now illustrate types A and B. For this, one introduces K as a differential field extension of Q defined by one generator θ and no relation, apart the implicit one δθ = 0. Observe that the field keyword is a placeholder. It will only be interpreted when passed as an argument to the Rosenfeld-Groebner algorithm:

Type C Example

Let us finish with type C. For this, one transforms K into a differential field extension of Q defined by one generator θ and a prime differential ideal of relations. This prime differential ideal is presented by a regular differential chain A : The Rosenfeld-Groebner algorithm now discards any case leading to a constraint on θ different from the one specified in A :

>

  with (DifferentialAlgebra): > p := theta*(theta-1)*(u[x]**2 -4*u); 2 p := theta (theta -1) (u[x] -4 u) > R := DifferentialRing (derivations = [x], blocks = [u,theta()]); R := differential_ring The Rosenfeld-Groebner algorithm produces four regular differential chains: > ideal := RosenfeldGroebner ([p], R); ideal := [regular_differential_chain, regular_differential_chain, regular_differential_chain, regular_differential_chain]

>

  K := field (generators = [theta]); K := field(generators = [theta]) The Rosenfeld-Groebner algorithm computes a regular differential chain decomposition of the differential ideal {p} in a differential polynomial ring obtained by changing the base field of R. It actually discards the two cases leading to constraints on θ: > ideal := RosenfeldGroebner ([p], R, basefield=K); ideal := [regular_differential_chain, regular_differential_chain] > Equations (ideal); 2 [[u[x] -4 u], [u]]

>

  with (Tools): > Aprime := PretendRegularDifferentialChain ([theta-1],notation=jet,R); Aprime := regular_differential_chain > K := field (generators = [theta], relations = Aprime); K := field(generators = [theta], relations = regular_differential_chain)

>

  ideal := RosenfeldGroebner ([p], R, basefield=K); ideal := [regular_differential_chain, regular_differential_chain] > Equations (ideal); 2 [[u[x] -4 u, theta -1], [u, theta -1]]

This was not the case in the early versions of the algorithm -as the name of the algorithm indicates it.


François Bonneau

Arnaud Botella

Antoine Mazuyer

Pierre Anquez

Benjamin Chauvin

Guillaume Caumon

Francois Bonneau

Guil- Laume Caumon

RINGMesh: An open source platform for Shared Earth Modeling

come

Introduction

A geological model aims at describing, representing and characterizing subsurface heterogeneities. It gathers raw data and their analysis in the light of the knowledge hold by an expert who has built the model. There are as many structural interpretations as people looking at the data [START_REF] Bond | What do you think this is ? "Conceptual uncertainty" in geoscience interpretation[END_REF][START_REF] Bond | Uncertainty in structural interpretation: Lessons to be learnt[END_REF]. Shared Earth Modeling (SEM) aims at sharing data, interpretations and uncertainties collected during the modeling chain [START_REF] Ringrose | Reservoir Model Design[END_REF]. It has to provide to any modeler/user (geophysicists, geologists, reservoir engineers and business experts) the level of knowledge they need to support decisions. As each subdiscipline tends to use its own tools, most initiatives since the 90's have aimed at developing standards both conceptually (ontology-based [START_REF] Perrin | Geological consistency : an opportunity for safe surface assembly and quick model exploration. In: 3D Modeling of Natural Objects, A Challenge for the 2000[END_REF][START_REF] Mastella | Knowledge Management for Shared Earth Modelling[END_REF]) and in standard file formats (RESCUE and then RESQML1 [START_REF] Morandini | Shared Earth Model Data Exchanges between Commercial Modelling Applications and In-House Developments, Demonstrated on Actual Subsurface Data[END_REF][START_REF] Morandini | RESQML Version 2.0.1 Makes it Easier To Update a Reservoir Model[END_REF]). Indeed, standards are a key component to enable communication between the various communities that build and use subsurface models (see [START_REF] Perrin | Shared Earth Modeling: knowledge driven solutions for building and managing subsurface 3D geological models[END_REF] for a review). These initiatives significantly facilitate linear workflows whereby each step provides inputs for the next one (Fig. 1 (a)). We suggest that more is needed to further foster innovation in subsurface modeling.

Fundamentally, modeling should support coupled physical processes in complex geological settings at the scale of interest. Both in bottom-up and top-down approaches, the feedback coming from physical simulations (reservoir production data, subsidence, acoustic sensors, etc.) must be re-injected spatially to better characterize geological features and reduce uncertainties (e.g., [START_REF] Irakarama | Appraising structural models using seismic data: problem and challenges[END_REF]). This is a highly challenging problem, which calls for a common geological data model that is able to represent knowledge at multiple scales and transfer information gained from various sources. We suggest that implementing such a data model and associated functionality in a common, extensible and open-source software library can significantly help academic and industry researchers to embrace these challenges. When several applications are involved, it can ensure the validity and the retro-compatibility of models when migrating from one software to another, which is a mandatory for editing interpretation performed at any stages of the modeling chain.

In this work, we present RINGMesh2 , which is a software platform supporting the alternative integrated subsurface modeling workflows described in Fig. 1 (b-c)). More precisely, RINGMesh propose an open-source geomodel data structure that can be used to implement operations on subsurface models. It allows for connecting several applications by converting the geomodels into various file formats. Finally, it enables the direct use of third-party mesh data structures which clear the path for a tight coupling of physical and geological modeling.

An open-source data model to deal with the geometry and the topology of geological objects RINGMesh implements a representation of subsurface organization and geometry using a set of geological and geometrical entities (Fig. 2). The GeoModel data structure combines geometrical and topological information with geological knowledge and interpretations in an abstract and flexible design (see section ). Each geometrical entity is associated with a mesh data structure that hold the entity discretization. [START_REF] Ringrose | Reservoir Model Design[END_REF]. This involves transfer of relevant information between distinct software applications, for which model validity on both ends is difficult to ensure. (b) RINGMesh provides features to check and fix some of these validity issues through a common and open data model implementation. (c) RINGMesh's abstract data structure gives the opportunity develop new methods or wrap existing implementations to address integrative modeling challenges.

The default mesh data structure is implemented in the geometric algorithm library "Geogram"3 . For instance, regions can be meshed with tetrahedra, hexaedra, prisms and pyramids; surfaces can be meshed with triangles or any other polygons; lines can be meshed with segments; and corners are isolated mesh nodes.

The "GeoModel" geometrical description is similar to the Gmsh4 (Geuzaine and Remacle, 2009) and MSTK5 (Garimella, 2004) representations. It is made of geometrical entities of various dimensions. Each region is bounded by surfaces, each surface is bounded by lines and each line is ended by two corners. A geometrical entity stores a list of upward and downwards entity adjacencies. The topology of the GeoModel can be investigated through "boundary/incident" relationships of geometrical entities.

A geological entity is a set of geometrical entities that shares the same geological feature (e.g., a horizon, a fault or a layer). For instance, a horizon is a set of surfaces. A "parent/child" relationship drives the access from the parent geological entity to its constitutive geometrical entities (children). Additionally, the conceptual data and knowledge collected during the modeling chain (such as the stratigraphic col-umn) are associated to geological entities. The current implementation needs to be improved to integrate an exhaustive description of the ontological based geological model [START_REF] Perrin | Shared Earth Modeling: knowledge driven solutions for building and managing subsurface 3D geological models[END_REF].

Together with the data model, RINGMesh provides an extensible set of test functions which check the geometrical and geological consistency of a particular GeoModel instance. Figure 2 RINGMesh implements the "GeoModel" to hold the geometry and the topology of geological models. The GeoModel (a) is made of geological entities (b). The organization, the geometry and the connectivity of geological objects can be accessed through the "parent/child" relationship defined between geological entities and geometrical entities (c). The connectivity of geometrical entities is explicitly defined by "boundary/incident" relationships.

A smart converter of geological model file format

RINGMesh gathers an array of features to load, export, repair and visualize "GeoModels". These converters are essential, as the main current use of RINGMesh is to communicate with external geomodeler, mesher and simulation software (Fig. 1 (b)). Classically, exchanging data between applications can be time-consuming, and errors are often detected quite late in the process. RINGMesh helps the early detection of transfer problems by performing validity checks, and helps their automatic resolution by proposing classical fixes for geometrical issues.

Ongoing development efforts aim at further supporting classical databases and standard file formats, brings the opportunity not only to transfer but to convert the information between software. These features are essential in a SEM project.

An abstract mesh data structure to bind software and libraries Since the first implementation of the GeoModel [START_REF] Pellerin | RINGMesh: A programming library for developing mesh based geomodeling applications[END_REF], an abstraction level has been added to manipulate purely topological GeoModel. This gives the opportunity to consider various topological representations without relying on a particular mesh data structure. This is achieved by an abstract mesh object supporting the RINGMesh features. The default mesh data structure implementation is provided by Geogram, but alternative data structures can be used (e.g., a discretized mesh or a mathematical parameterization (spline, NURBS)). This enlarges the representation capability of RINGMesh and may provide a direct wrapper to physical simulators, meshers and geomodelers (Fig. 1 (c)). We are currently working on a proof of concept with MFEM6 by specifying the RINGMesh abstract mesh with the MFEM data structure. Finite Element simulations can then be directly performed on RINGMesh entities. This clears the path for direct physical feedback on the geological model.

Conclusion

RINGMesh is an open-source platform to support R & D projects for integrative numerical geology. This initiative aims at providing some of the missing pieces needed for a true Shared Earth Model implementation.

(1) It enhances the communication between meshers, solvers and geomodelers to reduce subsurface uncertainties.

(2) It adapts input files to GeoModel requirements by performing validity tests. This smart converter needs to be enhanced both to diversify input and output file formats and to integrate more validity check and repairing functionalities.

(3) The GeoModel features use a generic mesh data structure. Specifying it allows to encapsulate various technologies in the heart of the data model. The idea is to enable flow and mechanical simulation by wrapping software directly in the RINGMesh platform. We see this as an essential component to perform joint inversion on realistic geological models. (4) RINGMesh proposes an extensive and mutable design. We encourage the community to develop their own applications and workflows (Botella, 2016;[START_REF] Anquez | Automatic sealing and simplification of 3D geological surface models using topology recovery[END_REF][START_REF] Chauvin | Validating novel boundary conditions for 3D mechanics-based restoration: an extensional sandbox model example[END_REF][START_REF] Mazuyer | Stress estimation in reservoirs using an integrated inverse method[END_REF] and to feed the project with generic tools and features.

(5) It is developed in C++11 and continuously tested using integration tools to ensure stability. Non regression and unit tests are run before and after any contribution to the main repository on github (https://github.com/ringmesh/RINGMesh). ( 6) As the project evolves quickly, an up-to-date list of features and tutorials are available on the website of the project (http://ringmesh.org).

Figure 1

 1 Figure1(a) Shared Earth Modeling is often based on the classical modeling chain[START_REF] Ringrose | Reservoir Model Design[END_REF]. This involves transfer of relevant information between distinct software applications, for which model validity on both ends is difficult to ensure. (b) RINGMesh provides features to check and fix some of these validity issues through a common and open data model implementation. (c) RINGMesh's abstract data structure gives the opportunity develop new methods or wrap existing implementations to address integrative modeling challenges.

http://www.energistics.org

http://ringmesh.org

http://alice.loria.fr/software/geogram

http://gmsh.info

https://github.com/MeshToolkit/MSTK 80 th EAGE Conference & Exhibition 2018 11-14 June 2018, Copenhagen, Denmark

http://mfem.org 80 th EAGE Conference & Exhibition 2018 11-14 June 2018, Copenhagen, Denmark

Acknowledgements

This research was performed in the frame of the RING project at the Université de Lorraine. We would like to thank for their support the industrial and academic sponsors of the RING-GOCAD Consortium managed by ASGA. We particularly thank Total S.A. for their active participation in RINGMesh development. We also acknowledge Inria for the Geogram library used in RINGMesh.


Jean J Moreau
email: moreau@lmgc.univ-montp2.fr

Michel Jean
email: mjean@imtumn-imt.fr

NUMERICAL TREATMENT OF CONTACT AND FRICTION: THE CONTACT DYNAMICS METHOD

Difficulties occur in the numerical treatment of mechanical systems when the unilateral constraints of mutual impenetrability of members are taken into account. At every time-step, algorithms have to be ready to handle the possible break.ing of some contacts or, on the contrary, their sudden introduction by collisions. Furthermore, the modeling of dry friction, generally present in the case of contact, rests on mathematically irregular relationships. It is shown in this paper how the numerical approach named 'Contact Dynamics' allows one to calculate evolutions in dynamic or quasi static regimes and also to analyze the possible equilibria. This method faces the essential nonsmoothness of the concerned problems without ressorting to the regularizing approximations applied by other commonly used techniques. It rests on time-discre tization schemes of the implicit type, allowing for coarser time steps, at the price of applying nonsmooth iteration procedures.

INTRODUCTION

Various engineering situations cal! for the mechanical analysis of collections of bodies which, instead of the articulations considered in the traditional mechanics of systems, are only submitted to the unilateral constraints of not interpenetrating each other. At a priori unknown instants, such bodies may enter into contact or separate out, but they can never overlap.

For example, gra11ulatematerials, studied from the static or the dynamic standpoints, are systems of this sort. Problems raised by the behaviour of a work of maso11ry submitted to a11 earthquake let themselvcs be formulated in a similar way, as well as various othcrs in Gcomechanics [START_REF] Jean | Frictional contact in collections of rigid or de formable bodies: numerical simulation of geomaterials[END_REF]. In these instances, the concerned bodies are primarily treated as perfectly rigid, but the numerical in vestigation of the formi11g of mate rials (drawing, rolling, etc.) show the same type of constrainl imposed to portions of deformable me dia (see e.g. [START_REF] Jean | Numerical methods for three dimensional dy namical problems[END_REF], Raous et al., 1995). Also in Robotics, some authors are currently taking inlo accounl, for simulation or control, the unilatcral constraints of impenetrability [START_REF] Brogliato | On the transition phase in robotics. Part I: Impact models and dynamics. Part 11: Control[END_REF] (Pierrol el al., 1994).

In terms of the configuration parameters of the system (in numerical computation, there is always a finite number of them, even if the considered bodies are deformable), impenetrability constraints are expressed by inequality conditions. Furthermore, when contact occurs, solving statical or dynamical problems re quires some phenomenological information about the co11tact forces that the bodies exert on each other.

Even in the ideal case of frictionless contact, due to the w1i lateralcharacter of the impenetrability constraints, the appraisal of contact forces in Dynamics when several contacts are simul taneously effective has long been known to need careful discussion.

Additional difficulties corne fromjriction at the contact locus, which is rarely negligible. Tribology describes it as a complex phenomenon so that the data needed for a quite accurate prediction of its effects are, in most engineering situations, unavailable.

In practice the Coulomb law of dry frictio11 is very commonly invoked, for it retains the most salient features of frictional contact. But even this fairly simple law raises analytical and numerical difficulties, since the relationship it formulates between the contact forces and the relative velocities of the bodies at the respective contact points is mathematically irregular. In the event of zero relative velocity, the la"• only requires of the contact force to belong to a certain cone (the cone of static friction). For nonzero relative velocity, the law connects the tangential component of the contact force with its normal component and with the direction of the velocity vector. This direction may vary sharply when the Yelocity approaches zero, while the normal component itself is not usually given, but depends on the whole dynamical analysis. AL equilibrium, familiar examples show that situations involving dry friction are frequently hyperslatic. For evolution problems, one should no more expect the uniqueness of solution. Furthermore, in dynamical evolutions subject to impenetrability constraints, collisio11s implying velocity jumps are likely to occur.

This constilules another sort of analytical irregularity white, from the phenomenological standpoint, collisions may be affected by many aspects of the material behaviour which make Lheir outcome dif f icult to predict with greal precision. ln Operation Research, when problems of Oplimization or Control in the presence of inequality conditions are treated, mathc-matical difficulties of the same nature are met, though somewhat milder. The analytical and numerical tools developed to solve the latter problems constitute the mathematical domain called No11smooth Analysis. By analogy, questions of the sort considered here, as well as various others arising from plasticity, shape optimi zation, etc. have been said to belong to No11smooth Mechanics (Moreau and Panagiotopoulos, 1988;Moreau, Panagiotopoulos and Strang, 1988).

A natural approximation procedure, for the numerical treatment of such problems consists in replacing the impenetrability constraints by some short range repulsion laws which enter into ac tion when two bodies corne close to each other. This is a me chanical analog to the penalization technique, classical in Constrained Optimization. The dynamics of the approximant mechanical system is governed by proper differential equations to which standard numerical schemes, usually of the explicit type, are applied. Similarly, the nonsmooth character of dry friction may be numerically overcome by replacing Coulomb Iaw by some regularized approximant.

After the pioneering work of P. [START_REF] Cundall | A computer mode) for simulating progressive large scale movements of blocky rock systems[END_REF],_many authors have developed numerical methods of this sort and used them to produce significant results. Depending on the respective authors environment and on technical specificities, such methods are named 'Distinct Element method', 'Discrete Element method' or 'Molecular Dynamics methods'. Sorne commercial programs rest on this numerical strategy. The main drawback is that the sake of precision in approximating impenetrability constraints requires of the repulsion laws to be stiff, so that the numerical integration of the corresponding differential equations needs fine time-steps. Furthermore, it may be difficult to distinguish the possible artefacts generated by such artificial elasticity from the effects of physical deformability. As a recent paper where these technical aspects are discussed, one may consult Drake and Walton (1995).

In contrast, the Colllacl Dynamics method, presented in this communication, faces nonsmoothness without resorting to any re gularizing approximation. The key feature is that the impenetrability inequalities are treated at the level of velocities. Clearly, if in some state of the system a certain contact is in effect, the corresponding impenetrability inequality requires of the consequent velocities (i. e. the time-derivatives of the configuration parameters on the right of the concerned instant) to verify a certain linear inequality. Thus, for each possible configuration of the system, the set V of the 'right admissible' velocity vectors is a convex polyhedral cone in R 0 . lt turns out that the information to be entered, conceming contact forces is also connected with the same cone. For instance, in the academic case of frictionless contact, the condition imposed to the element of R 0 which represents the contact forces is that of belon ging to the 'polar cone' of V. If friction is taken into account, contact forces have to be related to velocity by more complicated relationships but fitting very well into such a Convex Analysis fra mework. In Sec.3 below we introduce the concept of a complete contact lmv, a formalism which secures the mechanical and numeri cal consistency of the model while leaving widely open the possi bility of introducing phenomenological data.

The problem to be solved at each step of the time-discretization is to determine velocities and contact forces jointly, in order for the contact laws and the balance of mornenturn to be satistïed. The nurnerical schemes so constructed naturally turn out to be of the implicit type. This implicit character makes the algorithms ready, at each time-step, to face collisions, should one of them be detected as occurring during the time-interval of the step, a point developed in Sec.5 below.

2

The paper ends with a few examples of applications developed by the authors.

ANALYTICAL $ETTING

Let the configurations• of the investigated system be parametrized, at least Iocally, through generalized coordinates, say q = (q 1 , q 2 , ... , qn). If the system consists only of perfectly rigid bodies, this reduction to finite freedom usually results from some (bilateral) ideal constraints, such as the possible action of internai or external frictionless linkages. Also when deformable continuous bodies are involved, the numerical treatment rests on finite freedom approximations (modal representation or use of finite elements) so that the formalism presented in what follows applies as well.

After constructing the parametrization, one takes into account the constraints of impenetrability; their geometric effect is assumed expressed by a fini te set of inequalities fa(q) :S O, aE{l, 2, ... , k}

where f1, f2, ... ,fk are given fonctions. Equality fa=Ü corresponds to the occurrence of a contact between two members of the system or between one of these members and some externat obstacle. For brevity, the latter will be assumed fixed in the reference frame in use. The case of obstacles with prescribed motion would cal! for the time t being also a variable in fa and bring only notational compli cations to the sequel (see [START_REF] Moreau | Sorne numerical methods in multibody dynamics : application to granular materials[END_REF].

For every imagined motion t--+q(t) and for t such that the deriva tive q(t)ER" exists, the kinetic energy has a quadratic expression in q, say '.Eic( q, q) = t Aïj( q) qi qi

(2)

where A denotes a symmetric positive definite nxn-matrix.

As far as smooth, i.e. twice differentiable, motions are concer ned, the system Dynamics is governed by Lagrange's equations, written below as an equality in R" A( q) q = F(t, q. q) + L rU a (3)

The expression F comprises the q-dependent terms of Lagrange's equations and the covariant components, relative to the parametrization (q), of some applied forces supposed given as fonctions of time, position and velocity. The element rU of R 0 is made of the covariant components of the a priori unknown contact forces experienced by the system in case the contact fa = O holds. Their definition rests on the standard construction of the covariant components of forces located in physical space, connected as fol lows with the system kinematics.

Suppose that inequality fa:SÜ expresses the mutual impenetrabi lity of some pair of members of the system, say '.B and '.B', so that equality fa=O corresponds to these two bodies touching each other at some point of space denoted by Ma. One assumes it an isolated contact point, but other contacts, corresponding to different values of a, may also be in effect between the same bodies at the same instant. For every imagined motion t--+q(t) bringing the system through the considered contact position for some value of t, the ve locities 'V a and 'Il; of the respective partiel es of '.Band '.B' passing at point Ma let themselves be expressed as affine fonctions of the va-lue u of the derivative q. The same is thus true for the relative ve locity 'll a. ='ll a. -'11 � of '13 with respect to '13' at this point, say ( 4 ) where G a. : R n _...., R3 denotes a linear mapping, depending on t and q. No attention is paid al this stage to the imagined motion preserving contact or not.

Let !J?... O. denote the contact force that body '13 experiences at point Ma from body '13'; then '13' experiences from '13 the force -!f?...<X. Classically, the covariant components of this pair of forces are expressed by (5) with G�: R3--+R 0 denoting the transpose of G a. (the convention of implicit summation will never be applied to Greek indices).

Similar formulas hold if inequality f a. :S O represents the confinement of a part '13 of the system by some externat boundary. Then ra is made of the covariant components of !J?... a. alone. In both cases, the following relationship is found [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] to hold between the gradient af taq in R" and the normal unit vecteur n o. at point M a. to the two co�tacting bodies, directed toward '13 3"1x<!:Ü such that (6)

In ail the sequel, we shall assume that the rnapping G a. is 01110 from R 0 to R3; equivalently, its transpose G� is 011e-to-011e from R3 i11to R n . Only some special positions of certain linkages may give rise to 'wedging' effects which break this assumption.

CONTACT LAWS

To E.qs. (1), ( 3), (4), ( 5) one necessarily has to adjoin some phenomenological information regarding the contact forces. For instance, the ideal situation of frictio11/ess co/llact is described by asserting that the contact force !J?...O. (vanishing if contact a is not effective, i.e. if f a. (q)<O ) is a vector belonging to the half-line generated in R3 by na..

As an elementary description of dry friction, the law of Coulomb will be invoked in the sequel. Generally, let us cal! a contact law a relationship of the form law a. (t, q, 'll a. , !J?...O.) = true (7) imposed to the elements defined in the preceding Section.

An adequately designed formulation of contact laws greatly helps to face the difficulties raised by the unilateral character of the impenetrability constraints. The event of collisions, i.e. the sudden occurence of new contacts, which is expected to generate velocity jumps, is left for Sec.5. At the present stage, where only smooth motions are considered, one still has to face the possibility of the breaki11g of some contacts. The traditional approach consists in tentatively calculating the motion under the assumption that ail contacts present at a considered instant remain effective in the sequel. If the evaluation of contact forces in the course of such a motion yields, at a further instant, an unfeasible direction for some of them, one concludes that some contacts should break at this instant, so the continued motion has to be calculated differently. But easy examples show that the contacts which break are not nccessarily those for which unfeasible contact forces were just 3 found. In the ideal case of frictionless contacts, [START_REF] Delassus | Mémoire sur la théorie des liaisons finies unilatérales[END_REF] has elucidated the determination of instant accelerations under such circumstances and the same question was more recently investigated in terms of Convex Analysis [START_REF] Moreau | Quadratic programming in mechanics: dynamics of one-sided constraints[END_REF].

Strictly speaking, the matrix G a. and the ,•ector na. defined in the preceding Section make sense only in the case of contact, i.e. for a configuration q satisfying f a. (q)=O. In computation, as well as in existential studies, the definition of these elements has to be extended, in a smooth arbitrary way, to configurations q laying in a neighbourhood of the hypersurface f a. =0 of R". In particular, one should be ready to face a certain amount of violation of the impenetrability inequalities (one manages to keep it smalt).

Put

Xa(q) = { {'l1ER3: n o. .'ll<!:0} if �(q)<!:O R3 if �(q)<O
which may be called the set of the right-admissible values for the relative velocity of the two concemed bodies at the contact point.

Let LLS agree to say that a contact /aw such as ( 7) is complete if it implies 'll a. E.Xrt.i11 al/ circumsta11ces and !l<a=O if 'UaEinterior 'Jiw.. Through E.q. (4), the fact that a contact Law is complete is translated into a property of the element u of R 0 .0ne source of interest of the concept lies in the following k.inematical result [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF].

If. at some initial i11sta11t, ail impe11etrability i11equalities fag)

are ful/filled and if ail co11tacts possibly occurrùzg i11 the sequel are govemed by complete contact laws, thefu11ctio11 t--+q(t) co11structed by i11tegratio11 of the velocity jwictio11 t--+u(t), jullfil/s the impenetrability i11eq1uzlities i11 ail the subsequellf motio11.

The Contact Dynarnics method rests on time-discretization. The essential stage is to compute a discrete approximant of the velocity fonction t--+u(t), to which the position fonction t--+q(t) is simply connected, step by step, through integration. Starting from the value obtained for u at the beginning of an interval of the discretization, one has to assess the value at the endpoint. Discretizing the dynamical equation (3) amounts to express the momentum balance of the system over the interval, without explicitely referring to accelerations. If the contact forces which appear in this balance of momentum verify complete contact laws, the abo,•e kinematical result makes that the impenetrability inequalities ( 1) are autornatically taken care of.

Furthermore, the use of complete contact laws secures the correct handling of the possibility of some contact to get loose. Such is in fact the case for contact a if computation yields 'UaEinterior 'liw_ (i.e. na. 'U>O) hence !T<a=O. Clearly also the complete contact Law assumption implies that !l<a=Ü when fa(q)<O, since 'liw_(q) equals the whole of R3 in that case.

BASIC CONTACT DYNAMICS ALGORITHM

Let [t 1 , tF], tF = t 1 +h, denote an interval of the time discretization. Starting ,,•ith q 1 , u 1 , the approximate values of q and u at time tr, the objective is to calculate qF, uF, the approximate Yalues at the end point lF of this interYal.

Identification of contacts

By using qM = q 1 + h u 1 /2 as test position, the set of the contacts to be treated as active in the considered step is estimated to be

Discretization of the equation of Dynamics

Equation (3) may be dÎscretized in the form where S 13ER3 denotes the impulsion at contact 13, i.e. the integral of 9(11 over the time interval. In short (8)

Contact laws

In contact laws such as (7), which are assumed satisfied all over the time interval, the time-dependent quantities 'U a , !}( a have to be replaced by some values estimated to be typical of the interval. We choose here to have (7) hold between the impulsion s a and the . final �oe .( 9) Such a choice makes the present time-discretization scheme belong to the implicit sort, at least in what concems the velocity function u.

As a justification of entering impulsions instead of forces in (7) one may observe that Coulomb law states a positive/y homogeneous relationship between its two vector arguments.

Final position qp=q M + lh up 2

Iterative procedure By combining Eqs. ( 7) and ( 9), one obtains 'v'aEJ: law a (q M ,G a uF,s a )=true (10)

Solving the system of conditions ( 8) to (10) constitutes the heaviest part of the computation. Here is a relaxation technique, amounting to treat a succession of single-contact problems.

Let an e � timated solution u e Fu, S e � t i , j3 running throu _ gh J, . be obtamed, w1th (8) satisfied (the startmg guess may cons1st m taking zero values for the sll or in using values found at the preceding step). One attempts to construct a corrected estimate, say uc;r.s!r, by al � ering onl � s a , i.e. s J!r = s titi for 13 '°. a. The new esttmate 1s astramed to sattsfy (8), i.e. smce ffie old esttmate sat1sfies the same, éO!T =U esti + A -'o* (S a -s a .)

F F a corr esll (11)
and to satisfy (10), hence, after applying Ga to both members of (11 ),

VaEJ: law (G u esti + H (S a -s a .) S a )=true ( 12) a a F a corr esll ' corr where can be proved to be a symmetric positive definite 3x3 matrix. Solving ( 12) with regard to the unknown S c &r may, in some usual cases, be reduced to finding the zeros of a piecewise affine mapping. The above computation will then be iterated, with a 4 ranging cyclically through J. The decision of stopping iterations may be taken on observing the magnitude of Sc&r-s e iti and this tums out to be equivalent to checking the precision at which each pair 'll a .S a satisfies the corresponding contact Iaw. Observe that, provided this precision check is made, the operator H a in ( 12) may be replaced by any other mapping with zero limit at the origin: this may be used in tricks for accelerating convergence.

Clearly, this algorithm tolerates a certain amount of violation of the impenetrability inequalities. By adjusting the step-length and the stopping criterium, one may keep these errors arbitrarily small and prevent their accumulation.

The iterated calculation is very simple, but needs to be repeated many times in case of numerous contacts. Since the balance of momentum ( 8) is only preserved from one iteration to the next through the conservation condition ( 11), one should think of the possible ac�umulation of arithmetic errors. For safety, one may refresh u e F t I from time to _ time, by returni � g to ( 8) while keeping the approximate values obtamed for sll. This proves useful only for motions involving thousands of contacts.

Technically, let us also observe that in many usual applications, the nxn matrix A is constant and diagonal. o a is a 3 xn matrix, but only the elements corresponding to the two bodies involved in contact a are nonzero. So the treatment of large collections of bodies does not require the handling of large matrices.

. The convergence of this algorithm and the existence of a solution to the problem it addresses has only been proved in special cases [START_REF] Monteiro Marques | Dijferential Inclusions in Nonsmooth Mechanical Problems: Shocks and Dry Friction[END_REF]. Uniqueness cannot be expected to hold in general, since the mechanical problem of determining the reactions in a closely packed collection of rigid bodies (for instance a wall made of rectangular blocks) is usually hyperstatic.

The algorithm still works in analysing equi � ibrizun � ituations.

The balance of applied forces and contact forces stmply y1elds up = u 1 = 0 in such a case.

COLLISIONS

The sudden occurrence of new contacts produces velocity changes in some parts. The Jess deformable the colliding bodies are, the more brutal the phenomenon should be, involving large values of the contact forces at the impact locus. If the deformability of the bodies is full y taken into account, the problem is however not fundamentally different from others in the Dynarnics of Continua.

But if it has been decided to treat bodies as perfectly rigid, one has to face strict discontinuities in velocities, so that the smooth dynamic framework of Sec.2 cannot apply. In the majority of the papers devoted to the topic of 'Rigid Body Collisions' [START_REF] Brach | Mec/zanicallmpactDynamics(Rigid Body Collisions)[END_REF], it is attempted to formulate some collisio11 equations, connecting the values of velocities after the collision to the values they had before it. Traditionally, the intense effects which take place during a collision are assumed Iocalized in the vicinity of the impact locus. Then, a multiple scaling analysis of what happens m an 'infinitely small' domain, during the 'infinitely short' collisional episode allows one to take into account in more or Jess detail the material behaviour of the involved bodies.

Situations in which such a treatment is justified certainly exist, but in general the effect of a collision should not be local. For instance, material dissipation in the vicinity of the impact is not the only cause of the energy loss detected at the macroscopic level of obser v ation. Even if the concemed bodies are assumed perfectly elastic, energy conser v ation cannot be expected. Disturbances are t . o propagate from the collision locus to the whole system and also. ,r the latter is linked \\'ith some extemal support, to the outside world.

After contact recedes, vibrations are likely to persist somewhere. At the macroscopic observation level, this does not contradict the rigidity assertion, but the energy involved in microscopie agitation may not be negligible.

Also as a consequence of global deformation, a collision may, at the microscopie time-scale, split into several separate contact episodes: an example of such a double bouiice is calculated in closed form in (Timoshenko, 1948, Chap. 12). Finite element computation of the collision of two elastic bodies performed in our laboratory has shown the same. So the conception of a collision as consisting of a compression phase followed by a so-called restitution phase cannot be considered as general.

Ail this makes the outcome of a collision depend on many factors, in particular on the shape of the concerned bodies.

Still more severe difficulties arise from collisions being fre quently multiple, i.e .. several contact loci are involved at the same time. Such is the case if a colliding body is part of a cluster of objects already in contact. The propagation of disturbances into a cluster is a problem similar to that of sound in granular media One thus has to accept that any given mode! of collision can only have a limited scope. Every occasion of compa• ring its results with calibrated experiments should be seized, in order to estimate this scope as precisely as possible.

The algorithm described in Sec.4 is found to work consistently in the face of collisions, i.e. when the test position qM reveals some contacts which were not in effect at the preceding step. There only happens that the impulsive term in Eq.( 8) is no more the same order of magnitude as h. The elements u 1 and uF of Rn in this case are naturally interpreted as representing u-and u + , the values of the system velocity before and after the collision. Now, the decision of constructing a time-discretization algorithm of the implicit type bas led us to connect by Eq.( 10) each contact impulsion S a with the fi-1,al value uF. Contact laws are assumed complete in the sense preci sed in Sec.3. Clearly then S a can be nonzero only if na. 'll a p=O. Since 'll a F represents the after-collision value of the relative velocity at contact a, this means that all the contacts which take an active part in the collision process satisfy the traditional condition of zero restitu1io11. Using standard vocabulary, one may say that the algo rithm, as we just described it, treat ail contacts as i11elastic.

There remains to explain how situations of nonzero restitution can be handled in this framework. As a pragmatic way of overcoming the difficulties mentioned in the foregoing, one may decide to admit that a contact law of the form (7) holds between the contact impulsion s a and some formai velocity 'llaa constructed as a weighted average of the (known) vector 'llar and the ( unknown) vector 'll a F• The averaging may be effected separately for normal and tangential components, using different weights, say The numbers p and i; , with values in the interval (0, l], may be interpreted as tf�e nornZ'al and the tangential restitution coefficie11ts.

It turns out that, in the special case of the collision of two spherical objects, the velocity jumps deduced from this averaging trick coïncides with what has been proposed by several authors on the basis of a microscopie analysis and found in fairly good agreement with experimental data [START_REF] Foerster | Measurements of the collision properties of small spheres[END_REF]. The case p a =c a = 1, 5 with friction made equal to zero corresponds to energy preserving collisio11s [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF].

In the case of a two-body collision, p coïncides with the traditional restitution coefficient of Newto:r. But it should be stressed that the above formalism still yields plausible results in multiple collisions for which the Newton rule fails. For instance, what precedes applies to the problem, familiar in earthquake engineering literature, of the rocki11g of a slender rectangular block supported by a fixed horizontal plane. For simplicity, assume the lower edge slightly concave, so that contact can only occur through corners. Let the left corner remain in contact during an episode where the block rotates to the right, until the right corner collides. If at this time Newton's rule was applied to both contact points, this would yield zero normal velocity for the left corner, so no rocking could be found. On the contrary, our assumption of complete contact laws leaves the possibility for this normal velocity to be nonzero while the normal contact impulsion vanishes. The result of the calculation depends on the aspect ratio. of the block, in conformance with common observation.

For large collections of bodies such as samples of granular materials, one may have to handle many collisions in a single interval of the time-discretization even though they theorically are not simultaneous. This is the key of the efficiency of the method in the simulation of granular materials. Comparison with physical experiments on convection currents and size-segregation in shaken granulates have produced quite satisfactory results [START_REF] Moreau | Numerical experiments in granular dynamics: vibration-induced size segregation[END_REF]. By disclosing the values of experimentally inaccessible variables, the numerical simulation has perrnitted to understand the undelying mechanisms.

EXAMPLES OF APPLICATION

Ancient column submitted to ground shake

Force transmission in a sand pile

Sorne experimentalists have found that the pressure exerted upon the ground by a sand pile may present a local minimum near the center of the pile [START_REF] Smid | Pressure distribution under heaped bulk solids[END_REF]. The explanation of this effect is currently a malter of discussion among theoreticians. Figure 2 was produced by a numerical simulation intended, as part of these investigations, to explore the transmission of forces inside a pile. For clarity, the mode! is two-dimensional (corresponding to experiments on a pile of 'Schneebeli material', i.e. a collection of cylindrical objects). The pile construction has been computed as the successive addition of 3000 circular grains of dispersed sizes, falling down at a velocity of 25 cm/s, with slight lateral random dispersion. After that, the pile has been left to relax under gravity (actually, the observation described below is practically the same if made during the pile construction).

Two vertical rectilinear cuts are shown. For each of them, one draws the resultant of the contact forces exerted by the grains with centers on one side upon those with centers on the other side. This demonstrates that the section comprised between the two cuts has part of its weight counter-balanced by the action of the rest of the pile. However, with the present data, the effect is not large enough to produce a visible local minimum of ground pressure.

Two-dimensional stress-strain experiment

A two-dimensional mode! of granulate (or Schneebeli material), composed of 1024 randomly located disks as shown in Fig. 3, is being deformed between four walls. Horizontal walls are vertically loaded with opposite constant forces (average pressure o2=100 K.Pa). Vertical walls are moved inward with constant horizontal velocity v=5 cm/s. The granulate reacts on them by a pressure ol . Gravity is neglected in the computation.

The evolution investigated here makes the deformation E of the sample, with regard to its initial configuration, go from 2.625% to 5.250% .

On Fig. 3 are shown the displacement fields of the disk centers during two subintervals. The field on the left corresponds to some initial episode, E going from 2.625% to 3.280%, ahd the field on the right to a final episode, E going from 4.595% to 5:250%. In the initial episode, a band of quasi-stagnation in the 45 ° directio, n is visible, while in the final episode a similar feature appears with -45 ° direction. So a significant change in the de formation pattern occurs in the course of the investigated experiment. Such crises affecting the slow deformation of compact granular materials have been reported by experimentalists [START_REF] Meftah | Evidence of local 'seisms' of microscopie and 6 macroscopic stress fluctuations during the deformation of packings of grains[END_REF] and make laboratory tests harder to interpret.

The subsequent figures show how the crisis, which takes place near the middle of the evolution, reflects in various parameters.

Figure 1

 1 Figure 1 is part of a feasibility test of the Contact Dynarnics method in the simulation of the behaviour of a11cie11t momu11e111s under earthquakes. The column is made of 10 cylindrical blocks superimposed without mortar an of a cylindrical capitel. Instead of a proper earthquake, the ground is only affected by an oblique elliptic oscillation. (Development supported by the Commission of the European Communities, Environment programme, Contract N ° EV5V CT93 0300)

  Figure 4: 01/02 versus defonnation