Emilien Azéma

Farhang Radjaï

Robert Peyroux

Gilles Saussine

Force transmission in a packing of pentagonal particles

Keywords:

We perform a detailed analysis of the contact force network in a dense confined packing of pentagonal particles simulated by means of the contact dynamics method. The effect of particle shape is evidenced by comparing the data from pentagon packing and from a packing with identical characteristics except for the circular shape of the particles. A counterintuitive finding of this work is that, under steady shearing, the pentagon packing develops a lower structural anisotropy than the disk packing. We show that this weakness is compensated by a higher force anisotropy, leading to enhanced shear strength of the pentagon packing. We revisit "strong" and "weak" force networks in the pentagon packing, but our simulation data provide also evidence for a large class of "very weak" forces carried mainly by vertex-to-edge contacts. The strong force chains are mostly composed of edge-to-edge contacts with a marked zig-zag aspect and a decreasing exponential probability distribution as in a disk packing.

I. INTRODUCTION

Among singular features of granular media, force transmission has received particular interest during the last decade. The contact forces in model granular media, as observed by experiments and numerical simulations, are highly inhomogeneous and their probability density functions (pdf's) are wide [1,2,3,4,5,6,7,8,9]. The granular texture is generically anisotropic in two respects: 1) The contact normal directions are not random; 2) The force average as a function of contact normal direction is not uniform. The corresponding fabric and force anisotropies in shear are responsible for mechanical strength at the scale of the packing [10,11,12,13]. Another interesting aspect, first analyzed in Ref. [10] is the fact that the forces organize themselves in two distinct classes which contribute differently to fabric anisotropy, shear stress, and dissipation. In particular, the shear stress is fully transmitted via a "strong" contact network, materialized by force "chains". The stability is ensured by the antagonist role of "weak" contacts which prop strong force chains [10,14].

The force transmission properties have been for the most part investigated in the case of granular media composed of isometric (circular or spheric) particles. However, in various fields of science and engineering, the grains are seldom so "perfect". For example, elongated and platy shapes are encountered in biomaterials or pharmaceutical applications. Such shapes have unequal dimensions and induce thus a degree of anisotropy in the bulk behavior in addition to fabric and force anisotropies [15,16,17,18]. On the other hand, granular geo- * Electronic address: azema@lmgc.univ-montp2.fr materials are often composed of angular particles with plane faces as polyhedra. While rounded particles enhance flowability, angular shape is susceptible to enhance the shear strength, a factor of vital importance to civilengineering applications [18,19,[START_REF] Pena | Granular Matter In Press[END_REF]. The railway ballast is a well-known case where particle shape must be optimized to avoid excessive differential settlement under vertical loading [START_REF] Saussine | [END_REF]22,23]. In such circumstances, the analysis of force transmission is a key to improve performance.

In dealing with effects of particle shape, the issue is that a general quantitative description of particle morphology requires various shape parameters. For regular polygons in 2D, for instance, the only shape parameter is the number of sides (besides the diameter) whereas for irregular polygons more information is needed about the positions of the vertices in a reference system attached to the particle. In soil mechanics, angularity and roundedness are among basic parameters used to describe particle shapes [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. As far as force transmission is concerned, at least two parameters seem to be most relevant: 1) shape anisotropy (anisometry), which contributes to the anisotropy of stress transmission [15]; 2) facettedness, which allows for extended (face to face, edge to face and edge to edge) contacts between particles leading possibly to the formation of columnar structures within a granular assembly.

In this paper, we consider one of the simplest possible shapes, namely regular pentagons. Among regular polygons, the pentagon has the lowest number of sides, corresponding to the least roundedness in this category, without the pathological space-filling properties of triangles and squares. We seek to isolate the effect of edgeto-edge contacts on force transmission by comparing the data with a packing of circular particles that, apart from the particle shape, is identical in all respects (prepara-tion, friction coefficients, particle size distribution) to the pentagon packing. Both packings are subjected to biaxial compression simulated by means of the contact dynamics method. The presence of edge-to-edge contacts affects both quantitatively and qualitatively the microstructure and the overall behavior during shear. These contacts do not transmit torques, but they are able to accommodate force lines that are usually unsustainable in packings of disks.

This paper is organized as follows. We first present in Section II the numerical procedures and a brief technical introduction to the detection and treatment of edge-toedge contacts in the framework of the contact dynamics method. In Section III, we compare stress-strain and volume-change characteristics. Then, In Sections IV and V, we analyze the texture and force transmission features. In Section VI, we focus on the pentagon packing and we analyze the structure of force networks with vertex-to-edge and edge-to-edge contacts. The main results are summarized and discussed in Section VII.

II. NUMERICAL PROCEDURES

The simulations were carried out by means of the contact dynamics (CD) method [START_REF] Jean | Proceedings of Contact Mechanics International Symposium (Presses Polytechniques et Universitaires Romandes[END_REF][START_REF] Moreau | Novel approaches in civil engineering[END_REF]. The CD method is based on implicit time integration of the equations of motion and a nonsmooth formulation of mutual exclusion and dry friction between particles. This method requires no elastic repulsive potential and no smoothing of the Coulomb friction law for the determination of forces. For this reason, the simulations can be performed with large time steps compared to molecular dynamics simulations. We used the platform LMGC90 which is a multipurpose software developed in Montpellier, capable of modeling a collection of deformable or undeformable particles of various shapes [START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF].

A. Contact dynamics for polygons

The particles are rigid polygons exerting normal and shear forces, f n and f t , respectively, on each other. We attribute a positive sign to compressive normal forces. The relative normal velocity u n between two particles in contact is counted positive when they move away from each other. Then, the condition of geometrical contact between two particles is expressed by the following mutually exclusive alternatives:

f n 0 u n = 0 f n = 0 u n > 0.
(1)

In the same way, the Coulomb friction law involves three mutually exclusive conditions:

f t = -µf n u t > 0 -µf n f t µf n u t = 0 f t = µf n u t < 0 (2)
where u t is the sliding velocity at the contact and µ is the friction coefficient. The unknown variables are particle velocities and contact forces. These are calculated at each time step by taking into account the conservation of momenta, the constraints expressed by ( 1) and ( 2), and the dissipation of kinetic energy during inelastic collisions between particles (ref). We use an iterative research algorithm based on a nonlinear Gauss-Seidel scheme. The uniqueness is not guaranteed for perfectly rigid particles in absolute terms. However, by initializing each step of calculation with the forces calculated in the preceding step, the set of admissible solutions shrinks to fluctuations which are basically below the numerical solution.

Let us note that in molecular dynamics simulations, this "force history" is encoded by construction in the particle positions.

The research algorithm is applied to a set of potential contacts, identified or updated in each step. The contact detection between two bodies consists in looking for the overlaps of the portions of space they occupy. The treatment of the mechanical interaction requires additionally the identification of a common tangent plane (a line in 2D). Of course, contact may take place through a larger contact zone than a single point. Several algorithms exist for overlap determination between convex polygons [START_REF] Saussine | [END_REF][START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF]. In 2D simulations of the present paper, the detection of contact between two convex polygonal bodies was implemented through the so-called "shadow overlap method" devised by Moreau [START_REF] Saussine | [END_REF][START_REF] Dubois | Actes du sixième colloque national en calcul des structures -CSMA-AFM-LMS[END_REF], with reliability and robustness tested in several years of previous applications to various states of granular materials [18,[START_REF] Cholet | World Congress on Railway Research[END_REF][START_REF] Azéma | [END_REF].

In detection of contacts between two polygons, two situations arise: 1) If a single corner is found crossing an edge of the partner polygon, the direction of this edge is viewed as the tangent direction. By orthogonally projecting the intruding vertex onto the edge, one determines the penetration depth, while the nominal contact point is chosen at the center of this distance. Below, we will refer to this vertex-to-edge contact as "simple" contact. 2) In case of double intrusion, the common tangent line is fixed from as a mean between the two overlapping edges and a segment of this line is identified as the contact segment. The impenetrability between two particles at such an edge-to-edge contact is ensured by applying the contact laws (1) and (2) to only two points of the contact segment (Fig. 1). For this reason, we refer below to edge-to-edge contacts as "double" contacts. In practice, two forces are calculated at each double contact, but only their resultant and application point are material. In this respect, the choice of the two points representing a double contact does not affect the dynamics of the system.

B. Numerical samples

We generated two numerical samples. The first sample, denoted S1, is composed of 14400 regular pentagons of three different diameters: 50% of diameter 2.5 cm, 34% of diameter 3.75 cm and 16% of diameter 5 cm. The second sample, denoted S2, is composed of 10000 discs with the same polydispersity. Both samples were prepared according to the same protocol. A dense packing was first constructed following simple geometrical rules [30] and then compressed isotropically under a constant stress σ 0 = 10 4 Pa applied onto the right and top walls. The gravity was set to zero in order to avoid force gradients in the samples. The coefficient of friction was set to 0.4 between grains and to 0 with the walls. At equilibrium, both numerical samples were in isotropic stress state. The solid fraction was φ 0 = 0.80 for S1 and φ 0 = 0.82 for S2. The aspect ratio was h/l ≈ 2, where h and l are the height and width of the sample, respectively. Figure 2 displays snapshots of the two packings at the end of isotropic compaction.

The isotropic samples were subjected to vertical compression by downward displacement of the top wall at a constant velocity of 1 cm/s for a constant confining stress σ 0 acting on the lateral walls. The simulations were run up to a total cumulative vertical strain of 0.2 with a time step of 5.10 -4 s. The CPU time was 7.10 -4 s and 5.10 -4 s per particle and per time step on a G5 Apple computer. 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 [31] defined by

I = ε m p , (3)
where ε = ẏ/y is the strain rate, m is the total mass, and p is the average pressure. The quasistatic limit is characterized by the condition I ≪ 1. In our biaxial simulations, I was below 10 -3 .

III. STRENGTH AND DILATANCY

In this section, we compare the stress-strain and volume-change behavior between the packings of poly- gons (sample S1) and disks (sample S2). For the calculation of the stress tensor, we consider the "tensorial moment" M i of each particle i defined by [14,[START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF]:

M i αβ = c∈i f c α r c β , (4)
where f c α is the α component of the force exerted on particle i at the contact c, r c β is the β component of the position vector of the same contact c, and the summation is runs over all contacts c of neighboring particles with the particle i (noted briefly by c ∈ i). It can be shown that the tensorial moment of a collection of rigid particles is the sum of the tensorial moments of individual particles. The stress tensor σ for a packing of volume V is simply given by [14,[START_REF] Moreau | Friction, Arching, Contact Dynamics[END_REF]:

σ = 1 V i∈V M i = 1 V c∈V f c α ℓ c β , (5)
where ℓ c is the intercenter vector joining the centers of the two touching particles at the contact c. Remark that the first summation runs over all particles whereas the second summation involves all contacts in the volume

0 0.1 0.2 0.3 0.4 ε q 0 0.1 0.2 0.3 0.4 q/p S1 S2
FIG. 3: Normalized shear stress q/p as a function of cumulative shear strain εq for the samples S1 and S2.

V , with each contact appearing only once. We extract the mean stress p = (σ 1 + σ 2 )/2, and the stress deviator q = (σ 1 -σ 2 )/2, where σ 1 and σ 2 are the principal stresses. The major principal direction during vertical compression is vertical. The strain parameters are the cumulative vertical, horizontal and volumetric strains ε 1 , ε 2 and ε p , respectively. By definition, we have

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

ε p = dV V = ln 1 + ∆V V 0 , (7)
where V 0 is the initial volume and ∆V = V -V 0 is the cumulative volume change.

Figure 3 shows the normalized shear stress q/p for the samples S1 and S2 as a function of shear strain ε q ≡ ε 1ε 2 . For S2, we observe a classical behavior characterized by a hardening behavior followed by (slight) softening and a stress plateau corresponding to the residual state of soil mechanics [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF]. For S1, we observe no marked stress peak. The residual stress is higher for polygons (≃ 0.35) than for disks (≃ 0.28). This means that the polygon packing has a higher angle of internal friction ϕ defined by

sin ϕ = q p . (8)
Figure 4 displays the cumulative volumetric strain ε p for polygons and disks as a function of ε q . Both samples dilate and tend to isochoric deformation at large strains. It is remarkable that the polygon packing S1 initially dilates less than the disk packing S2. This behavior is reversed at larger strains with a crossover occurring after the peak state. Notice that the solid fraction is initially lower in S1 (0.80) than in S2 (0.82). This is because it is more difficult to obtain a compact packing with polygonal shapes by isotropic compression as a result of enhanced steric effects compared to disks. In other words, angular particles can form larger pores compared to rounded particles. The volumetric deformation can also be expressed in terms of the so-called "dilatancy angle" ψ defined by [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF] sin

ψ = ε p ε q . ( 9
)
The cumulative angle of dilatancy, i.e. during shear up to the residual state, is only slightly higher for the polygon packing than the disk packing. The plot of ψ as a function of ϕ, i.e. the so-called stress-dilatancy diagram, is shown in Fig. 5 for polygons and disks [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF]. Remarkably, both plots are parallel to the line ϕ = ψ with an offset ϕ 0 :

ϕ ≃ ϕ 0 + ψ. ( 10
)
The offset ϕ 0 is the friction angle at zero dilatancy. We have ϕ 0 ≃ 0.12 for disks and ϕ 0 ≃ 0.3 for polygons. This observation is in agreement with the arguments of Taylor [START_REF] Wood | Soil behaviour and critical state soil mechanics[END_REF][START_REF] Radjai | The Physics of Granular Media[END_REF] based on energy balance and recently revisited also in the case of cohesive granular media [START_REF] Taboada | [END_REF]. The higher level of ϕ for the polygon packing reflects the organization of the microstructure and the features of force transmission for each particle shape. This point is considered in more detail in the following section.

IV. GRANULAR TEXTURE

The granular texture, i.e. the organization of the particles and their contacts in space, is basically controlled by steric exclusions between the particles and force balance conditions [36]. The texture can be described in terms of various statistical descriptors pertaining to the forcebearing network of particles. At the lowest order, the compactness of the structure can be described in terms of both the solid fraction ρ and the coordination number z. The connectivity of the network can further be characterized by the fraction P (c) of particles having exactly c contact neighbors. These are scalar parameters or functions. At higher orders, the anisotropy of the texture is described by different "fabric tensors". We consider here these geometrical descriptors in order to identify the signature of particle shape.

A. Connectivity

The connectivity of the particles by force-bearing contacts is described at the lowest order by the average number z of contact neighbors per particle. The particles with no force-bearing contact are thus removed from the statistics. Note also that each double (edge-to-edge) contact for the polygons is counted once although double contacts are treated as two point contacts belonging to the contact segment (see section II). Fig. 6a displays the evolution of z for the pentagon packing (S1) and the disk packing (S2) as a function of ε q . The coordination number evolves to a steady-state value in both samples that is higher for S2 (≃ 3.85) than for S1 (≃ 3.75). The difference is, however, much less important than in the initial configuration (≃ 3.95 for S2 compared to ≃ 3.20 for S1) prepared by means of isotropic compaction.

It is also interesting to compare the two samples in terms of "contact lifetimes". Let us consider a reference configuration, e.g. the initial state of each sample. We follow the history of each contact listed in this state. In particular, we define γ as the fraction of persistent contacts of the initial list. During deformation, γ declines from 1 to 0 as an increasing number of initial contacts are lost due to particle rearrangements. Fig. 6b shows gamma as a function of ε q for S1 and S2. We see that, following a rapid initial falloff, γ decreases slowly in both samples but the rate of contact loss is globally higher for polygons than disks. We remark that even at ε q = 0.4, the contact list is renewed by only 50%.

The connectivity P (c) of the particles is plotted in Fig. 7 for S1 and S2 at ε q = 0.3. Interestingly, the two plots are nearly identical with a peak for c = 4. In both samples, the fraction of particles with 5 contacts is larger than that with 3 contacts. This shows that the connectivity does not reflect the difference in texture between the two packings although a qualitative difference exists as we shall see below by considering fabric anisotropy and force transmission.

B. Fabric anisotropy

The shear strength of dry granular materials is generally attributed to the buildup of an anisotropic structure during shear due to friction between the particles and as a result of steric effects depending on particle shape [37,[START_REF] Cambou | Powders and Grains 93[END_REF][START_REF] Radjai | Granular Materials: Fundamentals and Applications[END_REF]. Several methods have been used to quantify the fabric (structural) anisotropy of granular materials [13,[START_REF] Satake | Proceedings of the IUTAM symposium on deformation and failure of granular materials[END_REF][START_REF]Mechanics of Granular Materials[END_REF]. A common approach is to consider the probability distribution P (n) of the contact normals n which are generically nonuniform. In two dimensions, the unit vector n is described by a single angle θ, the orientation

S1 S2

FIG. 8: Polar representation of the probability density function P θ of the contact normal directions θ for the samples S1 and S2 in the residual state.

of the contact normal. The probability density function P θ (θ) of contact normals provides a detailed statistical information about the fabric. It is π-periodic in the absence of an intrinsic polarity for n.

Most structural information is generally condensed in the second moment of P θ , called fabric tensor [START_REF] Satake | Proceedings of the IUTAM symposium on deformation and failure of granular materials[END_REF]:

F αβ = 1 π π 0 n α (θ)n β (θ)P θ (θ)dθ ≡ 1 N c c∈V n c α n c β , (11)
where α and β design the components in a reference frame and N c is the total number of contacts in the control volume V . By definition, tr(F ) = 1. The anisotropy of the contact network is given the difference between the principal values F 1 and F 2 . We define the fabric anisotropy a by

a = 2(F 1 -F 2 ). ( 12
)
For fix coordinates, with the x-axis pointing along θ ′ , we define also a "signed anisotropy" a ′ by

a ′ = 2(F 1 -F 2 ) cos 2(θ F -θ ′ ), (13)
where θ F is the major principal direction of the fabric tensor. For θ ′ = θ F , we have a ′ = a. The signed anisotropy corresponds to the second term of the Fourier expansion of P θ (θ) and it is useful whenever the direction of anisotropy is not constant. Figure 8 displays a polar representation of P θ (θ) for the samples S1 and S2 at ε q = 0.3. We observe a nearly isotropic distribution for the pentagon packing in spite of shearing whereas the disk packing is markedly anisotropic. This is a surprising observation in view of the higher shear strength of the pentagon packing (Fig. 3). It is also counterintuitive as one expects that double contacts should allow a polygon packing to build more easily an anisotropic structure.

The evolution of a ′ is shown in Fig. 9 as a function of ε q for S1 and S2. The privileged direction of the contacts, corresponding to θ F , is vertical in both packings. In both cases, a ′ increases from 0 (as a result of the initial isotropic compression) to its largest value in the residual state. The anisotropy stays quite weak in the pentagon packing whereas the disk packing is marked by a much larger anisotropy, increasing to ≃ 0.3 and then relaxing to a slightly lower value in the residual state. As we shall see below, the low anisotropy of the pentagon packing results from a particular organization of the force network in correlation with the orientations of simple and double contacts in the packing (section VI). We will also show that the large shear strength of the pentagon packing is a consequence of a strong force anisotropy in this packing (see next section).

V. FORCE TRANSMISSION

In this section, we analyze the anisotropy and inhomogeneity of force networks in the packings of pentagons and disks. This leads us to consider the contributions of force and texture anisotropies to average shear stresses.

A. Force anisotropy

The angular distribution of contact forces in a granular packing can be represented by the average force f (n) as a function of the contact normal direction n. We distinguish the average normal force f n from the average tangential force f t formally defined by [13]

     f n (θ) = 1 Nc(θ) c∈S(θ) f c n , f t (θ) = 1 Nc(θ) c∈S(θ) f c t , (14)
where f c n and f c t are the normal and tangential forces, respectively, acting at the contact c (according to a sign convention attributing positive values to the normal forces), S(θ) is the set of contacts with direction θ ∈ [θ -∆θ/2, θ + ∆θ/2] for angle increments ∆θ, and N c (θ) is the number of contacts in S(θ). By definition, the two functions f n and f t are πperiodic. After sufficiently long monotonous shearing, these functions can be approximated by their Fourier expansions truncated beyond the second term [13,[START_REF] Radjai | Granular Materials: Fundamentals and Applications[END_REF]:

f n (θ) = f {1 + a n cos 2(θ -θ n )} f t (θ) = f a t sin 2(θ -θ t ), ( 15
)
where f is the average force, a n and a t represent the anisotropies of the normal and tangential forces, respectively, and θ n and θ t are their privileged directions. In Fig. 10, the functions f n (θ) and f t (θ) are displayed in polar coordinates at ε q = 0.3. The pentagon and disk packings show pronounced force anisotropy with a stronger anisotropy in the case of pentagons both for normal and tangential forces. These plots can be fitted by harmonic functions [Eq. (15)] in order to estimate the force anisotropies a n and a t . However, it is more convenient to estimate the anisotropies through the following "force tensors":

       H (n) αβ = π 0 f n (θ)n α n β dθ, H (t) αβ = π 0 f t (θ)n α n β dθ. ( 16
)
It is easy to see that tr(H (n) ) = tr(H (t) ) = f , and by identification with (15) we have

     a n = 2 H (n) 1 -H (n) 2 H (n) 1 +H (n) 2 , a t = 2 H (t) 1 -H (t) 2 H (t) 1 +H (t) 2 , ( 17
)
where the subscripts 1 and 2 refer to the principal values of the tensors. Figure 11 shows the evolution of a n and a t with ε q in samples S1 and S2. We see that, in contrast to fabric anisotropies (Fig. 9), the force anisotropies in pentagon packing remain always above those in the disk packing. This means that the aptitude of the pentagon packing to develop large force anisotropy and strong force chains is not solely dependent on the global fabric anisotropy of the system. In section VI, we will show that the force anisotropy of the pentagon packing stems from the high anisotropy of the sub-network of double contacts and strong activation of friction forces. Indeed, due to the geometry of the pentagons, i.e. the absence of parallel sides, the strong force chains are mostly of zig-zag shape, as observed in Fig. 13b, and the stability of such structures requires strong activation of tangential forces. This explains, in turn, the large value of a t for pentagons, very close to a n , whereas in the disk packing a t is nearly half of a n .

The anisotropies a, a n and a t are interesting descriptors of granular microstructure and force transmission as they underlie the shear stress. Indeed, it can be shown that the general expression of the stress tensor Eq. ( 5) leads to the following simple relation [13,[START_REF] Radjai | Granular Materials: Fundamentals and Applications[END_REF]:

q p ≃ 1 2 (a + a n + a t ), ( 18
)
where the cross products between the anisotropies have been neglected and it has been assumed that the stress tensor is coaxial with the fabric tensor Eq. ( 11) and the q/p, S1 q/p, S2 0.5 (a+a n +a t ), S1 0.5 (a+a n +a t ), S2

FIG. 12: Evolution of the normalized shear stress q/p for the samples S1 and S2 with εq together with the corresponding predictions from its expression as a function of fabric and force anisotropies [Eq. ( 18)].

force tensors Eq. ( 16). Fig. 12 shows that Eq. ( 18) holds quite well both for pentagons and disks.

A remarkable consequence of Eq. ( 18) is to reveal the distinct origins of shear stress in pentagon and disk packings. The fabric anisotropy provides a major contribution to shear stress in the disk packing (Fig. 9) whereas the force anisotropies are more important for shear stress in the pentagon packing (Fig. 11). In this way, in spite of the weak fabric anisotropy a, the larger force anisotropies a n and a t allow the pentagon packing to reach higher levels of q/p compared to the disk packing.

B. Force distributions

The strong inhomogeneity of contact forces is a wellknown feature of granular media. It has been investigated mostly for spherical or cylindrical particles both by experiments and numerical simulations [1,3,4,5,6,7,8,9]. The probability density function (pdf) of normal forces is characterized by two features which seem to be specific to granular media: 1) The pdf is roughly a decreasing exponential function for forces above the mean, 2) In the range of weak forces below the mean, the pdf does not decline to zero with the force. The relative scatter of data reported by different authors for weak forces shows the sensitivity of the pdf in this range to the details of the microstructure. But, the common observation that there is a large number of contacts transmitting very weak forces, is a straightforward signature of the arching effect. From this point of view, one expects that angular particle shape will influence mainly the distribution of weak forces by enhancing the arching effect.

Figure 13 displays maps of normal forces in a portion of each of the samples S1 and S2 at a large cumulative strain. We observe the strong anisotropy of normal forces in the pentagon packing compared to the disk packing (as discussed in section V) as well as the zig-zag form of the strong force chains. The normal force pdf's are shown in

f n ∝ e -α1fn/ fn in S1, f n ∝ e -α2fn/ fn in S2, (19)
with α 1 ≃ 0.74 and α 2 ≃ 1.4. The smaller value of α 1 means that the distribution is wider for pentagons compared to disks. The distribution is nearly uniform in the whole range of weak forces (f n < f n ) in S2. In the pentagon packing S1, we observe a uniform distribution only in the range 0.1 f n < f n < f n . Nearly 30% of forces are in this range. The number of "very weak" forces in S1 in the range f n < 0.1 f n increases faster than a power law as f n tends to zero. A fraction ≃ 30% of contacts belong to this range. The presence of numerous "very weak" forces in the pentagon packing is a clear signature of enhanced arching effect that can be characterized, as we shall see below, by the respective roles of simple and double contacts with respect to force transmission.

C. Bimodal character of stress transmission

The genuine organization of contact forces in granular media, involving strong force chains propped by weak forces, was first analyzed by Radjai et al. by means of contact dynamics simulations for packings of circular and spherical particles [10]. This analysis proceeds by considering the subset of contacts which carry a force below a cutoff force ξ normalized by the mean force. This subset is referred to as the "ξ-network". The variation of a quantity evaluated for the "ξ-network" as ξ is varied from 0 to the maximal force in the system, provides its correlation with the contact force. Here, we apply this same approach to S1 and S2 samples for the stress ratio q(ξ)/p, defined as stress deviator q(ξ) (normalized by the total pressure p of the sample) in the ξ-network, and for a(ξ), defined as the fabric anisotropy in the ξ-network.

The plot of q(ξ)/p is shown in Fig. 15 for S1 and S2 in the residual state. In both samples, the stress deviator is nearly zero for ξ < 1, i.e. for the normal forces below the average force. This means that the shear stress is almost totally sustained by the "strong" contact network ξ > 1 for the pentagon packing as well as for the disk packing. Fig. 16 shows the fabric anisotropy a ′ (ξ) as a function of ξ in the samples S1 and S2. By definition, a positive value of a ′ corresponds to the principal stress direction whereas a negative value corresponds to the orthogonal direction. We see that the direction of anisotropy is orthogonal to the principal stress direction (a ′ < 0) for weak forces (small ξ). This "orthogonal" anisotropy of the weak forces is more important in the pentagon packing compared to the disk packing, and, as shown in the inset to Fig. 16, it is mainly due to "very weak" forces. When ξ is increased beyond f n , a ′ becomes less negative and finally changes sign, showing that the strong contacts are preferentially parallel to the principal axis. These strong contacts are less than 40% of all contacts, but their positive contribution to a ′ overcompensates the negative contribution weak contacts. For large ξ, the partial anisotropy approaches the fabric anisotropy of the whole system. These data demonstrate the bimodal character of stress transmission also in the pentagon packing in spite of a very different particle geometry. The mean force plays a particular role in differentiating strong contacts from weak contacts. However, the force pdf's (Fig. 14) and the anisotropy of weak forces (Fig. 16) provide also evidence for the existence of a class of very weak forces, corresponding approximately to the range f n < 0.1 f n , within the weak network. This class is strongly anisotropic with a privileged direction which is orthogonal to the major principal stress direction, and the corresponding force pdf diverges as the force tends to zero. Fig. 17 displays a tricolor map of the contact network representing very weak, intermediate (0.1 f n < f n < f n ) and strong contacts in the pentagon packing. Large cells of strong contacts are composed of zig-zag chains. The anisotropy of strong contacts is reflected in the elongated shape of these cells along the major principal stress direction. Both intermediate and very weak forces prop these cells.

VI. SIMPLE VERSUS DOUBLE CONTACTS

In this section, we focus on the organization of simple and double contacts in the pentagon packing. The double contacts, i.e. the side-sharing polygons, are generally assumed to be at the source of the higher strength of polygon packings. For the texture, we would like also to investigate the proportions of simple and double contacts and their respective contributions to the overall anisotropy of the pentagon packing. It is also important to identify the role of double contacts in force transmission. The general expression [Eq. 5] of the stress tensor σ allows us to perform a unique additive decomposition of the stress into two parts:

σ = σ s + σ d , (20)
where σ s is obtained from the expression (5) by restricting the summation to simple contacts, and σ d is the complementary tensor involving only double contacts. The respective stress deviators q s and q d normalized by the mean stress p are shown in Fig. 18 as a function of strain ε q . The strength q d /p of double contacts varies from two to three times that of simple contacts during shear deformation of the pentagon packing. We see that k d increases with strain but remains below k s . On the other hand, initially we have f d = f s = 0.5, reflecting the isotropic state of the packing prepared by isotropic compaction. However, f d increases with shear up to f d ≃ 1.5f s in the residual state. This means that the larger shear stress carried by double contacts in the residual state is due to the larger mean normal force of double contacts despite their smaller proportion in the packing.

The growth of the number of double contacts shown in Fig. 19 represents the gradual consolidation of the sample. In Fig. 20 we plot the cumulative proportions ∆γ s→d and ∆γ d→s of simple contacts turning to double and vice versa, respectively. Although transformation between the two contact types occurs at each step in both directions s → d and d → s, the consolidation involves on average a net fraction of simple contacts transforming into double contacts.

The connectivity of the pentagon packing by simple and double contacts can be represented by the proportion P (m s , m d ) of particles with exactly m s simple contacts and m d double contacts. Fig. 21 shows a grey level map of this function for the pentagon packing in the residual state. The row m d = 0 corresponds to particles with only simple contacts (nearly 2% of the total number of particles) whereas the columnm s = 0 represents the particles with only double contacts (nearly 6%). On average, a particle has more simple contacts than double contacts but the maximum occurs at m s = m d = 2.

We now consider the fabric tensor decomposed in a similar way as the stress tensor [Eq. [START_REF] Pena | Granular Matter In Press[END_REF]] into two partial tensors:

F = F s + F d , (21)
where F s and F d are defined as F in Eq. ( 11) by simply restricting the summation to simple and double contacts, respectively, and by dividing the sum by the total number N c of contacts. The respective anisotropies a ′ s and a ′ d of simple and double contacts are displayed in Fig. 22 as a function of ε q . The interesting observation here is that the simple contacts have a negative anisotropy which, according to Eq. ( 13), means that simple contacts are mostly oriented perpendicular to the major principal fabric direction θ F . In other words, most simple contacts belong to the weak network. In contrast, the double contacts have an increasing positive anisotropy which is larger than the mean anisotropy a of the sample. This is consistent with the fact that the double contacts take over larger forces and they contribute more to the shear stress than simple contacts.

The normal force pdf's for simple and double contacts are shown in Fig. 23. Both contact types are involved in weak and strong networks and the pdf's have the same functional form. But the contribution of simple contacts is more important in the range of very weak forces. Once again, as for anisotropy, the very weak contacts appear to be related to the particular geometry of the pentagons. At large strains, about 32% of all contacts belong to the very weak force network with 25% simple contacts against 7% double contacts. A snapshot of the normal force network is shown in Fig. 24 where the line widths are proportional to the line width with different colors (or grey levels) for simple and double contacts. The remarkable feature of this map is the network of very strong zigzag force chains composed mostly of double contacts and occasionally mediated by simple contacts.

The proportions k S s and k W s of strong (S) and weak (W) simple (s) contacts, respectively, as well as the proportions k S d and k W d of strong and weak double (d) contacts are plotted in Fig. 25 as a function of ε q . We see that in the strong network (f n > f ) the proportion k S d of double contacts is nearly the same as the proportion k S s of simple contacts in the initial (isotropic) state, but during shear k S s declines down to k S s ≃ 0.5k S d in the residual state, in agreement with the impression left by Fig. 24. We have an inverse situation for the weak network composed of two times more simple contacts than double contacts, i.e. k W s ≃ 2k W d in the residual state. It is also interesting to remark that the fraction of weak contacts, i.e. k W s + k W d ≃ 0.58 in the residual state is very close to that (0.62) in the case of the disk packing.

VII. CONCLUSION

The objective of this paper was to isolate the effect of particle shape on force transmission in granular media by means a detailed comparison between two similar packings with different particle shapes: pentagons vs. disks. We observed enhanced shear strength and force inhomogeneity in the pentagon packing. But, unexpectedly, the pentagon packing was found to develop a lower structural (fabric) anisotropy compared to the disk packing under shear. This low fabric anisotropy, however, does not prevent the pentagon packing from building up a strong force anisotropy that underlies its enhanced shear strength compared to the disk packing.

This finding is interesting as it shows unambiguously that the force anisotropy in a granular material has two distinct sources: (1) Fabric anisotropy, with a maximum value depending on particle shapes; (2) Particle shapes. The first mechanism is crucial for the disk packing so that the force anisotropy, and the shear stress as a result, vanishes in an isotropic disk packing (e.g. when the friction coefficient is set to zero). The second mechanism may be the predominant source of strength for "facetted" particles that can give rise edge-to-edge (in 2D) or faceto-face (in 3D) contacts allowing for strong force localization along such contacts in the packing. Since the fabric anisotropy is low in a pentagon packing, the role of force anisotropy and thus the local equilibrium structures or arching are important with respect to its overall strength properties. The pentagons analyzed in this work provide thus the first counter-example of a system where the role of fabric anisotropy in shear strength is marginal.

Another shape-related effect was the observation of zigzag force chains mostly composed of edge-to-edge contacts in steady shearing. The vertex-to-edge contacts belong thus mainly to the weak force network or a class of "very weak" forces that can be considered as a signature of enhanced arching or screening effect of forces in the presence of edge-to-edge chains. These "very weak" forces can also be observed, though to a lesser extent, in a disk packing with high coefficients of friction [8] or on experimental pdf's of normal forces acting on the walls of a container [5]. Let us recall that a "very weak phase" was also evidenced by considering the correlation between friction mobilization and the anisotropy of granular texture in a disk packing at the stability limit [14].

By focusing on pentagon packings, we were able to demonstrate the nontrivial phenomenology resulting from the specific shape of particles as compared to a disk packing. Although general features of force transmission (pdf's, bimodal character, etc) seem to be robust, the details of force transmission (relative importance of force and fabric anisotropy, the role of edge-to-edge contacts, etc) seem to be strongly shape-dependent. Cur-rently, we work to elucidate this issue for regular polygons (hexagons and higher number of sides) as well as polyhedral particles in three dimensions. We warmly thank Fréderic Dubois for assistance with the LMGC90 platform used for the simulations. This

FIG. 1 :

 1 FIG.1: Representation of simple (vertex-to-edge) and double (edge-to-edge) contacts between two pentagons.

FIG. 2 :

 2 FIG. 2: Snapshots of a portion of the samples S2 (a) and S1 (b) composed of circular and pentagonal particles, respectively.

FIG. 4 :FIG. 5 :

 45 FIG.4: Cumulative volumetric strain εp as a function of cumulative shear strain εq for the samples S1 and S2.

FIG. 6 :FIG. 7 :

 67 FIG.6: The coordination number z (a) and the fraction γ of persistent contacts (b) as a function of cumulative shear strain εq for the samples S1 and S2.

FIG. 9 :

 9 FIG. 9: Evolution of the anisotropy a ′ with cumulative shear strain εq for the samples S1 and S2.

FIG. 10 :

 10 FIG. 10: Polar representation of the angle-averaged normal (a) and tangential (b) forces fn (θ) and ft (θ) for the samples S1 and S2 in the residual state.

FIG. 11 :

 11 FIG. 11: Evolution of force anisotropies an (a) and at (b) as a function of cumulative shear strain εq in samples S1 and S2.

FIG. 13 :

 13 FIG. 13: (color online) Snapshots of normal forces in samples S2 (a) and S1 (b). Line thickness is proportional to the normal force.

Fig. 14

 14 Fig. 14 in log-linear and log-log scales at large strains. The forces are normalized by the mean normal force f n in each sample. In both samples, the number of strong forces (above the mean f n ) falls off exponentially:

FIG. 14 :

 14 FIG. 14: Probability density functions of normal forces in samples S1 and S2 in log-linear (a) and log-log scales (b).

FIG. 15 :FIG. 16 :

 1516 FIG.15: Partial shear stress q(ξ)/p as a function of force cutoff ξ (normalized by the mean force) for the samples S1 and S2 in the residual state.

FIG. 17 :FIG. 18 :

 1718 FIG. 17: (color online) Tricolor map of the contact network composed of very weak (blue), intermediate (green) and strong (red) contacts in the pentagon packing.

FIG. 19 :

 19 FIG. 19: Proportions ks and k d of simple and double contacts, and the corresponding relative force averages fs and f d , as a function of cumulative shear strain εq.

  The proportions k s and k d of simple and double contacts are shown in Fig. 19 as a function of ε q . The same figure displays the relative force averages f s = k s f n d / f n and f d = k d f n s / f n , where f n s and f n d are the mean normal forces of simple and double contacts, respectively.

FIG. 20 :

 20 FIG. 20: Cumulative proportion ∆γ of simple contacts turning to double (s → d) and vice versa (d → s).

FIG. 21 :

 21 FIG. 21: Grey level map of the connectivity function P (ms, m d ) of the pentagon packing in the residual state.

FIG. 22 :

 22 FIG.22:The anisotropy a ′ of simple (s) and double (d) contacts as a function of cumulative shear strain εq in the pentagon packing.

FIG. 23 :

 23 FIG. 23: Probability density function of normal forces for simple (s) and double (d) contacts in log-linear (a) and loglog scales (b).

FIG. 24 :

 24 FIG. 24: (color online) Color map of the normal force network in the residual state with simple contacts (s) in blue and double contacts (d) in red. Line thickness is proportional to the normal force.

FIG. 25 :

 25 FIG. 25: Proportions k S s and k W s of strong (S) and weak (W) simple (s) contacts, respectively, as well as the proportions k S d and k W d of strong and weak double (d) contacts as a function of cumulative shear strain εq in the pentagon packing.

work was funded by the French Railway Society, the SNCF, and the Région Languedoc-Roussillon of France.


David Néron
email: neron@elmt.ens-cachan.fr

Hachmi Ben Dhia

Régis Cottereau

A decoupled strategy to solve reduced-order multimodel problems in the PGD and Arlequin frameworks

Keywords: reduced-order modeling, PGD, multimodel, multiscale, Arlequin, LATIN

In this paper, we investigate the coupling of reduced models for the simulation of structures involving localized geometrical details. Herein, we use the Arlequin method, originally designed to deal with multimodel and multiscale analyses of mechanical problems, to mix reduced models built using the Proper Generalized Decomposition. Instead of solving the global coupled problem in a monolithic way, the LATIN strategy is used to propose a decoupled algorithm. The numerical examples demonstrate the feasibility of the approach and in particular its potentiality in terms of flexibility.

Introduction

In engineering design, the numerical simulation of very large multiscale models is becoming increasingly important because of the need to describe realistic scenarios and to derive tools to facilitate the virtual design of new structures. However, the incredible evolution of computing resources over the last decade is balanced by the increasing complexity of the models that engineers want to address in their efforts to design, control and optimize innovative products. Solving several problems with very large number of degrees of freedom, with the presence of multiscale and multiphysics aspects, or with the need to take into account uncertainties or parameter variations cannot be handled by standard solution techniques. In this context, model reduction methods have a huge potential to develop innovative tools for intensive computation and allow a "real time" interaction between the user and the simulations, which gives the opportunity to explore a large number of scenarios in the design office.

Several approaches of reduced-order modeling have been proposed, many of them relying on the assumption of a separated form of the unknowns. A first family of algorithms starts by a learning phase (performing some preliminary computations, called snapshots) to build a reduced-order basis that allows to capture the principal characteristics of the solution of the problem to be solved. It is the case of the methods based on the Proper Orthogonal Decomposition (see e.g. [START_REF] Maday | The reduced-basis element method: application to a thermal fin problem[END_REF][START_REF] Barrault | An "empirical interpolation" method: Application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Lieu | Reduced-order fluid/structure modeling of a complete aircraft configuration[END_REF][START_REF] Gunzburger | Reduced-order modeling of timedependent pdes with multiple parameters in the boundary data[END_REF]) or of the Reduced Basis method (see e.g. [START_REF] Patera | Reduced Basis Approximation and A Posteriori Error Estimation for Parametrized Partial Differential Equations[END_REF][START_REF] Rozza | On the stability of the reduced basis method for Stokes equations in parametrized domains[END_REF]). The latter adds an automatic selection of the snapshots by a greedy algorithm based on some efficient error indicators. Another family follows a different path as it builds progressively an approximate separated representation of the solution, without assuming any basis or selecting any snapshots. It is the case of the Proper Generalized Decomposition (PGD, see e.g. [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics-New Approaches and Non-Incremental Methods of Calculation[END_REF][START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids: Part II: Transient simulation using space-time separated representations[END_REF] or [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF] for a review of the method) which was considered in our former works, in particular for the analysis of elastic-viscoplastic problems [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF], multiscale problems [START_REF] Cremonesi | A PGD-based homogenization technique for the resolution of nonlinear multiscale problems[END_REF] or parametrized studies [START_REF] Néron | Time-space PGD for the rapid solution of 3D nonlinear parametrized problems in the many-query context[END_REF]. This method has also been widely developed by Chinesta and coauthors who proposed very efficient implementations of the approach and a large variety of applications (see [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF] for an overview). To focus only on some recent examples, one can cite real-time simulations in surgery [START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF][START_REF] Alfaro | Real-time in silico experiments on gene regulatory networks and surgery simulation on handheld devices[END_REF], real-time monitoring of thermal process [START_REF] Aguado | Real-time monitoring of thermal processes by reduced-order modeling[END_REF] or the simulation of viscoelastic models [START_REF] Ammar | On the space-time separated representation of integral linear viscoelastic models[END_REF].

However, to move a step forward, a major issue which is addressed herein is the coupling of reduced-order models to mix several models within the same simulation. This point is exemplified by Figure 1, which represents a holed square structure with different shapes. The brute force classical approach to evaluate the influence of the shape of the hole is to mesh the structure for each geometry and run the simulation, possibly using a reduced modeling strategy that generates modes specific to this particular geometry.

The approach which is explored in this paper is different. We propose to superpose two reduced models: one for the global structure and one for the local detail (see Figure 2). This would allow a considerable increase in the flexibility of the solution and provide new opportunities in terms of design, such as:

• the ability to consider a complete system as an assembly of components, each represented by its own model;

• the possibility of handling reduced models arising from several actors involved in the design of the same product;

• the ability to take into account the critical phenomena, which occur locally in time and space, with a finer numerical model that could be adapted automatically during the simulation;

• the flexibility to change some local data (geometry, topology, material, architecture ...) in order to change the global structure. The issue of coupling reduced models is essentially the possibly different natures of the mathematical models to "marry", that is the management of incompatibilities between the models. This question was addressed in [START_REF] Néron | A computational strategy for poroelastic problems with a time interface between coupled physics[END_REF][START_REF] Néron | A computational strategy for thermo-poroelastic structures with a time-space interface coupling[END_REF] for multiphysics problems such as thermo-poroelasticity. In these works, based on the concept of "interface between physics", the LATIN method was used to solve in a decoupled manner the problems corresponding to the different physics, and the PGD was applied to reduce the computational cost of the models.

In this paper, we propose to use the Arlequin technique [START_REF] Dhia | Multiscale mechanical problems: the Arlequin method (in French)[END_REF]. This technique has been designed to deal with problems in which several "zones" of interest can be distinguished and require different levels of analysis. The term "zones" should be understood in the broad-sense, as it concerns different numerical models whose fields can be mixed and glued together. For that, a superposition technique based on a weak formulation, in which the energy is distributed between the various models is used. This technique allows to deal with incompatible models, including the ones defined on different meshes. The Arlequin method has been the subject of many developments (see e.g. [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF][START_REF] Cottereau | A stochasticdeterministic coupling method for continuum mechanics[END_REF]) and showed its capabilities to locally refine models to introduce an essential local modification in models.

The Arlequin modeling framework has already been used in [START_REF] Nazeer | Arlequin based PGD domain decomposition[END_REF] in the context of the PGD. In this work, the technique was implemented to split a complex physical geometry in rather simple ones, in order to allow a separation of the geometrical coordinates. The different simple zones were then solved using the PGD approximation introduced in the global formulation, allowing a strong decrease of the computational cost.

Herein, the aim is different. It consists in transforming the classical weak-Arlequin volume coupling into a strong (pointwise) coupling. This local coupling between models allows for the use of the LATIN algorithm as a decoupled solver and the implementation of the PGD technique for each of the overlapping models. This approach leads to an interesting flexibility as it will permit at the end, to use different PGD solvers, dedicated to the specificities of the models that are considered (linear/nonlinear, deterministic/stochastic, atomistic/continuum ...) as well as, in the future, resorting to precomputed reduced-models stored in a library.

The present paper is organized as follows. In Section 2, the framework of the Arlequin strategy is recalled in the case of linear elasticity and reformulated to be solved by a decoupled algorithm. In Section 3, the LATIN algorithm is introduced to solve the problem in a decoupled manner. The PGD approach is then used to solve the problems corresponding to the different models. Finally, some numerical examples are proposed in Section 4 to illustrate the behavior of the strategy.

Arlequin formulation

Only the main features of the approach are recalled herein. For further details, especially on the theoretical aspects, the interested reader can refer to [START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF]. In order to recall the Arlequin formulation, we consider the following representative model.

Classical linearized elasticity problem

Let us consider the evolution over the time interval I = (0, T ) of an elastic body occupying the closure of a bounded regular domain Ω 1 included in R d , with d = 1, 2 or 3 in practice. It is submitted to a field of volume density of forces f , a field of surface density of forces F on Γ F , a part of its boundary, and clamped on Γ u , a non zero measured other part of the boundary (see Figure 3).

f F Ω 1 Γ F Γ u Figure 3: Reference problem
For a given space V , we denote L 2 (I; V ) = {v :

I → V ; I ||v|| 2 V dt < +∞} and ||•|| V a norm on V . Then, we introduce V 0 = {v ∈ H 1 (Ω 1 ) ; v = 0 on Γ u } and we define W 0 = L 2 (I; V 0 ), F 0 = L 2 (I; L 2 (Ω 1 )) and G 0 = L 2 (I; L 2 (Γ F )).
If we assume that f ∈ F 0 and F ∈ G 0 , the displacement field u 0 is a function of (x, t) ∈ Ω 1 × I and the weak primal (monomodel) formulation of this problem reads as follows:

Find u 0 ∈ W 0 ; ∀ v 0 ∈ W 0 , a 0 (u 0 , v 0 ) = 0 (v 0 ) (1)
where the virtual works of the internal and external forces respectively read:

∀ (u 0 , v 0 ) ∈ W 0 × W 0 , a 0 (u 0 , v 0 ) = Ω1 σ(u 0 ) : ε(v 0 ) dΩ (2) ∀ v 0 ∈ W 0 , 0 (v 0 ) = Ω1 f • v 0 dΩ + Γ F F • v 0 dS (3)
In these equations, ε(v 0 ) and σ(v 0 ) denote the linearized strain and stress tensors, associated to field v 0 , which are assumed to be connected through a Hooke's law.

Lagrange multiplier-based Arlequin formulations of the elasticity problem

Let Ω 2 be a non-zero measured given bounded regular domain overlapping Ω 1 . Let Ω 12 be the overlap. For clarity and with no major restrictions, it will be assumed that Ω 2 is strictly embedded in Ω 1 , leading to Ω 12 = Ω 2 . The overlap where Ω c 12 is the models gluing zone (see the grey zone in Figure 4). Now, to define Arlequin formulations of the elasticity problem, we denote by

Ω 1 Ω 1 Ω 2 Ω 2
W 1 = L 2 (I; V 1 ) with V 1 = {v ∈ H 1 (Ω 1 ) ; v = 0 on Γ u } and W 2 = L 2 (I; V 2 ) with V 2 = H 1 (Ω 2 )
. Weighted internal and external virtual works are defined by:

∀ (u 1 , v 1 ) ∈ W 1 , a 1 (u 1 , v 1 ) = Ω1 α 1 σ(u 1 ) : ε(v 1 ) dΩ (4) ∀ (u 2 , v 2 ) ∈ W 2 , a 2 (u 2 , v 2 ) = Ω2 α 2 σ(u 2 ) : ε(v 2 ) dΩ (5) ∀ v 1 ∈ W 1 , 1 (v 1 ) = Ω1 α 1 f • v 1 dΩ + Γ F α 1 F • v 1 dS (6) ∀ v 2 ∈ W 2 , 2 (v 2 ) = Ω2 α 2 f • v 2 dΩ ( 7
)
where weight parameter functions (α 1 , α 2 ) are defined respectively in Ω 1 and Ω 2 and satisfy (see Figure 5) for i = 1, 2:

α i 0 in Ω i , α i = 1 in Ω i \ Ω 12 , α 1 + α 2 = 1 in Ω 12 (8)
The dual volume coupling-based continuous Arlequin formulation of the model elasticity problem reads: [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF] for details on the choice of these functions and the value of ε ≈ 10 -2 , which allows to avoid the ill-conditioning of the operators of the method)

Find (u 1 , u 2 , Φ d ) ∈ W 1 × W 2 × M d ; ∀ v 1 ∈ W 1 , a 1 (u 1 , v 1 ) + c d (Φ d , v 1 ) = 1 (v 1 ) (9) ∀ v 2 ∈ W 2 , a 2 (u 2 , v 2 ) -c d (Φ d , v 2 ) = 2 (v 2 ) (10) ∀ Ψ d ∈ M d , c d (Ψ d , u 1 -u 2 ) = 0 (11) 0 ε ε 1- 1 α α 2 1 α 1 Ω2 = Ω12 Ω Ω1 12 f Figure 5: Example of weight parameter functions α i , (i = 1, 2) (see
where M d is the dual mediator space of M = W 1 /Ω c 12 = W 2 /Ω c 12 = L 2 (I; V c 12 ) with V c 12 = H 1 (Ω c 12 ) and c d (• , •)
is the volume coupling operator:

∀ (Ψ d , v) ∈ M d × M , c d (Ψ d , v) = Ψ d , v M d , M (12)
where • , • M d , M stands for the duality bracket. The dual volume coupling is a natural mechanical coupling operator, in the sense that, interpreting the Lagrange multiplier field Φ d as a density of forces, it has to be in the dual space M d of the space of the displacements in Ω c 12 .

By using the Riesz-Fréchet representation theorem, the natural H 1 -scalar product of the space M can be substituted to the duality bracket which leads to:

Find (u 1 , u 2 , Φ) ∈ W 1 × W 2 × M ; ∀ v 1 ∈ W 1 , a 1 (u 1 , v 1 ) + c(Φ, v 1 ) = 1 (v 1 ) (13)
∀ v 2 ∈ W 2 , a 2 (u 2 , v 2 ) -c(Φ, v 2 ) = 2 (v 2 ) (14)
∀ Ψ ∈ M , c(Ψ, u 1 -u 2 ) = 0 ( 15
)
where the coupling operator, denoted by c(• , •), is defined by:

∀(Ψ , v) ∈ M × M , c(Ψ, v) = Ω c 12 κ(ε(Ψ) : ε(v) + 1 e 2 Ψ • v)dΩ ( 16
)
where κ is a positive parameter of the order of magnitude of the rigidity of the material (typically the Young modulus E) and e is homogeneous to a length (typically e c 12 , the width of the gluing zone, see Figure 4). The choice of the norm and of this width is not at the core of this paper and the guidelines of [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Dhia | Further insights by theoretical investigations of the multiscale arlequin method[END_REF] will be followed. Lagrange multiplier Φ ∈ M appearing in [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF][START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF] is homogeneous to a displacement field.

After discretization in space, Problem [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF][START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF][START_REF] Alfaro | Real-time in silico experiments on gene regulatory networks and surgery simulation on handheld devices[END_REF] can be rewritten using ob-vious notations:

∀t ∈ I,   K 1 0 C T 1 0 K 2 -C T 2 C 1 -C 2 0     u 1 (t) u 2 (t) Φ(t)   =   F 1 (t) F 2 (t) 0   (17)
The solution of ( 17) is denoted s ref .

It is usually obtained using a "monolithic" approach which consists in solving directly the previous system that couples the models. Despite the simplicity of implementation, the monolithic approach is not the most appealing approach as it does not allow to use appropriate solvers for each of the super-imposed models. In particular, one can think of the coupling between models such as linear/nonlinear, deterministic/stochastic, atomistic/continuum ... That is why some works have been proposed in [START_REF] Ben Dhia | Multimodeling of multialterated structures in the arlequin framework. solution with a domaindecomposition solver[END_REF] to circumvent this issue by mean of the FETI approach. In the present paper, we propose to use another solver, the LATIN method, to ensure the decoupling of the solution process but also to use reduced-order modeling to solve more efficiently the various models. For that purpose, the Arlequin problem is reformulated as follows.

Reformulation of the Arlequin problem

Problem [START_REF] Chinesta | A short review on model order reduction based on proper generalized decomposition[END_REF][START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF][START_REF] Cremonesi | A PGD-based homogenization technique for the resolution of nonlinear multiscale problems[END_REF] is rewritten by introducing new dual Lagrange multipliers (Φ d1 , Φ d2 ) ∈ M 2 d in the formulation, such that (9,10) reads:

∀v 1 ∈ W 1 , a 1 (u 1 , v 1 ) -c d (Φ d1 , v 1 ) = 1 (v 1 ) ∀v 2 ∈ W 2 , a 2 (u 2 , v 2 ) -c d (Φ d2 , v 2 ) = 2 (v 2 ) ( 18
)
provided that we enforce:

Φ d1 + Φ d2 = 0 (19)
which can be interpreted as the equilibrium in the volume Ω c 12 of the densities of forces. Coming back to the primal formulation [START_REF] Chinesta | The Proper Generalized Decomposition for Advanced Numerical Simulations: a Primer[END_REF][START_REF] González | Computational vademecums for the real-time simulation of haptic collision between nonlinear solids[END_REF][START_REF] Alfaro | Real-time in silico experiments on gene regulatory networks and surgery simulation on handheld devices[END_REF], new primal Lagrange multipliers (Φ 1 , Φ 2 ) ∈ M 2 are introduced such that (13,14) can be rewritten:

∀v 1 ∈ W 1 , a 1 (u 1 , v 1 ) -c(Φ 1 , v 1 ) = 1 (v 1 ) ∀v 2 ∈ W 2 , a 2 (u 2 , v 2 ) -c(Φ 2 , v 2 ) = 2 (v 2 ) (20)
with:

Φ 1 + Φ 2 = 0 (21)
We also introduce new displacement fields (w

1 , w 2 ) ∈ M 2 such that ∀Ψ ∈ M , c(Ψ, w i ) = c(Ψ, u i ) (22)
Equation ( 22) allows to define the restriction operators Π i , (i = 1, 2), that transfer a field defined on space W i to the mediator space M :

Π i : u i ∈ W i → w i = Π i u i ∈ M (23)
Using these new fields, equation ( 15) which corresponds to the accommodation between the two models can simply be rewritten as the equality of w 1 and w 2 :

w 1 -w 2 = 0 (24)
Then, the Arlequin formulation that will be considered in the next section is:

Find i = (1, 2) (u i , w i , Φ i ) ∈ W i × M × M ; ∀v i ∈ W i , a i (u i , v i ) -c(Φ i , v i ) = i (v i ), w i = Π i u i ( 25
)
w 1 -w 2 = 0 and Φ 1 + Φ 2 = 0 (26)
3. Decoupled resolution using the LATIN method

The LATIN method as a decoupled solver

Problem [START_REF] Ben Dhia | Multimodeling of multialterated structures in the arlequin framework. solution with a domaindecomposition solver[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF] is solved using the LATIN framework. Only the basics of the method are presented in the following and the interested reader can refer to [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF] for coupling material subdomains or [START_REF] Néron | A computational strategy for thermo-poroelastic structures with a time-space interface coupling[END_REF] for coupling multiphysics models. Equation ( 26) plays the role of an "interface" between the models, which is similar to the interfaces between subdomains and between physics used in these works.

To find the solution [START_REF] Ben Dhia | Multimodeling of multialterated structures in the arlequin framework. solution with a domaindecomposition solver[END_REF][START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF], the idea is a partitioning of the equations by introducing two subsets of elements of

s ref = {(u i , w i , Φ i )} i=1,2 of problem
W 1 × M × M × W 2 × M × M .
The first one is denoted A d and corresponds to the set of solutions to problem [START_REF] Ben Dhia | Multimodeling of multialterated structures in the arlequin framework. solution with a domaindecomposition solver[END_REF]. The second, denoted Γ corresponds to the set of solutions to problem [START_REF] Ladevèze | The LATIN multiscale computational method and the Proper Generalized Decomposition[END_REF]. The solution s ref = A d ∩ Γ is build using an alternated-direction-scheme which consists in finding alternatively an element of Γ and an element of A d . The advantage of this partitioning is mainly that A d involves equations that are independent between the two models and Γ involve equations that couple the two models but which are defined locally in the accommodation space. This algorithm is sketched in Figure 6.

At iteration n + 1, an element s n ∈ A d is assumed to be already computed. The algorithm consists of two stages, called coupled and decoupled. These stages lead to ŝn+1 ∈ Γ and s n+1 ∈ A d and work as follows (subscript i, corresponding to the two models, has to be interpreted as i = 1, 2; subscript n + 1 is skipped to simplify the notations):

• • • -→ s n ∈ A d coupled stage --------→ ŝn+1 ∈ Γ decoupled stage ----------→ s n+1 ∈ A d Iteration n + 1 -→ ŝn+2 -→ • • •
Coupled stage. -Knowing a solution s ∈ A d , the coupled stage consists in finding a solution ŝ ∈ Γ such that a linear search direction E + is fulfilled:

Φi -Φ i = k i ( ŵi -w i ) ( 27
)
where k i ∈ R + are parameters of the method that will be discussed in the following. Recalling the compatibility conditions on the "interface":

ŵ1 -ŵ2 = 0 and Φ1 + Φ2 = 0 (28)
the solution at this stage is defined explicitly by:

ŵi = 1 k 1 + k 2 (ϕ 1 + ϕ 2 ) Φi = k i ŵi -ϕ i ( 29
)
where ϕ i = k i w i -Φ i are known quantities.

Decoupled stage. -Knowing a solution ŝ ∈ Γ, the decoupled stage consists in finding a solution s ∈ A d such that a linear search direction E -is fulfilled:

Φ i -Φi = -k i (w i -ŵi ) (30)
Recalling the equilibrium equations of the models:

∀v i ∈ W i , a i (u i , v i ) -c(Φ i , v i ) = i (v i ), w i = Π i u i ( 31
)
and using search direction [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF] to substitute Φ i , one obtains:

∀v i ∈ W i , a i (u i , v i ) -c( Φ i -k i w i -αi = -k i Π i u i -αi , v i ) = i (v i ) ( 32
)
where αi = -k i ŵi -Φi are known quantities. One obtains:

∀v i ∈ W i , a i (u i , v i ) + c(k i Π i u i , v i ) = i (v i ) -c( αi , v i ) ( 33
)
and the decoupled problems to be solved at this stage can be rewritten formally:

Find u i ∈ W i ; ∀v i ∈ W i , a i (u i , v i ) = l i (v i ) (34) Then (w i , Φ i ) ∈ M × M ; w i = Π i u i and Φ i = -k i w i -αi (35)
where

a i (•, •) = a i (•, •) + c(k i Π i •, •) and l i (•) = i (•) -c( αi , •).

Some remarks on the algorithm

Initialisation.

-The previous algorithm is initialized by assuming that there is no link between the two models. For that purpose, a first decoupled stage is performed to find a solution s 0 ∈ A d by assuming:

ŵ1 = ŵ2 = 0 and Φ1 = -Φ2 = 0 ( 36
)
Convergence indicator. -The convergence of the algorithm is controlled by the indicator η LATIN computed at each iteration:

η LATIN = ||ŝ -s|| Ω c 12 || 1 2 (ŝ + s)|| Ω c 12 ( 37
)
with the norm ||s|| 2

Ω c 12 = ||w 1 || 2 k1 + ||w 2 || 2 k2 = I×Ω c 12 w 1 • k 1 w 1 dSdt + I×Ω c 12 w 2 • k 2 w 2 dSdt
, and the iterations are stopped when the indicator is less than a threshold η (in practice η = 10 -2 in the following numerical examples). Note that, following for example [START_REF] Ladevèze | Nonlinear Computational Structural Mechanics-New Approaches and Non-Incremental Methods of Calculation[END_REF], this indicator can also be linked to the error with respect to the reference monolithic solution s ref but this aspect is out of the scope of this first work on the subject and will be studied in some further works.

Search direction parameters. -In the LATIN method, the choice of the search direction parameters k i influences the convergence rate of the algorithm and the order of magnitude is chosen following the works on domain decomposition in this framework (see, e.g. [START_REF] Néron | Proper Generalized Decomposition for multiscale and multiphysics problems[END_REF] or [START_REF] Allix | Virtual delamination testing through non-linear multi-scale computational methods: Some recent progress[END_REF] for the case of nonlinearities). During the decoupled stage, the search direction, written in terms of the density of force Φ di and of the displacement field w i , should be written

Φ di -Φdi = -k di (w i -ŵi ),
where k di is of the order of magnitude of the rigidity of the material (typically the Young modulus) divided by the square of a characteristic length of the coupling area (typically the width of the gluing zone). This leads to k di which is chosen of the order of E/e c2 12 .

Recalling that the primal field Φ and the dual field Φ d are linked by ( 16), Φ d is of the order of κ/e 2 Φ. Finally, using the fact that κ is chosen equal to E and e to e c 12 , the search direction parameters used in the numerical examples will be k 1 = k 2 = 1. The influence of this value will be discussed in the following.

Discretization of the Arlequin formulation.

-Classical finite element discretizations are used for spaces W i with basis functions that are continuous and polynomial by parts over each of the elements. The meshes corresponding to the two models are not a priori compatible. A finite element discretization is also introduced for the mediator space M . In practice, the latter will often be a restriction of the coarser discretized space, and hence incompatible with the finer discretized mesh, but this is not a requirement. When considering these incompatible meshes, the integration of the term c(k i Π i •, •) in (33) is not classical because it involves functions that may not be polynomials over the elements of any of the meshes. Quadrature methods can therefore not be applied in a straightforward way. Several strategies can be used, for instance by adding adaptively Gauss points to the integration, as done in [START_REF] Ben Dhia | Application of the arlequin method to some structures with defects[END_REF]. In the present work, one prefers to construct the geometrical intersection of the two meshes, as done in [START_REF] Ben Dhia | The Arlequin method as a flexible engineering design tool[END_REF], so that one falls back on classical integration schemes.

The coupled stage (see Equation ( 29)) is solved in a strong form, while the decoupled stage (see Equation ( 34)) is solved in a weak form. The question of the stability of this scheme, in particular when the functional spaces W i are discretized over incompatible meshes, will be the subject of some theoretical investigations in further works and, herein, we limit ourselves to the apparent stability of the numerical experiments that will be described in Section 4.

Proper Generalized Decomposition at a glance

During the decoupled stage, the problem to be solved for each model i = (1, 2) is defined by (35), that is the space weak formulation:

Find u i ∈ W i ; ∀v i ∈ W i , a i (u i , v i ) = l i (v i ) (38)
where we recall that u i is a function of (x, t)

∈ Ω i × I, W i = L 2 (I; V i ) and V 1 = {v 1 ∈ H 1 (Ω 1 ) ; v 1 = 0 on Γ u } and V 2 = H 1 (Ω 2 )
. We also introduce I = L 2 (I, R) the space of square integrable functions on I.

An approximation of solution u i is sought using the Proper Generalized Decomposition (PGD) framework. The aim is to find a low-rank approximation of a space-time field as a sum of products of spatial and temporal functions. In our case, space W i is approximated by S i = V i ⊗ I, that is:

u i (x, t) ≈ u m i (x, t) = m k=1 Λ k i (x)λ k i (t) where (Λ k i , λ k i ) ∈ V i × I (39)
Many choices are possible to build the PGD (Galerkin orthogonality criteria, minimal residual criteria ... [START_REF] Nouy | A priori model reduction through proper generalized decomposition for solving time-dependent partial differential equations[END_REF]). The algorithm which is used herein is based on the Galerkin orthogonality and, for that purpose, a space-time weak formulation of (38) is introduced using the bilinear form A i and the linear form L i defined by:

A i (u i , v i ) = I a i (u i , v i )dt and L i (v i ) = I l i (v i )dt (40)
The new problem reads:

Find u i ∈ S i ; ∀v i ∈ T ui (S i ), A i (u i , v i ) = L i (v i ) ( 41
)
where T ui (S i ) is the tangent linear space to S i at u i . The construction of the solution is done using a greedy procedure. Assuming that a decomposition u m-1 i of order (m -1) is already known:

u m-1 i = m-1 k=1 Λ k i λ k i where (Λ k i , λ k i ) ∈ V i × I (42)
the solution of ( 41 of order (m -1) is known (42), a recombination of the reduced-order basis is searched by computing an approximation ūm-

1 i = u m-1 i + m-1 k=1 Λ k i λk
i where only the time functions λk i have to be searched. For that, a Galerkin orthogonality formulation of (41) is used:

Find { λk i } m-1 k=1 ∈ I m-1 ; ∀{λ l i } m-1 l=1 ∈ I m-1 , A i (u m-1 i + m-1 k=1 Λ k i λk i , m-1 l=1 Λ l i λ l i ) = L i ( m-1 l=1 Λ l i λ l i ) ( 43
)
This problem is rewritten:

Find { λk i } m-1 k=1 ∈ I m-1 ; ∀{λ l i } m-1 l=1 ∈ I m-1 , A i ( m-1 k=1 Λ k i λk i , m-1 l=1 Λ l i λ l i ) = L i ( m-1 l=1 Λ l i λ l i ) -A i (u m-1 i , m-1 l=1 Λ l i λ l i ) (44)
or, in a more compact form:

Find { λk i } m-1 k=1 ∈ I m-1 ; ∀{λ l i } m-1 l=1 ∈ I m-1 , A i ( m-1 k=1 Λ k i λk i , m-1 l=1 Λ l i λ l i ) = B m i ( m-1 l=1 Λ l i λ l i ) ( 45
)
where one introduces

B m i (•) = L i (•) -A i (u m-1 i , •).
It is straightforward to derive that (45) leads to a small system of order (m-1) which can moreover be decoupled provided that {Λ k i } m-1 k=1 is orthonormalized with respect to the scalar product associated to a i , that is a i (Λ k i , Λ l i ) = δ kl . The new approximation ūi ≈ ūm-1 i being known, the corresponding wi and Φi can be computed using the restriction operator and the search direction.

A criterion is then used to check if this recombination of the previous reduced-order basis was sufficient or if an enrichment is required. Reusing the subscripts introduced in Figure 6 for the consecutive iterations of the LATIN algorithm, this criterion is based on the value of ξ i = || wi,n+1 -ŵi,n+1 || ki /||w i,nŵi,n || ki . If ξ i is less than a threshold (chosen as 0.9 in the following examples), the enrichment step is skipped and the next iteration of the LATIN algorithm is performed, jumping directly to the next coupled stage. For more details on this aspect, one can refer for example to [START_REF] Heyberger | Multiparametric analysis within the proper generalized decomposition framework[END_REF].

Enrichment step.

-Assuming now that a new decomposition ūm-1 i of order (m -1) is known, the reduced-order basis is enriched by seeking an approximation u m i = ūm-

1 i + ∆u = ūm-1 i + Λ i λ i of order m. For that, a new pair (Λ i , λ i ) is added solving: Find (Λ i , λ i ) ∈ V i × I ; ∀(Λ i , λ i ) ∈ V i × I, A i ( ūm-1 i + Λ i λ i , Λ i λ i + Λ i λ i ) = L i (Λ i λ i + Λ i λ i ) ( 46
)
This problem is rewritten in two equations:

Find (Λ i , λ i ) ∈ V i × I ; ∀Λ i ∈ V i , A i (Λ i λ i , Λ i λ i ) = L i (Λ i λ i ) -A i ( ūm-1 i , Λ i λ i ) (47)
∀λ i ∈ I, A i (Λ i λ i , Λ i λ i ) = L i (Λ i λ i ) -A i ( ūm-1 i , Λ i λ i ) (48)
or, in a more compact form:

Find (Λ i , λ i ) ∈ V i × I ; ∀Λ i ∈ V i , A i (Λ i λ i , Λ i λ i ) = Bm i (Λ i λ i ) ( 49
)
∀λ i ∈ I, A i (Λ i λ i , Λ i λ i ) = Bm i (Λ i λ i ) (50)
where one introduces Bm

i (•) = L i (•) -A i ( ūm-1 i , •).
The details of the resolution of nonlinear system (49,50) can be found for example in [START_REF] Néron | Proper Generalized Decomposition for multiscale and multiphysics problems[END_REF] but, roughly, a fixed-point method initialized by λ i (t) = 1 is used. Knowing time function λ i , Equation (49) allows to compute a space function Λ i . Then, knowing space function Λ i , Equation (50) allows to compute a time function λ i . The process can be performed until convergence but, in practice, only 2 iterations are used. Before adding the new pair to the approximation, a Gram-Schmidt algorithm (with respect to the scalar product associated to

a i (•, •)) is performed. Again, a new approximation u i ≈ u m
i being known, the corresponding w i and Φ i can be computed using the restriction operator and the search direction.

Details of the coupled-decoupled strategy

The details of the whole procedure are given in the Algorithm 1.

Algorithm 1: Summary of the procedure, based on a succession of coupled problems to accommodate the models and decoupled problems solved by PGD

1 begin 2 Initialization construction of s 0 ∈ A d 3 for i = 1, 2 do 4
Solve decoupled problem (34) for model i to find (u i , w i , Φ i ), assuming ŵi = 0 and Φi = 0 Both examples concern proportional loadings that are, in spite of their simplicity, probably ones of the most common loadings encountered in structural mechanics analysis. In that case, as this first demonstration concerns linear static elasticity, the solution of the monolithic formulation is also proportional and a PGD approach that would apply directly of the coupled monolithic formulation must build a solution with one pair. However, the decoupled treatment of the models leads to a different behavior. During the initialization of the LATIN algorithm, the two models are solved independently, without taking into account of their coupling, which leads to a first PGD mode for each model. The coupling is then recovered iteratively and, for that purpose, new modes are generated. This behavior is similar to that of a classical domain decomposition approach in which the materials subdomains are glued by the mean of interface conditions. The motivation of this choice of proportional loading in the examples is to separate the behavior of the PGD to represent complex solutions from the effect of using a decoupled algorithm to build the different reduced models.

5 s 0 = {(u i , w i , Φ i )} i=1,
The case of nonlinear problems, such as the coupling between a global linear domain and a local viscoplastic patch, is under investigation. Following our former works on viscoplasticity in the context of the LATIN method, it will only impact the coupled stage of the algorithm, in which the nonlinear evolution law will be also solved. The iterations (and then the modes that will be generated) will allow to satisfy both the nonlinear behavior and the coupling between the models. This aspect will be the subject of a further publication.

Open-hole plate

In this first example, we come back to the illustration of the introduction, which consists of an open-holed plate simulated under the assumption of plane stress. The two geometries that are considered are presented in Figure 7. The overall dimensions are 100 mm × 100 mm. The diameter of the round hole and the side length of the square hole are 20 mm. No body force f is considered, but a sinusoidal-cycle linear force F with an amplitude of 10 MPa is prescribed on the right side. The time interval I = (0, T ) with T = 1 s is discretized using 100 time steps. The material is linear elastic and isotropic with a Young modulus E = 210 GPa and a Poisson coefficient ν = 0.3. Figure 9 illustrates the two possible computational approaches. On the upper part, two different meshes have been generated, to compute the solutions with a monomodel approach. On the lower part, the multimodel simulation is performed using a coarse global model of dimensions 100 mm × 100 mm and two refined patches of dimensions 50 mm × 50 mm corresponding to the geometrical details. Apart from the gluing zone which allows to adapt the models to each other, these patches correspond to the zone of interest of the structure (delimited by a red line in the figure) as the stress will be maximum around the hole.

F Ω 1 Γ F Γ u F Ω 1 Γ F Γ u time (s)
force (MPa) T 0 10 The LATIN algorithm is stopped when the indicator η LAT IN < η = 10 -2 . Figure 8 shows the convergence of the algorithm by plotting the evolution of the indicator along the number of iterations. One can notice that the value of the search direction parameters k 1 = k 2 = k influences only slightly the convergence rate around the value k = 1 which is selected herein. In future works, in particular including some nonlinearities, further investigations will be derived from the numerous works on this subject in the LATIN method. As already mentioned, the behavior of the algorithm is similar to that of a domain decomposition approach and the advanced techniques (multiscale features, preconditionners ...) that are proposed in the literature to accelerate the convergence of this type of algorithms will be investigated carefully and will certainly lead to an improvement of the convergence rate. The four first PGD modes generated by the algorithm for each geometry are presented in Figure 10 whereas Figure 11 shows the upper-left quarter of the plate in order to study the zone of interest around the hole at t = T /4. One can notice a very good agreement of the Mises stress obtained by the monomodel and the multimodel approach.

Gamma-shaped structure

The second example consists of a Gamma-shaped structure with four different geometries in the inside corner, presented in Figure 12. The overall The same procedure than for the previous case is followed and Figure 13 shows that the various local patches can be embedded in the simulation to study the influence of different shapes and, for example, chose a design which satisfies some prescribed conditions such as the limit of elasticity.

Conclusion

In this paper, we investigated the coupling of reduced models for the simulation of structures involving localized geometrical details. For that purpose, we used the Arlequin method to setup a convenient framework to deal with multimodel problems and the LATIN strategy to solve the models independently. This approach allows a very flexible solution of the initial problem. Indeed, different codes can be used to solve the various models and different approaches can be mixed. Herein, the idea was to use the Proper Generalized Decomposition to reduce the computational cost of the models and the promising first The aim of this first work was to assess the feasibility of the technique and its flexibility, that is why no deep study of the computational gains has been performed, which is the purpose of further works in progress. One of the advantages of the approach is that it permits to use different PGD solvers, dedicated to the specificities of the models that are considered (linear/nonlinear, deterministic/stochastic, atomistic/continuum ...). The coupling of elastic and viscoplastic reduced-order models (following [START_REF] Relun | A model reduction technique based on the PGD for elastic-viscoplastic computational analysis[END_REF]) to capture the nonlinear behavior of material in specific zones is currently under investigation.

F Ω 1 Γ F Γ u F Ω 1 Γ F Γ u F Ω 1 Γ F Γ u F Ω 1 Γ F Γ u time (s
In the previous numerical examples, the reduced-order bases corresponding to the two models are built along the iterations of the LATIN algorithm. An interesting approach, that will be investigated in the next future concerns the coupling of different models, some of which being precomputed and stored as PGD modes in a library. This would allow to increase the performances of the strategy and to optimize an industrial structures by exploring the influence of various designs. In that context, some questions will need to be answered such as

Figure 1 :Figure 2 :

 12 Figure 1: Remeshing global structure involving a local geometrical detail

Figure 4 :

 4 Figure 4: Arlequin models and the coupling zone (in grey)

Figure 6 :

 6 Figure 6: The coupled and decoupled stages of the LATIN method at Iteration n + 1

2 6 n ← 0 7 9 11 Decoupled 18 ( 19 s 20 n ← n + 1 4.

 279111819201 while error indicator η LAT IN > threshold do 8 begin Coupled stage solve coupled problem (29) to find {( ŵi , Φi )} i=1,2 ∈ Γ 10 begin stage solve decoupled problems (38) to find s n+1 ∈ A d using a PGD approach in two steps for each model: 12 for i = 1, 2 do 13 Update step: solve (45) for model i, reusing the previous reduced-order basis, to find ( ūi , wi , Φi ) 14 if skipping criterion ξ i > threshold then 15 Enrichment step: solve system (49,50) for model i to add a new pair and find (u i , w i , Φ i ) = ( ūi , wi , Φi ) + (∆u i , ∆w i , ∆Φ i ) 16Orthonormalization of the reduced-order basis of model i17 else u i , w i , Φ i ) = ( ūi , wi , Φi ) n+1 = {(u i , w i , Φ i )} i=1,2 Numerical examplesIn this section, to illustrate the capabilities of the strategy, we consider two numerical examples of coupling a geometrical local detail in a global structure. The first one is an open-holed plate where two shapes of the hole are considered: round or square. The second is a gamma-shaped structure where the inside corner can have four different designs.

Figure 7 :Figure 8 :

 78 Figure 7: Open-holed plate (the zone of interest is delimited in red)

e 12 cFigure 9 :

 129 Figure 9: Meshes of 3-nodes triangles with linear interpolation (top: monomodel simulations with 1,696 and 1,644 DOFs; bottom: multimodel approach with a same global mesh with 505 DOFs and two local meshes with 1,016 and 944 DOFs, witdh of the gluing zone e c 12 = 5 mm) and the zone of interest delimited by a red line

Figure 10 :

 10 Figure 10: Four first PGD modes for local model 2 {Λ k 2 } k=1...4 , top: for the round hole; bottom: for the square hole, (using the symmetry, only the upper part is represented herein.

Figure 11 :

 11 Figure 11: Influence of local geometrical details by embedding local reduced models (Mises stress at t = T /4; top: monomodel approach; bottom multimodel simulation)

Figure 12 :

 12 Figure 12: Gamma-shaped structure with different designs in the inside corner (the zone of interest is delimited in red)

Figure 13 :

 13 Figure 13: Influence of the local geometrical details by embedding local reduced models (Mises stress at t = T /4)

  ) is searched in two steps. The first (called update) consists in seeking a solution ūm-1

	i	of order (m -1) in the reduced-order basis {Λ k i } m-1 k=1
	without any enrichment but only by computating new time functions. The
	second (called enrichment), performed only if the update step is not sufficient,
	consists in adding a new pair to the decomposition to compute a new solution
	u m

i of order m. These steps are described hereafter.

Update step.

-Assuming that a decomposition u m-1 i

Acknowledgments

The authors want to thank the French National Centre for Scientific Research (CNRS) and the French National Research Agency (projet number ANR-14-CE07-0007 CouESt) for the fundings provided to this work.

The numerical simulations were performed using the routines CArl, freely available at https://github.com/cottereau/CArl.