Karol Marty

Pierre Bordet
email: pierre.bordet@grenoble.cnrs.fr

Virginie Simonet

Mickael Loire

Rafik Ballou

Céline Darie

Jakob Kljun

Pierre Bonville

Olivier Isnard

Pascal Lejay

Sr, Ca, B=Ta, Nb, Sb A=ba

Si K D=ge

Marty

B Zawilski

C Simon

Magnetic and dielectric properties in the langasite

Keywords: numbers: 75.25.+z, 71.27.+a magnetic frustration, helical magnetic order, langasite, dielectric properties

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

Multifunctional materials, in which coexist different physical properties (optical activity, laser properties, piezo-and/or ferroelectricity, ferroelasticity, magnetic ordering, etc), are the subject of an intense research effort owing to their strong potential for various industrial applications. To date, the basic mechanisms at the origin of the coupling between these properties inside a single compound are still poorly understood. Experimental realizations of such materials are rare, and this is particularly the case in the field of multiferroics, i.e. compounds which display both ferroelectric and magnetic ordering 1 . Proper multiferroics are defined as ferroelectric compounds which undergo a magnetic ordering at lower temperature, as for example BiFeO 3 2 . In these compounds, there is a large temperature gap between the ferroelectric and magnetic transition, and the coupling between the two properties is generally weak. For improper multiferroics (for example TbMn 2 O 5

3 ), the electric polarization is induced by the onset of a complex magnetic order within a paraelectric and generally centrosymmetric phase. Such materials therefore present a naturally strong coupling between magnetization and electric polarizations thus giving the opportunity to study the microscopic mechanisms driving such couplings.

The langasite family, the prototype of which is La 3 Ga 5 SiO 14 , has been widely studied for the striking piezoelectric and non-linear optical properties of its members [START_REF] Mill | Proc. 2000 IEEE Int. Freq. Contr. Symp[END_REF][START_REF] Bohm | [END_REF]6,7 , related to the non-centrosymmetric nature of its Ca 3 Ga 2 Ge 4 O 14 structure type 8 . These compounds have generated a strong interest for applications in bulk acoustic waves (BAW) and surface acoustic waves (SAW) devices, as well as in non-linear optics and electrooptics [START_REF] Mill | Proc. 2000 IEEE Int. Freq. Contr. Symp[END_REF]9,10 .

The langasite structure belongs to the trigonal non-centrosymmetric P 321 space group. The general formula is A 3 BC 3 D 2 O 14 , thus containing four different cationic sites. The decahedral A site and the octahedral B site form a layer at z=0, while the two tetrahedral sites C and D are located on the plane z=1/2 (see Figure 1) [START_REF] Mill | Proc. 2000 IEEE Int. Freq. Contr. Symp[END_REF]15 .

Due to its atomic arrangement, this structure type is able to accommodate a large number of different cations with various sizes and valences, leading to a wide variety of isostructural compounds [START_REF] Mill | Proc. 2000 IEEE Int. Freq. Contr. Symp[END_REF]11,13 . Among these, some contain magnetic cation sublattices. Therefore, this family might provide interesting examples of coexisting magnetic order and optical or electrical properties. Recently, we have undertaken the general investigation of the magnetic langasites 14 . In this article, we present the study of the magnetic and dielectric properties of several Fe-containing langasite compounds with formula A 3 BFe 3 D 2 O 14 , with A=Ba, Sr, Ca, B=Ta, Nb, Sb and D=Ge, Si.

II. EXPERIMENTAL A. Synthesis and structural characterization

Powder samples were synthesized by solid state reactions. The initial components of the powder mixtures were high purity oxides or carbonates in stoichiometric quantity. The mixtures were ground, compacted and heated up to 1150 -1250˚C in Al 2 O 3 crucibles in air for 7 hours. The pellets were heated at a rate ranging from 200˚C/h to 300˚C/h. The samples containing GeO 2 were kept at 1000˚C for 2h before reaching the final temperature to prevent the loss of the relatively volatile germanium oxide. The synthesis conditions were chosen taking into account some properties of the initial and final compounds, such as the decomposition of the carbonates, the oxidation of Sb 2 O 3 , the fusion temperatures of GeO 2 , SrO and of the final product. The products were analyzed with x-ray diffraction using a D8 Brucker diffractometer in transmission geometry with the CuKα1 radiation, to check the presence of possible impurities. As some of the products showed minor impurities (max. 3 %), they were ground, compacted again and annealed for 10-50h, up to three times in order to decrease the impurity content.

Rietveld refinements with the program Fullprof 17 were carried out on the final products diffractogramms. In all cases, the expected stoichiometry and the full occupation 18 . This technique has the advantage of being fast and exempt of contamination from a crucible. The single crystals were cut out of the grown ingot with a size roughly of 3 mm in diameter and 8 mm in length for the neutron diffraction measurements. The crystallographic quality were checked by Laüe photographs. No impurity phase could be detected by x-ray powder diffraction of crushed crystal pieces. A small crystal piece of each compound was used for single crystal x-ray diffraction analysis using a Kappa ApexII Bruker diffractometer equipped with graphite monochromatized AgKα radiation. This confirmed the crystal quality and the langasite type structure and allowed to determine the handedness of each of these chiral crystals by use of anomalous scattering sensitivity. In all the compounds under study, the magnetic Fe 3+ cations form a peculiar array of triangles on a triangular lattice, in the plane at z=1/2 (figure 2). Two consecutive planes are separated from each others by a layer of A 2+ and B 5+ cations. The Fe 3+ triangles are equilateral by symmetry. The distance between the Fe 3+ cations is 3.662(8) Å for BNFS. Two equivalent super-exchange interactions J1 between two nearest neighbor magnetic cations on one triangle are each mediated via an O3 anion. The Fe-O3-Fe angle 100.1(9)˚and the Fe-O3 distances are 1.73(2) and 2.93(2) Å for BNFS. One Fe 3+ cation of a triangle interacts with four neighbors in the nearest triangles within the plane by the super-superexchange paths J2 through two O2 anions forming an edge of a DO4 tetrahedron (figures 1 and 2). The Fe-O2-O2 angles are 144(1)˚and 140(1)˚, and the Fe-O2 and O2-O2 distances are respectively 1.93(2) and 2.59(4) Å in BNFS. The shortest distance between two Fe 3+ cations from neighbouring planes is the c parameter (5.2353(1) Å), with also two O2 oxygens in the exchange path forming an edge of the BO 6 octahedron. But one needs to distinguish three different inter-plane interactions, each one with multiple possible super-super-exchange paths. J4 is the interaction between a Fe 3+ cation and the cation just above or below in the adjacent plane, while J3 and J5 are the diagonal interactions between this Fe 3+ cation and the two other Fe 3+ cations of the neighbouring triangle (figure 2). It is worth noticing that due to the structural chirality, J3 and J5 are obtained through different exchange paths and expected not to be equal.

B. Magnetization measurements

The magnetostatic properties of BNFS and BTFS single crystals and of all the powder title compounds were investigated under magnetic fields up to 5 T in the temperature range from 2 K to 300 K with a commercial Quantum Design MPMS SQUID magnetometer, and up to 10 T from 1.6 K to 300 K with a purpose-built magnetometer using the anti-Helmholtz two-coil axial extraction method. Figure 3 illustrates for BNFS the thermal evolution of the magnetization divided by the magnetic field (i.e. equal to the linear susceptibility χ for temperatures greater than ≈ 26 K), and of its inverse, measured under a magnetic field of 5 T applied both perpendicular to the c-axis and parallel to it. A cusp in the M/H measurements suggested the onset of a magnetic transition at 26 K, further corroborated by a sharp peak in the specific heat 19 . Above this temperature, no significant anisotropy is observed within the accuracy of the measurements, which is expected for Fe 3+ ions without orbital contribution. Below the magnetic transition, whereas M/H is almost constant for H parallel to the c axis, the in-plane susceptibility drops, although not to zero as would be the case for a canonical antiferromagnetic transition. This implies a more complex magnetic arrangement. This is also illustrated on the magneti-zation curves M (H) shown on figure 3. At high temperature (275 K), the magnetization is the same for the two orientations of the applied field. On the contrary, at low temperature (2 K), an anisotropy of the magnetostatic properties comes out: the magnetization measured with an applied field parallel to the c axis is no longer linear and shows a slight curvature at small fields. This curvature would indicate the rise of a small component along the c axis. At high temperatures (150 K to 300 K), a linear fit of the inverse of M/H yields, in a Curie-Weiss model χ = C/(T -θ), an effective moment µ eff = g S(S + 1)=5.95 µ B , compatible with the value of the Fe 3+ free ion (3d 5 , J=5/2, µ eff =5.92 µ B , M sat =5 µ B ), and a negative Curie-Weiss temperature of -171 K, characteristic of antiferromagnetic interactions.

The magnetic behavior of a BTFS crystal was found essentially similar to that of BNFS, and this was also the case for all the investigated compounds prepared as powders : all of them undergo an antiferromagnetic ordering transition in the 24-35 K range. As an example, the magnetic susceptibility for the Ba 3 SbFe 3 Si 2 O 14 powder compound is shown in Figure 4 where the magnetic transition is clearly visible at T N =35 K. The values of T N , θ and Fe 3+ magnetic effective moments obtained in a Curie-Weiss fit of the inverse magnetic susceptibilities are summarized in Table III. All compounds, whatever the cations on sites A, B and D, present a Néel temperature around 25 K, except those containing Sb 5+ cations on site B. These latter display a magnetic transition around 35 K, and at the same time a noticeably smaller absolute value of the Curie-Weiss temperature. In order to understand the microscopic origin of this difference, and more generally to further our understanding of the magnetic order in these compounds, microscopic probes such as neutron scattering and Mössbauer spectroscopy were used.

C. Neutron diffraction study

Neutron powder diffraction (NPD) measurements were performed at the Institut Laüe Langevin on samples of the compounds Ba 3 NbFe 3 Si 2 O 14 (BNFS), Ba 3 SbFe 3 Si 2 O 14 (BSFS), Sr 3 NbFe 3 Si 2 O 14 (SNFS) and Sr 3 SbFe 3 Si 2 O 14 (SSFS) using the D1B instrument (λ=2.52 Å). These four compounds were chosen in order to track, via the structural and magnetic arrangements determination, the influence of large cell parameter variation (Sr versus Ba) and the microscopic origin of the different transition temperatures (Nb versus Sb). A BNFS single crystal was also studied on the D15 diffractometer of the Institut Laüe Langevin (λ=1.1743 Å). The thermal evolution of the diffraction patterns were measured on D1B using an orange cryostat by increasing the temperature at a rate of 0.42 K/min with 75 scans of 10 minutes each. In order to obtain sufficient statistics at 30 K (above the magnetic transition), 20 K (below the magnetic transition) 1.8 K (for BSFS, SNFS and SSFS) and 1.5 K (the lowest reachable temperature, for BNFS), 6 scans of 20 minutes were summed at each temperature. The magnetic structure was solved from refinement of the D15 and D1B diffractograms measured on BNFS 19,20 . The magnetic structures of the three other compounds were refined using the same ordering scheme (see figure 5). The magnetic peaks could be indexed using an incommensurate propagation vector (0, 0, τ ), with τ ranging from 0.14 for SNFS to 0.19 for BSFS, and the magnetic structure was found to be a helical spin arrangement propagating along the c axis from equal moments lying in the (a,b) plane at 120˚from each other within each triangle (see figure 6). This 120˚arrangement results from the usual compromise of frustrated Heisenberg spins on a triangle-based lattice. All spins rotate around c by ±2πτ between to consecutive Fe 3+ layers. It is worth noting that a helical magnetic order often results either from the presence of frustration of interactions in the helix direction [START_REF] Yoshimori | [END_REF] , or from Dzyaloshinskii-Moriya interactions in a ferromagnetic structure 29 . In the present case, it was shown that the helical stucture is generated by the twist in the magnetic super-super-exchange paths from plane to plane. Furthermore, the in-plane 120˚magnetic arrangement and the helical propagation both possess a definite sense of rotation defined as chirality, which is related to the structural chirality. A thorough investigation of this interesting aspect has been previously reported in ref. 20 and will not be discussed further here.

The main results of the NPD refinements carried out for the four samples at 1.8 K are summarized in Table IV. The magnetic moment values are close to 4 µ B and they are smaller for Nb compounds than for Sb ones. This is markedly smaller than the expected value for a saturated Fe 3+ magnetic moment (5 µ B ). An explanation for this behavior could be the strongly covalent nature of the Fe-O bonds for Fe 3+ cations in tetrahedral coordination, which is known 23 to favor spin transfer towards the oxygen orbitals and therefore decrease the apparent magnetic moment on the Fe 3+ cations. Indeed, in order to obtain accurate refinements of the NPD data, we had to introduce an ad-hoc "magnetic Debye-Waller factor", similar to an overall isotropic atomic displacement term. The refined values of the corresponding B mag parameters are given in Table IV. The effect of this term can be viewed qualitatively as a modification of the magnetic form factor with increased spin density at larger distances from the center of the electron cloud, as would be expected in the presence of spin transfer to the oxygen orbitals. Another qualitative information is given by the diffractograms measured at temperatures up to 20 K above the magnetic transition, where some diffuse scattering centered on the first nuclear peak is visible. This signal is characteristic of the onset of short range mag-netic correlations in the (a, b) planes. The width of the signal decreases, i.e. the correlation length of this shortrange order increases, down to the magnetic transition where it disappears at the expense of the magnetic Bragg peaks characteristic of the 3-dimensional long-range order (see figure 7).

As seen in Table IV, the τ values of the antimony compounds are higher (between 0.177 and 0.196) than those of the niobium compounds (around 0.14). The τ values appear also slightly larger for the barium compounds than for the strontium ones, especially in the Sb compounds. However, the τ values are much less sensitive to the Sr/Ba substitution than to the Sb/Nb one, although the former have a much larger effect on the lattice parameter due to the big difference of their ionic radius (1.56 Å for Ba versus 1.4 Å for Sr). After examination of the refined structures of the four compounds, and especially focusing on the lengths and bonds angle of the magnetic exchange paths, no systematic trends could be extracted from the data that would explain the peculiarity of the Sb based compounds magnetic properties on pure structural grounds. The influence of the Sb cation in the Felangasites is therefore more probably ascribable to the ion electronic state than to induced structural distortion. We will come back to this aspect in the discussion.

D. Mössbauer measurements

Mössbauer absorption spectra on the isotope 57 Fe (E 0 = 14.4 keV) were recorded for the BNFS compound at room temperature and in the temperature range 4.2 to 30 K, in zero external magnetic field (see figure 8). At 27 K and above, the spectrum consists of a quadrupolar hyperfine doublet characteristic of the paramagnetic phase. At 27 K, the isomer shift with respect to the α-Fe reference is δ = 0.35(1) mm/s (1 mm/s corresponds to 11 MHz) and the quadrupolar hyperfine splitting is ∆E Q = |eQV ZZ | /2 = 1.29(1) mm/s, where Q is the quadrupole moment of the nuclear excited state and V ZZ is the principal compo-FIG. 8:

57 Fe Mössbauer absorption spectra in Ba3NbFe3Si2O14 at selected temperatures. The spectra at 4.2, 23 and 26 K (circles) are fitted (black lines) to the helical magnetic structure model determined by neutron scattering. At 26 K, a quadrupolar subspectrum contribution (represented in red) overlapping the magnetic signal suggests that the magnetic transition is first order. nent of the electric field gradient (EFG) tensor at the 57 Fe nucleus site. These values are typical of a Fe 3+ ion, although the ∆E Q value is somewhat larger than generally found in ferric insulators, indicating a strongly distorted Fe site, compared to a perfect tetrahedral environment, in agreement with crystal analysis. The point symmetry at the Fe site is low (a two-fold axis along the crystal a axis), and thus the directions of the EFG tensor axes cannot be determined from symmetry considerations alone.

Below 27 K, a magnetic hyperfine structure appears, with 6-line characteristic spectra showing the presence of a magnetic hyperfine field acting on the nuclear moment. This is the fingerprint of a long range magnetic order. The 4.2 K spectrum (see figure 8) corresponds to a single saturated hyperfine field H hf = 440(2) kOe; considering the hyperfine proportionality constant typical for Fe 3+ insulators of 110 kOe/µ B , this yields a unique Fe 3+ magnetic moment value of 4 µ B , in agreement with the neutron derived value. Close inspection of the 4.2 K spectrum reveals slight inhomogeneous line broadenings, most clearly seen when comparing the widths of the two central lines (labelled 3 and 4 in figure 8) and of the two intermediate lines (2 and 5). This spectral effect is due to a distribution of one or more hyperfine parameters. The magnetic structure determined by neutron diffraction being an incommensurate transverse spiral with propaga-tion vector along c, then the angle between the hyperfine field, proportional to the Fe 3+ moment, and the principal axis OZ of the EFG tensor, linked to the crystal axes, shows a distribution. The second order shifts of the hyperfine energies due to the quadrupolar coupling then lead to a broadening which is different for each line 12 . The line fits in the magnetic phase in Fig. 8 were performed using this model, the key parameter determining the size of the broadenings being the angle Θ between OZ and the propagation vector, which is found to be Θ ≃ 36 • . On heating, the lines become broader and the overall magnetic splitting decreases, indicating that the mean hyperfine field decreases and that a distribution in hyperfine field values appears. The broadenings are significant near T N (see the spectra at 23 and 26 K in figure 8). At 26 K, a small intensity paramagnetic doublet (13% relative weight) is present together with the magnetic order spectrum. This coexistence of paramagnetic and magnetically ordered regions near T N is usually associated with a certain degree of first order for the magnetic transition. At 27 K, the sample as a whole is paramagnetic.

E. High resolution X-ray diffraction

The remarkable magnetic properties of the Felangasites underlined in the previous sections may also have consequences on their dielectric behaviour. Indeed, the non-centrosymetric P 321 space group symmetry of the langasites does not allow the existence of an electric polarization. However, we have established that the Fe-langasite compounds undergo a transition towards a magnetic helical order whose symmetries are compatible with ferroelectricity: preservation of the 3-fold axis and loss of the 2-fold axes. It is therefore interesting to determine whether a structural phase transition and a lowering of symmetry authorizing a ferroelectric order appear in the temperature region of the magnetic transition. The D1B data have not revealed any such modification in the four samples studied. However the limited Q-range and spatial resolution make this instrument moderately suitable for the detection of subtle structural anomalies.

In order to follow the structure as function of temperature with a much higher sensitivity, we have used high resolution synchrotron powder diffraction at the ID31 beam line of the ESRF, equipped with a Ge(111) multianalyzer stage. This experiment was performed on the Ba 3 SbFe 3 Si 2 O 14 compound. The sample contained in a 1 mm diameter quartz capillary was placed inside a helium flow horizontal cryostat allowing sample rotation. Data collections were carried out in the range 2θ=5˚to 35˚at 26 temperatures between 250 K and 10 K at an incident energy of 30 keV (λ=0.39816(2) Å). The diffractograms were analyzed by Rietveld refinement using the Fullprof program. The background was interpolated and the Bragg reflection shapes were described using a Thomson-Cox-Hastings profile function and a uniaxial description of anisotropic peak broadening about the [001] direction. All positional parameters in the P 321 space group description were refined, and all atoms were given isotropic atomic displacement parameters (a.d.p.), those of all oxygen atoms being constrained to be equal. The Rietveld plot for the 250 K data refinement is shown in figure 9. The refined parameters at 250 K and 10 K are given in table V and do not exhibit any anomalous difference, the only marked variations being the expected cell contraction and the decrease of the a.d.p.'s with temperature. Figure 10 a, b, c andd display the temperature evolution of the cell parameters, quadratic term U of the Gaussian peak width function, and a.d.p.'s of the cations. An obvious anomaly is detected in the a-parameter variation below 35 K (i.e. the magnetic ordering transition temperature for this compound). However the amplitude of the anomaly is quite small since the cell parameter change between 40 K and 10 K is only about 0.0015 Å. The U parameter depends on the instrumental resolution (which does not vary with temperature) and the presence of strain in the sample. This parameter, related to Gaussian strain, increases on cooling and then starts decreasing below about 40 K. This could be related to the building up of strain in the compound followed by a relaxation at lower temperature due to a structural reorganisation associated to a phase transition. In this context, the increase of the cation isotropic a.d.p.'s could be due to the accommodation of atomic displacements associated to this structural transition.

If a structural phase transition coupled to the onset of magnetic order was present, then the symmetry of the low temperature phase should keep the three-fold axis and loose the two-fold axis. The expected symmetry would then be P 3, which allows an electric polarization. In this space group, the silicon and each of the oxygen sites are split into two positions. Refinements of the 10 K synchrotron data with the P 3 symmetry were attempted, but did not lead to conclusive results. However, the departure from the P 321 symmetry may be very small leading to hardly detectable atomic displacements. As a conclusion, the results of this temperature dependent high resolution structural study might indicate the presence of a structural phase transition coupled to the onset of magnetic order. Single crystal x-ray diffraction experiments are underway to determine the nature of the low temperature phase.

F. Dielectric measurements

In order to test the possible influence of the magnetic ordering on the dielectric properties of Fe-langasites, dielectric measurements were performed at the Néel Institute by measuring the complex impedance of our samples using a commercial HP 4284A LCR-meter. The purposebuilt experimental setup consists of a sample holder protected by a dewar for a direct immersion in liquid helium of an helium container. Four coaxial cables linked to the electronic devices allow the complex impedance measurement. A model of capacitor and resistance in parallel was selected to extract the dielectric permittivity constant proportionnal to the capacitance of the sample, designed to approach the ideal geometry of a parallel plate capacitor. The metallic electrodes are realized by sputtering with a silver target. Systematic check of the amplitude and frequency dependence were undertaken to establish the better conditions of measurements ruling out extrinsic effects. The chosen conditions were finally an amplitude voltage of 1V and a frequency within the range of 10-100 kHz. The measurements under magnetic field were performed at the CRISMAT laboratory, with a similar experimental setup designed to work in a commercial Quantum Design PPMS and allowing to reach magnetic fields of 9 T.

The powder samples of BNFS, BTFS and BSFS were prepared by pressing thin pellets which were annealed at 1100˚C during 78 hours. Single crystal samples of BNFS and BTFS were oriented and cut within the rods obtained in the image furnace in order to have two different orientations for the applied electric field: along the c axis and along the a* axis.

The thermal variation of the dielectric permittivity constant of all the powder samples shows a small cusp near the temperature associated with the magnetic transition (24 K for BNFS and BTFS, 30 K for BSFS). This cusp is a contribution from the dielectric permittivity within the (a,b) plane. It is indeed clearly visible in the single crystal samples with the electric field applied along the a* axis (also reported in ref. 32 ), while it is absent for the other orientation perpendicular to the (a,b) plane which only shows a broad variation at higher temperature without any slope discontinuity (see figure 11 for BNFS). For a powder sample (see BSFS measurement on figure 12), the observed cusp is obviously smaller than for a single crystal, as crystallites of every orientation, even those not presenting this effect, contribute to the signal.

The variation of the relative increase of the dielectric permittivity as a function of the magnetic field applied in the direction of the electric field, perpendicular to the c axis, is illustrated for different temperatures by the figure 13. It shows mainly a quadratic dependence with the field. As the dielectric permittivity is the second partial derivative of the free energy with respect to the electric field, the magneto-electric coupling terms measured here correspond to the terms of second order in electric field within the free energy development:

- 1 2 ǫ 0 ǫ ij E i E j - 1 2 γ ijk E i E j H k - 1 2 η ijkl E i E j H k H l -...
The results obtained for the Fe-langasite means that only η ijkl contributes to the measured magneto-electric effect. It is worth noting that this coupling term is always present in the free energy development, as it is of second order in both electric and magnetic field, so no particular symmetry conditions determine its existence, unlike for γ ijk . This quadratic coefficient however varies with the temperature, being maximum close to the Néel temperature. This variation reflects the critical behavior of the magneto-electric effect describable in the frame of the fluctuation-dissipation theory 31 . Other magneto-electric coupling terms are accessible when measuring the magnetic dependence of the electric polarisation (which is the first partial derivative of the free energy with respect to the electric field). Such experiments are underway.

III. DISCUSSION

A. Magnetic properties

Our investigations of the Fe-langasites magnetic properties have revealed an original magnetic order at low temperature associated to the triangular array of triangles of magnetic Fe 3+ . The magnetic moments are orientated at 120˚from each other within the triangles and propagate helically perpendicular to the triangular planes. There are at least two energy scales in this system since the intra-triangle super-exchange is expected to be much stronger than the in-plane and out-of-plane inter-triangle interactions which are mediated by two oxygens. The Curie-Weiss temperature mainly reflects this intra-triangle interaction which also gives rise to magnetic fluctuations above the ordering temperature (see figure 7). The Curie-Weiss temperature is larger than the Néel temperature by almost an order of magnitude. The 3 dimensional ordering sets in when the temperature becomes of the order of the super-super exchange paths within the planes and between them. We have previously shown 20 that this helical arrangement results from a twist of the interplane interaction: the strongest interaction among the three that link one Fe 3+ of a triangle to the three Fe 3+ of the triangle of the next plane is diagonal (J3 or J5). If, moreover, the two other interactions are zero, it is easily shown that the z-component of the propagation vector would be exactly τ =1/6=0.167, i.e. a rotation of the moments of 60˚from plane to plane 20 . The departure of the propagation vector from 1/6 arises when the two other inter-plane interactions are no more negligible.

In the present study, we have focused on the effect on the magnetic properties of different non-magnetic cations on sites A, B and D. The most striking result is that all the magnetic characteristics are roughly identical for all compounds even when the substitution produces a large variation of the cell parameters, except for those containing Sb (in place of Nb or Ta). The presence of this cation systematically increases the Néel temperature, decreases the Curie-Weiss one, increases the value of the ordered magnetic moment, and increases the τ value (i. e. shortens the helical modulation period). The origin of this effect could be structural, the Sb 5+ radius being slightly smaller than the Ta 5+ or Nb 5+ ones. However, no systematic trend on exchange path geometry was revealed by the structural investigation performed with powder neutron scattering on chosen compounds. The microscopic mechanism at play should therefore be of electronic origin. It could be related to the larger Pauling electronegativity of Sb 5+ (2.05) with respect to Nb 5+ (1.6) and Ta 5+ (1.5) 30 . The departure of this value from the oxygen electronegativity (3.44) gives the degree of ionicity/covalency and polarisation of the cation-oxygen bond involved in the magnetic exchange paths. This should therefore give a more covalent character to the Sb-O bond. The electronic structure of the Sb 5+ cation is also different from the Nb/Ta ones (different column of the periodic table) with in particular a full outermost 4d shell.

Remembering that two consecutive Fe 3+ triangles along the c axis are coupled via oxygen anions forming the octahedral coordination of the D site cation (see figure 2), it is not surprising that a change of the electronic affinities of this cation has consequences onto the magnetic coupling within and between the planes. For instance, the increased covalence of the Sb-O bond seems in turn to decrease the strength of the J1 Fe-O-Fe bond involving the same oxygen ions. The consequence of this is the observed lowering of the Curie-Weiss temperature and smaller spin transfer from the Fe 3+ to the oxygen (hence a larger ordered moment).

The influence of the electronic properties of Sb on the Néel temperature and on the helix propagation vector is more difficult to track since it implies more complicated exchange paths between two planes. It is however interesting to note that the angle of rotation of the spin is larger than 60˚for Sb compounds and smaller than 60˚for the Nb/Ta ones. This implies a contribution of the two weaker interplane interactions, whose ratio would be inverted between Sb compounds and Nb/Ta ones.

B. Dielectric properties

It is of course of interest to address the question of the influence of this increased Sb-O bond covalence on the dielectric properties of the Fe-langasites. The dielectric measurements on all powder samples have not shown significant difference though, with a similar cusp of the dielectric permittivity close to T N , suggesting the same kind of magneto-electric coupling for the different compounds. This comparative study will be carried on with electric polarisation measurements.

The dielectric properties reported here for the Felangasites are strikingly similar to those of YMnO 3 . The magnetic structure of this compound presents some resemblance with the Fe-langasite one : a triangular array of in-plane magnetic moments oriented at 120˚from each other in the ordered phase. The dielectric measurements performed on this multiferroic compound show a comparable anomaly in the thermal variation of the (a,b) in-plane dielectric permittivity while nothing is visible for the dielectric permittivity measured along the c axis, and a quadratic dependence of the dielectric permittivity close to this anomaly with the magnetic field (note however that in the YMnO 3 study, the magnetic field was oriented perpendicular to the electric one). This anisotropic dielectric behaviour of YMnO 3 has been described in a phenomenological Landau description of ferroelectric antiferromagnets including a fluctuative part 31 . Unlike the langasites, YMnO 3 is already ferroelectric at temperatures well above the magnetic transition (proper multiferroic). Anyway, the magnetic ordering still induces an additional contribution to the overall electric polarisation. The rise of electric polarisation at the Néel temperature has also been claimed very recently in the Fe-langasite by Zhou et. al. 32 from polarization measurements performed along the c-axis. This would imply that, at variance with YMnO 3 , these langasite compounds should be classified as improper multiferroics, the magnetic transition driving a paraelectric-ferroelectric transition. Such behaviour has already been reported in another compound containing an Fe triangular lattice where the electric polarisation is induced at the Néel transition in the c direction along the axis of the helix and perpendicular to the plane of the moments oriented at 120˚from each other 33 .

IV. CONCLUSION

The Fe Langasites are the first compounds of this family to evidence a magnetic ordering whose peculiar ar-rangement reflects magnetic frustration and the chirality of the structure. This rich family offers several Fe-based compounds with different cationic substitutions among which we have highlighted the singular role of the Sb cation on the magnetic properties. The langasite structure is non-centrosymmetric and a further loss of symmetry elements is achieved in the magnetic structure that leave an overall polar structure authorising the appearance of a spontaneous electric polarisation. This would then be associated to a structural phase transition as suggested by high resolution x-ray diffraction. Dielectric measurements reveal the presence of magneto-electric coupling through an anomaly at the Néel temperature in the thermal variation of the dielectric permittivity and its quadratic dependance with the magnetic field. The Fe-langasite could then be a new example enlarging the class of magnetoelectric/multiferroic triangle-based antiferromagnets.

FIG. 1 :

 1 FIG. 1: Polyhedral representation of the langasite structure with its four different cation sites on two layers.

  I: Refined structural parameters obtained from xray powder diffraction on BNFS with a D8 Brucker diffractometer. a. d. p. is the atomic displacement parameter. The agreement factor of the fit is given by RBragg = 3.2 ; Rwp = 27.0 ; χ 2 = 1

FIG. 2 :

 2 FIG. 2: Magnetic network of Fe 3+ cations projected along the c-axis (left), with the different intra-plane (left) or inter-plane (right) exchange paths (dashed lines).

FIG. 3 :FIG. 4 :

 34 FIG.3: Magnetization measurements on a Ba3NbFe3Si2O14 single crystal under a magnetic field applied parallel and perpendicular to the c-axis. Top: Magnetic isotherms at 2 K and 275 K. Bottom: Thermal variation of the M/H and its inverse under an applied magnetic field of 5 T. The continuous line represents the linear fit of H/M in a Curie-Weiss model.

14 FIG. 5 :

 145 FIG.5: Rietveld plot for Ba3NbFe3Si2O14 (top left), Ba3SbFe3Si2O14 (top right), Sr3NbFe3Si2O14 (bottom left) and Sr3SbFe3Si2O14 (bottom right) neutron powder data at 1.8 K (1.5 K for BNFS) from D1B, including nuclear and magnetic phases.

FIG. 6 :

 6 FIG.6: Magnetic structure of the Fe-langasites: within each triangle of the (a,b) plane, the Fe 3+ magnetic moments are rotated by 2π/3 from each other. The global phase factor, which cannot be determined by neutron diffraction, was arbitrary chosen equal to zero in this representation (one moment pointing towards the center of the triangles). Along the c-axis, the helix propagates with a period of 1/τ in cell unit.

FIG. 7 :

 7 FIG.7: Powder neutron diffraction patterns of the compound Ba3NbFe3Si2O14, below the magnetic transition, at 1.5 K, and above the magnetic transition, at 30 K.

FIG. 9 :

 9 FIG. 9: Rietveld plot for Ba3SbFe3Si2O14 x-ray powder data at 250 K from ID31; RBragg = 2.58, Rwp = 23.4, χ 2 = 3.6.

FIG. 10 :

 10 FIG.10: Thermal variation of some refined parameters from ID31 x-ray powder diffraction data on Ba3SbFe3Si2O14: a) cell parameter a ( Å). b) cell parameter c ( Å). c) quadratic term U of the Gaussian peak width function. d) a.d.p.'s of Ba 2+ , Sb 5+ and Fe 3+ cations.

FIG. 11 :

 11 FIG.11: Thermal variation of the electric capacity (proportional to the dielectric permittivity) of BNFS single crystal for an electric field applied perpendicular to the c axis (left) and parallel to it (right).

FIG. 12 :FIG. 13 :

 1213 FIG. 12: Thermal variation of the electric capacity (proportional to the dielectric permittivity) of a pellet of BSFS powder. The red line is a linear extrapolation of the data above 30 K, shown to emphasize the cusp.

TABLE

TABLE II :

 II Structural parameters and agreement factors of Fe-langasites as refined from x-ray powder diffraction data at room temperature with a D8 Brucker diffractometer.

	Langasite	a ( Å)	c ( Å)	RBragg Rwp χ 2
	Ba3SbFe3Ge2O14 8.6174(1) 5.2708(1)	3.9 22.7 3.2
	Ba3TaFe3Ge2O14 8.6166(1) 5.2618(1)	6.0 22.5 3.9
	Ba3NbFe3Ge2O14 8.6073(2) 5.2686(2)	6.1 22.9 2.5
	Ba3SbFe3Si2O14 8.5085(1) 5.25024(8) 3.1 20.2 2.8
	Ba3TaFe3Si2O14 8.5231(2) 5.2354(2) 13.5 33.6 2.2
	Ba3NbFe3Si2O14 8.5227(2) 5.2353(1)	3.2 27.0 1.5
	Sr3SbFe3Si2O14 8.2888(2) 5.1445(2)	5.8 30.6 2.1
	Sr3TaFe3Si2O14 8.2777(2) 5.1283(1)	4.3 27.0 2.0
	Sr3NbFe3Si2O14 8.2563(2) 5.1306(1)	3.2 27.6 1.6
	Ca3SbFe3Si2O14 8.1118(2) 5.0570(2)	5.7 26.0 2.0

of all sites was found, without any sign of intersite substitution. As an example, the results for Ba 3 NbFe 3 Si 2 O 14 (hereafter referred to as BNFS) are reported in Table

I

. Table II summarizes the principal structural parameters for the investigated samples. Single crystals of Ba 3 NbFe 3 Si 2 O 14 and Ba 3 TaFe 3 Si 2 O 14 (BTFS) were grown by the floatingzone method in a Cyberstar image furnace, under a 99% Ar + 1% O 2 atmosphere, at a growth rate of 10 mm/h

TABLE III :

 III Magnetic parameters obtained from magnetization measurements of Fe-langasites: Néel temperatures and Curie-Weiss temperatures.

	Langasite	TN (K) θ (K)
	Ba3NbFe3Si2O14	26	-173
	Ba3NbFe3Ge2O14	24	-177
	Ba3TaFe3Si2O14	28	-160
	Ba3TaFe3Ge2O14	24	-165
	Ba3SbFe3Si2O14	35	-114
	Ba3SbFe3Ge2O14	34	-111
	Sr3NbFe3Si2O14	26	-242
	Sr3SbFe3Si2O14	36	-128

TABLE IV :

 IV Refined lattice parameters, magnetic moment value, overall isotropic atomic displacement parameter, zcomponent of the propagation vector and goodness of the fit parameters obtained from neutron powder diffraction data of Fe-langasites at 1.8 K (1.5 K for BNFS) on D1B.

	Langasite a ( Å) m (µB) Bmag	τ	RBragg χ 2
		c ( Å)		Rmag
	BNFS 8.5026(6) 4.04(5) 3.9(9) 0.1429(2) 0.7 28.0
		5.2136(7)		5.8
	SNFS	8.297(1) 3.87(9) 5(2) 0.1398(3) 0.9 35.5
		5.155(1)		9.2
	BSFS 8.5048(5) 4.26(3) 2.1(4) 0.1957(1) 1.2 44.3
		5.2449(4)		4.8
	SSFS 8.2896(9) 4.31(5) 4.6(7) 0.1769(2) 1.8 20.3
		5.140(1)		14.2

TABLE V :

 V Refined structural parameters obtained from high resolution x-ray powder diffraction on BSFS on ID31, at 250 K (RBragg = 2.58 ; Rwp = 23.4 ; χ 2 = 3.6) and 10 K (RBragg = 2.35 ; Rwp = 22.0 ; χ 2 = 0.95). a.d.p. and U are respectively the atomic displacement parameter and the quadratic term of the Gaussian peak width.

		Structural parameters refined at
			250 K
			10 K
	Atom Wyckoff	x	y	z	a.d.p.
	Ba	3e	0.43747(7)	0	0	1.29(1)
			0.43686(6)		1.000(8)
	Sb	1a	0	0	0	1.07(2)
						0.92(2)
	Fe	3f	0.7527(2)	0	1/2	1.27(3)
			0.7523(1)			1.11(2)
	Si	2d	1/3	2/3	0.481(1) 1.22(7)
					0.479(1) 1.00(7)
	O1	2d	2/3	1/3	0.790(2) 1.39(8)
					0.783(2) 1.02(6)
	O2	6g	0.4739(9) 0.2966(8) 0.649(1) 1.39(8)
			0.4742(9) 0.2956(7) 0.647(1) 1.02(6)
	O3	6g	0.2174(7) 0.1001(6) 0.219(1) 1.39(8)
			0.2161(7) 0.0986(6) 0.223(1) 1.02(6)
		a ( Å)	c ( Å)		U
	8.51009(1)	5.251660(7)	0.0071(3)
	8.50665(1)	5.247587(8)	0.0122(4)

N. A. Hill, J. Phys. Chem. B104, 6694 (2000).

I., T. Peterlin-Neumaier and E Steichele,J. Phys. C: Solid State Phys.15, 4835 (1982)

N. Hur, S. Park, P. A. Sharma, J. S. Ahn, S. Guha and S-W. Cheong, Nature 349, 392(2004)

Acknowledgments

This work has benefited from the financial support of the ANR contract BLAN06-1 140756. We would like to thank J. Marcus, J. Debray and B. Kundys for their help during the dielectric measurements, and E. Ressouche for his contribution to the single crystal neutron diffraction experiment.


Simon Le Berre
email: simon.leberre@cea.fr

Isabelle Ramiere
email: isabelle.ramiere@cea.fr

Jules Fauque
email: jules.fauque@cea.fr

David Ryckelynck
email: david.ryckelynck@mines-paristech.fr

Isabelle Ramière

Condition number and clustering-based efficiency improvement of reduced-order solvers for contact problems using Lagrange multipliers

Keywords: Hyper-reduction, contact mechanics, Lagrange multipliers, solvability condition, condition number, clustering, local reduced order models

Introduction

Variational inequalities (VI) arising from contact mechanics problem are still nowadays computationally expensive to solve. The size of the discrete full-order model (FOM) becomes rapidly prohibitive when accurate solutions are required. The issue of high computational cost has been widely studied for variational equality problems, and a common solution consists in using reduced order modeling (ROM) techniques. But in case of contact VI problems, the application of such techniques is not straightforward. Indeed the non-smoothness of the contact constraints, especially the non-negativity of the contact pressure, remains a challenge for the adaptation of ROM methods.

A number of numerical methods, such as the penalty method or the Augmented Lagrangian (AL) formulation, consist in building unconstrained minimization problems for solving contact. In this particular framework, as the main challenges related to ROM adaptation are circumvented, various reduced-order methods have been successfully applied, see [START_REF] Balajewicz | Reduced order models for pricing european and american options under stochastic volatility and jump-diffusion models[END_REF]2,3] for the penalty method and [4, 5, 6] for the AL (or AL-like) formulation. However these kinds of contact formulations require numerical parameters adjustments to reach convergence. Concerning the penalty method, another disadvantage is that practically this method leads to some approximation of contact conditions (cf. interpenetration).

In counterpart, the Lagrange multipliers approach enable to satisfy precisely surface contact conditions. The resulting contact problem is a saddle-point problem with primal (displacement) and dual (contact pressure) variables. The Lagrange multipliers physically represent the contact pressures (opposite of the contact forces) and must remain non-negative. This constraint avoid the use of the famous Proper Orthogonal Decomposition (POD) snapshots method for the reduction of the dual basis. Hence, very few ROM methods have been applied to such problems. Their features mainly concern the construction of a reliable low-rank dual basis. We can distinguish approaches related to quasi-static contact problems to those related to dynamic contact problems. Concerning the reduction of (quasi-)static contact problems, we can cite the extension of the reduced bases method [7], where the dual reduced basis is obtained by snapshots selection. This may lead to low compressed reduced basis for large parametric space. In [8], a projected based method has been proposed with a Non-Negative Matrix Factorization (NNMF) to obtain the reduced dual basis. The NNMF satisfies the positivity constraint but doesn't give a high compression rate [9]. In [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] the dual reduced basis is built through a Cone-Projected Greedy algorithm resulting in a better compression rate. However, this approach has not yet been verified outside the training set. Another strategy has been proposed in [9], lying on Hyper-Reduction techniques. Such kind of techniques [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF]13,[START_REF] Ballani | A component-based hybrid reduced basis/finite element method for solid mechanics with local nonlinearities[END_REF]have shown their efficacy to reduce high-dimensional reconstruction required by non-linearity evaluations. Hence, in [9] an Hyper-Reduction method based on a reduced integration domain (RID) [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF] was adapted to contact problems. The low-rank dual basis is then simply defined as the restriction of the full-order dual basis to the RID, ensuring naturally the non-negativity constraint and leading to precise dual forces on the RID. Moreover, this method ensures the verification of the full contact conditions inside the RID. Concerning dynamic contact problems, the problem of building a reduced dual basis respecting the nonnegativity has been recently bypassed [16] through the use of the reduction of the primal variable to compute the Lagrange multipliers in an incremental way or [17] through the use of the Linear Complementarity Programming method where only the primal variable is reduced.

In this work we propose to improve the efficiency of the so-called Hybrid Hyper-Reduced (HHR) method, the hyper-reduced method based on the RID introduced in [9], especially for large parametric variation of the contact zone. The consistency of the HHR approach has been proven in [9], lying on the solvability condition. Since high primal basis compression by POD is performed, this condition may be not initially respected. It has been proposed in [9] to enrich the primal reduced basis with primal finite element (FE) shape functions while not respecting the solvability condition. However, the respect of the solvability condition is generally not sufficient to obtain accurate dual forces. No guideline was proposed in [9] for the choice of the FE shape functions to be added. Practically, for D-dimensional problems (D > 1) it was proposed to add all primal FE shapes function lying on the contact zone in the reduced primal basis. For high dimensional problems with large parametric variations, this strategy may compromise data compression. In order to automatically and efficiently extend the reduced primal basis, we first highlight in this contribution the strong relation between the projected contact rigidity matrix condition number and the precision of the dual solutions obtained by the HHR method. Taking this relation into account, we propose two primal basis enrichment strategies aiming at improving the condition number of the projected contact rigidity matrix while respecting the solvability condition. These enrichment strategies are both based on a greedy algorithm. There are very few elements in the literature studying the relation between condition number and precision of a ROM method. To the best of our knowledge, only [13] uses the condition number as a selection criteria for the sampling in the Missing Point Estimation (MPE) method. Secondly, as the enrichment strategy may increase significantly the size of the primal reduced basis, we propose a clustering strategy derived from [18] to treat contact problems with large parametric spaces. The clustering enables us to build local ROMs and hence to reduce high nonlinear (and non-smooth) problems in a piece-wise low-rank manner. The interested reader can refer to [18,[START_REF] Amsallem | Fast local reduced basis updates for the efficient reduction of nonlinear systems with hyper-reduction[END_REF]20,[START_REF] Peherstorfer | Localized discrete empirical interpolation method[END_REF]22] for examples of local ROMs applied to Hyper-Reduction schemes. Here, the local ROMs are obtained by applying on each cluster the HHR method for contact problems with condition number based enrichment introduced in this paper. The proposed strategies result in highly accurate dual (and primal) solutions with reduced basis of limited sizes even for large variations of the contact zone.

The rest of the article is structured as follows. A brief presentation of the FOM for Lagrange multipliers-based contact problems is witten in Section 2. The HHR model for contact mechanics is also recalled (and few adapted) in this section. In Section 3 we present the correlation between condition number of the projected contact rigidity matrix and precision of dual solutions. We then propose the two primal reduced basis enrichment strategies. In Section 4 the clustering framework for contact mechanics is introduced. Finally in Section 5 numerical experiments are performed allowing to appreciate the efficiency of the proposed strategies.

Contact problem, FE and HHR models

Static unilateral contact probem

Ω 1 Ω 2 d t 1 n 1 t 2 n 2 Γ C 1 Γ C 2 g Γ N g Γ N u 0 Γ D u 0 Γ D Figure 1: Contact model problem.
We are interested in the solution of elastostatic frictionless unilateral contact problems. The notations we consider are presented in Figure 1. Moreover, we classically denote by u N the relative normal displacement on the potential contact surface Γ C , d the initial gap and F N the normal force component. Then, the Signorini's law for unilateral contact is given on Γ C by: the non-penetration condition u N ≤ d, the non-adhesion condition F N ≤ 0 and the complementarity slackness (u Nd)F N = 0. Under the small strain hypotheses, the strong formulation of the elastostatic frictionless unilateral contact problem can be expressed as [23,[START_REF] Kikuchi | Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF][START_REF] Wriggers | [END_REF]:

                                 -div σ = f in Ω, σ = C : ε in Ω, ε = 1 2 grad u + (grad u) T in Ω, u = u 0 on Γ D , σn = g on Γ N , σn = F N n on Γ C , u N ≤ d, F N ≤ 0, (d -u N )F N = 0 on Γ C , (1)
with u the displacement field, σ the stress field, ε the strain field, f the prescribed volume density of forces, u 0 the imposed displacement on Γ D , g the imposed forces (surface density) on Γ N and n the external unit normal.

The tensor C denotes the fourth-order elasticity tensor.

The weak formulation derived from problem (1) induces a variational inequality (VI) [23]. This VI, equivalent to a minimization problem, can be classically treated by Lagrange multipliers [START_REF] Kikuchi | Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF]. The resulting weak formulation is then a saddle-point problem, cf. Eq (2).

       Find (u, λ) ∈ V × W such that a(u, v) + b(v, λ) = l(v) ∀v ∈ V 0 , b(u, q -λ) ≤ d, q -λ H 1 2 (Γ C ),H -1 2 (Γ C ) ∀q ∈ W, (2)
with λ the Lagrange multiplier,

W = {q ∈ H -1 2 (Γ C ) | q ≥ 0}, H -1 2 the dual of H 1 2 , b(v, q) = v N , q H 1 2 (Γ C ),H -1 2 (Γ C ) and •, • H 1 2 ,H -1 2 the dual product, V = {v = r(u 0 ) + ũ | ũ ∈ V 0 }, r(u 0 ) the lifting of u 0 over Ω, V 0 = {v ∈ (H 1 (Ω)) D | v = 0 sur Γ D } and a(u, v) = Ω ε(v) : C : ε(u) dΩ, l(v) = Ω f T v dΩ + Γ N g T v dΓ.

Finite element discretisation

The standard FE discretization [START_REF] Wriggers | [END_REF] of saddle-point problem (2) results in system (3).

           Find (U, Λ) ∈ R N × (R + ) N λ such that KU + B T Λ = F, BU ≤ D, Λ T (D -BU) = 0, (3)
with U being the discretization of the displacement solution in the primal FE basis, Λ the discretization of the Lagrange multipliers in the dual FE basis, N is the number of primal DoFs and N λ the number of potential contacts. The matrices or vectors

K ∈ R N×N , B ∈ R N λ ×N , F ∈ R N , D ∈ R N λ
are respectively the rigidity matrix, the contact rigidity matrix, the source vector and the initial gap vector. The discrete saddle-point problem (3) has to respect the classical Ladyzenskaia-Babuška-Brezzi (LBB) [START_REF] Babuška | The finite element method with lagrangian multipliers[END_REF]27] or inf-sup condition implying stability of the solution and hence solvability of the problem. For node-to-node contact with linear (or bilinear) finite elements approximation, which are of interest here, this condition is ensured [28].

Hybrid Hyper-Reduced model

2.3.1 Reduced integration domain for contact problems.

Ω A Ω B Ω A ∩ Ω B Γ C A Γ C I Γ C 1 Γ C 2 Ω 1 Ω 2 (a) Extended RID interface (Γ C I )
.

Ω A Ω B Ω A ∩ Ω B Γ C A Γ C 1 Γ C 2 Ω 1 Ω 2 (b) Extended RID (Ω A ).
Figure 2: Modification of the RID definition for the contact problems.

The RID [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF], that will be the support of the hyper-reduced discretized problem, is built from a set of DoFs associated to nodes located of the whole domain Ω. The DoFs selection method usually consists of a sampling method such as the DEIM [START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF] or an arbitrary selection for areas of interest. For contact problems, the RID definition must be adapted to ensure computable gaps on the RID. We focus here on node-to-node contact algorithm. In [9], an extension of the RID interface has been proposed, see Figure 2a.

The unpaired DoFs on the potential contact surface of the RID were hence removed from the set of available balance equations. In this paper, we propose an extension of the RID to treat contact on each node of the potential contact zone included in the RID, see Figure 2b. Denoting by L the list of node indices used to build the RID, the complementary contact list L reads

L = {node j on Γ C ; j / ∈ L | ∃i ∈ L such that (i, j) is a pair of contact nodes}. ( 4
)
The RID is now defined as:

Ω A = i∈L supp(ϕ i ) j∈L supp(ϕ j ) , (5)
with ϕ k , k = 1, ..., N the FE shape functions. We denote Ω B the counterpart of Ω A such that Ω = Ω A ∪ Ω B . The complementary list L plays the role of a zone of interest as defined in [START_REF] Giacoma | A multiscale large time increment/fas algorithm with time-space model reduction for frictional contact problems[END_REF]. This construction variation allows to deal with more potential contacts, but slightly enlarges the size of the RID.

Reduced bases

Concerning the primal basis reduction, the standard POD-snapshot method [START_REF] Fauque | Hybrid hyper-reduced modeling for contact mechanics problems[END_REF] is used. The size l of this reduced basis (RB) is automatically obtained according to the approximation error P OD on the snapshot matrix. The primal RB is then the matrix V ∈ R N×l containing the l first left singular vectors of the singular value decomposition (SVD) of the snapshots matrix. The discrete displacement is then projected in this RB U = Vγ, γ ∈ R l being degrees of freedoms (DoFs) associated to the primal RB , l N.

In order to avoid a non efficient compression of the dual basis, we use the RID concept to simply reduce the number of dual DoFs by restriction of the dual FOM basis to the RID. This strategy enable precise dual forces to be obtained on the RID, as shown in [9]. In counterpart, the reduced dual basis is of local support only. If the solution outside the RID is of interest, a reconstruction strategy based on the solution of a non-negative least squares (NNLS) problem using the dual snapshots as well as the ROM primal and dual solutions has been proposed in [9].

Hybrid Hyper-Reduced formulation

Complexity reduction of hyper-reduced models accounts for the sparsity of finite element matrices. We introduce the set A of primal DoFs not connected to Ω B :

A = i ∈ [[1,N]] | Ω B ϕ T i ϕ i dΩ = 0 , (6)
and the set I of DoFs connected to A by the rigidity matrix K (interface nodes in the discretization scheme sense):

I = j ∈ [[1,N]], j / ∈ A | ∃i ∈ A s.t. (K) ij = 0 . ( 7
)
Concerning the reduction of the contact conditions, we define the set A c of DoFs of A localised on the potential contact zone and the set A λ of dual DoFs connected to A c by the contact matrix B.

A c = {j ∈ A | ∃i ∈ [[1,N λ ]] s.t. (B) ij = 0}, (8)
A λ = {i ∈ [[1,N λ ]] | ∃j ∈ A c s.t. (B) ij = 0}. ( 9
)
The cardinal of A λ , denoted N C λ in the sequel, corresponds to the number of potential contact points in the RID and hence to the size of the dual reduced basis. The Hybrid Hyper-Reduced (HHR) model, introduced in [9], reads: [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] with γ the primal solution and Λ[A λ ] the dual solution of the reduced model. As shown on Eq. ( 10), the HHR formulation respects all the Signorini's conditions inside the RID.

           Find (γ, Λ) ∈ R l × (R + ) N C λ such that V[A, :] T K[A, A ∪ I]V[A ∪ I, :]γ + V[A c , :] T B[A λ , A c ] T Λ[A λ ] = V[A, :] T F[A], B[A λ , A c ]V[A c , :]γ ≤ D[A λ ], Λ[A λ ] T (D[A λ ] -B[A λ , A c ]V[A c , :]γ) = 0,
As Problem (3), Problem ( 10) is a saddle-point problem. It has been proven in [9] that the HHR formulation is consistent with the FE formulation if the solvability condition of Problem ( 10) is respected. In the setting of ( 10) it leads to have B[A λ,active , A c ]V[A c , :] of full row rank, A λ,active being the set of DoFs associated to an active contact (equality constraint), which is a priori unknown.

Solvability condition and enrichment of the primal reduced basis

Extended solvability condition

The set of active contact A λ,active being a priori unknown, we have to consider a stronger solvability-like condition. We impose the projected contact rigidity matrix including all potential contacts, i.e. B[A λ , A c ]V[A c , :], to be of full row rank. Since B contains as many lines as there are potential contact points (node-to-node contact), the extended solvability condition on the reduced model [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] imposes

rank(B[A λ , A c ]V[A c , :]) = N C λ . (11)
A necessary condition to verify Eq. ( 11) is that N C λ ≤ l, which is not always respected due to the RID construction. In order to respect condition [START_REF] Barrault | An 'empirical interpolation' method: application to efficient reduced-basis discretization of partial differential equations[END_REF] without modifying the RID, the primal RB has to be enriched.

In [9] it was proposed to enrich the primal RB with FE shape functions related to the nodes on the contact surface. In the case of Lagrange FE, the extended primal RB denoted V then reads:

V = V I N :, L C , ( 12
)
with I N the identity matrix of size N, L C ⊂ A C the node indices of FE shape functions to be added. The primal reduced basis dimension is then l = l + card(L C ). This enrichment is advantageous to accurately solve the solution on the contact surface, but the choice of L C may be not trivial. In [9] the enrichment was either done manually with an arbitrary selection (for one-dimensional problems) or it was proposed to add all primal FE shape functions of the potential contact zone in order to simply ensure the full row rank condition, good stability, and accurate solutions on the contact zone. In this latter case L C = A C , which may lead to a significant expansion of the primal RB.

One of the main contribution of this paper is to propose strategies to automatically and wisely obtain the set L C in order to optimize the ratio between the primal RB dimension l and the precision of the dual reduced solutions.

Importance of the projected contact rigidity matrix condition number

In the FE model, the condition number of the contact rigidity matrix B, written κ(B), is equal to 1. Accordingly to this observation, the more the condition number of the projected contact rigidity matrix of the HHR model, denoted κ(B[A λ , A c ]V[A c , :]), is close to 1, the more precise are the dual reduced solutions. On the other hand, the primal reduced solutions are less sensitive to this condition number. An illustration of this behaviour of dual and primal HHR solutions is reported in Figure 3, from the test case of Section 5.2. The relative errors considered are presented in Section 5.1, Equations ( 17) (primal) and ( 18) (dual).

FE enrichment of the primal RB

In order to make the condition number

κ(B[A λ , A c ]V[A c , :]) closer to κ(B)
, it seems natural to add to the primal RB V some FE shape functions located on the potential contact zone of the RID, cf. Eq [START_REF] Bader | An empirical interpolation approach to reduced basis approximations for variational inequalities[END_REF]. When

L C = A C , we then get κ(B[A λ , A c ]V[A c , :]) = κ(B). We propose to automatically obtain the set L C ⊂ A C through a greedy algorithm based on κ(B[A λ , A c ]V[A c , :]), see Algorithm 1. In this algorithm, A C (p)
V L C = ∅ ; k = 0 ; G = [[1,card(A C )]] ; κ min = κ(B[A λ , A c ]V[A c , :]) ; V = V; while (κ min > κ max or rank(B[A λ , A c ]V[A c , :]) < N C λ ) and k < card(A C ) do k = k + 1 ; for p ∈ G do V tmp = V I N [:, A C (p)] ; if κ(B[A λ , A c ]V tmp [A c , :]) < κ min then κ min = κ(B[A λ , A c ]V tmp [A c , :]) ; p min = p ; end end G = G \ {p min } ; L C = {L C , A C (p min )} ; V = V I N [:, A C (p min )] ; end
It is worth noting that even when the extended solvability condition is respected, the set L C could continue to be enlarged in order to improve the dual solutions.

POD modes and FE shape functions enrichment of the primal RB

The FE shape functions enrichment is a way to extend the primal RB, but it naturally comes to mind to use an enrichment strategy that uses more POD modes (less truncated SVD) for a given RID. In this way, more global information (from snapshots) can be included in the extended primal RB. Practically, this strategy is based on a second POD threshold P OD2 < P OD to obtain l 2 supplementary primal POD modes. The primal POD basis considered in the HHR model [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]

is now V ∈ R N×(l+l 2 ) .
The extension of the primal POD RB through the order of singular vectors does not guarantee neither the extended solvability condition nor a satisfying condition number of the projected contact rigidity matrix. This strategy has hence to be combined with FE shape functions enrichment (see Section 3.3 and Algorithm 1) to obtain accurate dual solutions.

The main advantage of this strategy is that it allows us to use a more precise POD RB without impacting the size of the RID. Indeed, in most cases, the RID construction is based on a selection method relying on the primal RB (running the DEIM algorithm on the primal RB for example). Hence, there is a strong connection between the size of the primal RB and the size of the RID. Using the proposed strategy, the RID construction is still based on the POD threshold P OD while the HHR model is based on the extended POD threshold P OD2 . Let us underline that P OD2 must remain large enough in order to avoid overfitting. Hence this strategy is mainly of interest for large parametric spaces.

Remark For both enrichment strategies, one has to verify the basis property of the extended matrix V on the RID. The verification has to be done in the offline part of the method, without impacting the online part.

Clustering

In case of highly nonlinear problems such as non-smooth contact problems, clustering strategies may offer the possibly to treat the nonlinearity in a linear piece-wise manner through low-rank (ROM) approximations on different clusters partitioning the parametric space [18]. For large parametric spaces, we then expect to improve the accuracy of the reduced solution for a given dimension of RBs. We then consider a partition of P composed of N s elements relying on the discretization P tr . On each element j = 1, ..., N s , an orthonormal basis V j is obtained from POD on the elementwise snapshots matrix. This local snapshots matrix is composed by snapshots associated to points in P tr belonging to element j.

From regular snapshots

Regular

As suggested in [18,[START_REF] Ryckelynck | Mechanical dissimilarity of defects in welded joints via Grassmann manifold and machine learning[END_REF], the Grassmann distance d Gr(∞,∞) between local RBs is then used for the clustering the simulated data for model order reduction. The dissimilarity matrix δ ∈ R M s ×M s writes:

∀j 1 , j 2 ∈ [[0,N s ]] 2 , (δ) j 1 ,j 2 = d Gr(∞,∞) (V j 1 , V j 2 ), (13)
with

d Gr(∞,∞) (A 1 , A 2 ) = π 2 4 |a 1 -a 2 | + min(a 1 ,a 2 ) ∑ i=1 α 2 i 1 2 , ( 14
)
where a 1 and a 2 are the dimension of bases A 1 and A 2 respectively. The angles α i are related to the singular values Σ ii , 1 < i ≤ min(a 1 , a 2 ), of A T 1 A 2 by Eq. ( 15)

Σ ii = cos(α i ). ( 15
)
The smaller (δ) j 1 ,j 2 the closer the subspace spanned by bases V j 1 and V j 2 separately. A clustering algorithm based on the dissimilarity matrix δ is finally applied to build K clusters from the initial partition of P into N s elements. In this work we use the k-medoids algorithm [31], to have a focus on compactness.

On each cluster k, k = 1, ..., K the local snapshots matrix gathers the snapshots associated to points of the training set P tr included into the cluster. From this local snapshots matrix, a local HHR model is built following the strategy described in Section 2.3. Especially, a local RID is considered in each cluster k. As the clustering aims to reduce the range of the parameter variations, an accurate low-rank approximation of the primal basis can be obtained. Hence, in this case only the enrichment strategy based on FE shape functions (see Section 3.3) has to be considered in order to ensure both the respect of the local solvability condition as well as accurate reduced dual forces.

Concerning the classification, a trivial strategy based on the clustering partition of the parametric space P is used. This classification is used online for the local ROM selection associated to given input data (generally outside the training set P tr ).

From irregular snapshots

The regular-based clustering strategy introduced in Section 4.1 may become prohibitive for high-dimensional parametric spaces. In such cases, an irregular distribution of snapshots in the parametric space has to be considered for the training set P tr , generally thanks to greedy algorithms. Then, it is no more suitable to rely on an element-wise partition of the parametric space for the clustering. A simple adaptation of the procedure of Section 4.1 is proposed here. The dissimilarity matrix is now based on the Grassmann distance between each normalized snapshot S j , j = 1, ..., M s , with M s = card(P tr ) the number of snapshots. In this case, Eq. ( 14)

reduces to d Gr(∞,∞) (S j 1 , S j 2 ) = arccos(S T j 1 S j 2 ), ∀j 1 , j 2 ∈ [[0,M s ]] 2 . ( 16
)
The clustering strategy returns in this case in labelling the snapshots.

Local ROMs are built similarly to Section 4.1 from the cluster-based local snapshots matrix.

The simplification of the clustering phase comes with adding complexity in the classification construction required for the online local ROM selection. Advanced classifiers, see for example [18], have to be considered as no clustering partition of the parametric space is available. In this work, we propose to build offline a supervised learning classification method trained on the snapshots and their cluster label. In the numerical examples, the Random Forest classification method [START_REF] Breiman | Random forests[END_REF]33] has been performed.

Numerical results

In all the following numerical experiments, the FOM solution is obtained thanks to the Cast3M finite elements software (2020 version) [START_REF]Cast3M[END_REF]. Lagrange P1 or Q1 finite elements have been used, with 3 or 4 integration points respectively for strain/stress calculation.

Test case a: Half disks of Hertz, unidimensional parametric space

This test case is described in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF]. It considers two half disks in a two-dimensional setting, each disk having the same radii R = 1 m. The two disks are separated by an initial gap being at minimum d = 0.1 m. Both disks have the same material characteristics: a Young modulus E = 15 Pa and a Poisson coefficient ν = 0.35. Dirichlet boundary conditions are imposed on the horizontal edges of the disks, with u x = 0 and u y = ±µ with an upward movement for the lower disk and a downward movement for the upper disk. The range values of µ defines the parametric space. A representation of the initial state of the test case can be seen in Figure 5a.

The Finite Element full-order model (FOM) is built on a mesh with 678 nodes. The circular part of both disks, constituting the potential contact zone, contains 55 nodes on each side. The deformation obtained with the FOM for µ = 0.3 is presented in Figure 5b. We consider a parametric space P = [0.15, 0.45] (m) for µ and build the training set P tr = {0.15 + 0.01i 0 ≤ i ≤ 30} (m). Hence card(P tr ) = 31. To assess the accuracy of the methods on this test case we will use three kind of errors, the primal relative error in 2-norm Eq. ( 17), the dual relative error in 2-norm Eq. ( 18) and the error on the minimum energy used in [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF] Eq. ( 19)

e u (µ) = ||u ROM (µ) -u FOM (µ)|| 2 2,Ω ||u FOM (µ)|| 2 2,Ω 1 2 , ( 17
)
e λ (µ) =   ||λ ROM (µ) -λ FOM (µ)|| 2 2,Γ C A ||λ FOM (µ)|| 2 2,Γ C A   1 2 , ( 18
)
e ener (µ) = 1 2 |a(u ROM (µ), u ROM (µ)) -f (u ROM (µ)) -a(u FOM (µ), u FOM (µ)) + f (u FOM (µ))| , (19)
with (u ROM (µ), λ ROM (µ)) and (u FOM (µ), λ FOM (µ)) denoting the ROM (here obtained with the HHR model) and FOM solutions (primal and dual) respectively. Let us underline that in the case of the HHR method, the error on the Lagrange multipliers (18) is naturally evaluated on the RID. The extension from Γ C A to Γ C could however being performed through the reconstruction strategy introduced in [9]. On the counterpart, the extension to Ω of the displacement field is obvious thanks to the POD basis. Hence, the primal error will always be evaluated on the whole domain Ω.

This test case has the particularity to present a double symmetry (horizontal and vertical) for the dual and the primal solutions. This symmetry added to the relatively small contact surface and the one-dimensional parametric space leads to a test case relatively easy to treat for ROMs. Using the HHR, a small RID is enough to give very good results.

In the following results, the RID has been built from the set L obtained applying the DEIM to the primal POD reduced basis. This set is supplemented by two nodes on the top and bottom surfaces to be able to treat Dirichlet boundary conditions. As the parametric space is quite small, an accurate error threshold of P OD = 10 -8 is used resulting in a primal RB with 9 vectors (l = 9). The RID is reported in Figure 6a. It results in 16 nodes on the contact surfaces (8 on each side), hence N C λ = 8. The solvability condition (cf. Eq. ( 11) is in this case directly respected. Moreover, the condition number is also initially very good κ(B[A λ , A c ]V[A c , :]) = 3.4. Hence, this test case does not require any enrichment of the primal RB to obtain accurate results with the HHR method. In Figure 6b, the Lagrange multipliers for the FE and HHR models are reported for different values of µ / ∈ P tr . This To further evaluate the accuracy the method we compute the errors e u (µ), e λ (µ) and e ener (µ), either on the RID and on Ω (with the reconstructed dual solution) on a set of 101 points in P (including the 11 points of the training set P tr ). The results are reported in Figure 7. The primal solution is always approximated with accuracy, with relative errors in 2-norm below 0.05% (5.10 -4 ) on all the parametric space. The dual solution is also approximated accurately ranging from 0.1% to a maximum of 0.8%. The error on the minimum energy is very also satisfying, with errors below 10 -4 on Ω, even outside the training set.

In [START_REF] Benaceur | A reduced basis method for parametrized variational inequalities applied to contact mechanics[END_REF], the CPG (Cone-Projected Greedy) algorithm has been applied on the same problem with the same training set and the same POD truncation (l = 9 POD modes). The CPG method yields to errors ranging from 10 -6 to 10 -3 for the error on the minimum energy, and from 10 -5 to 4.10 -3 for the error on the displacement field. Those errors were obtained on parametric points of the training set only. No errors on Lagrange multipliers are reported. Compared with the errors shown in Figure 7 of the present paper, it highlights the robustness of the proposed HHR method, that is to say the accuracy of the results on all the parametric space even outside the training set.

Test case b: 2D axisymmetric contact -bidimensional parametric space

We consider here a two-dimensional axisymmetric test case introduced in [START_REF] Liu | On the coupling of local multilevel mesh refinement and zz methods for unilateral frictional contact problems in elastostatics[END_REF], representing a simplified contact problem between a nuclear fuel pellet and its surrounding cladding. The test's geometry is sketched in Figure 8a. The two materials follow an elastic isotropic behaviour with a Young modulus and Poisson coefficient of E = 190 GPa, ν = 0.3 for the fuel, and E = 78 GPa, ν = 0.34 for the cladding. The uniform initial gap between the two solids is d = 2 µm. The test case consists in applying a pressure P 1 on a length h 1 on the outside part of the cladding. The range values of parameters P 1 and h 1 form the parametric space P = [5, 200] MPa×[0.6, 6.0] mm. This configuration leads to a low displacement along the z axis, allowing a good modeling of the contact with a node-to-node configuration. The parametric variations are challenging for ROM models as it includes strong variations of the contact surface as well as the contact pressure. The FOM is built on a mesh with 7070 nodes, 6161 for the mesh representing the fuel and 909 for the cladding, see Fig. 8b. There are 101 potential contacts (N λ = 101).

In order to build the HHR model, 500 snapshots regularly distributed in the parametric space are used for the training set P tr . Due to the large variations of the parametric space, a POD threshold of P OD = 10 -6 is used to initially limit the dimension of the primal RB to l = 36 (see the decrease of approximation error on the snapshot matrix in Figure 16). Using the set of points L obtained from the DEIM applied to the primal RB results in the RID illustrated in Figure 9a where a region of interest corresponding to the Neumann surface has been added. The number of potential contact is the RID is N C λ = 12. This RID gives ride to unsatisfactory primal solutions. We then expand the set L using a second DEIM on a dual POD RB (with a threshold = 4 × 10 -4 ). It worth underlying that this dual POD RB is only used for the RID expansion. The expanded RID is presented in Figure 9b, having also the Neumann surface as zone of interest. In this case N C λ = 34. The advantage of using the DEIM on a dual POD RB is that the selected DoFs will be located on the contact surface. In the following the results presented focus on the expanded RID (l = 36, N C λ = 34) and on parameters outside of the training set P tr . In this case, the solvability condition is not initially respected. Hence, the primal RB have to be enriched to solve the HHR problem.

Enrichment strategies

We focus here on the algorithm presented in Section 3.3 for the enrichment of the primal RB with FE shape functions. Figure 10 show the output at the different iterations of the greedy algorithm 1. Although the solvability condition isn't verified initially, a numerical solution can be obtained. However this latter exhibits dual relative errors. Figure 10 illustrates that the strict respect of the solvability condition (indicated by the dashed red vertical lines) isn't enough to obtain a satisfying accuracy on the dual solutions. With FE shape functions enrichment, the solvability condition is verified for l = 40 (4 supplementary FE shape functions), but with a condition number of B[A λ , A c ]V[A c , :] superior to 1000, lying to dual relative errors greater than 70%. The further enrichment of the primal RB based on the improvement of this condition number enables to drastically improve the dual relative errors, up to 2% for the parametric point (170 MPa, 3.2mm) considered in Figure 10a. The dual relative error becomes significantly lower for a condition number around 10, see Figure 10b. For κ max = 10, the greedy algorithm results in card(V)=55 and a dual relative error of 6.3%. The enrichment strategy of the primal RB with POD and FE modes, introduced in Section 3.4, is now compared to the enrichment strategy based on FE modes only, and previously illustrated. We set κ max = 1 (maximal enrichment) in Algorithm 1 for both approach. We report in Table 1 the results obtained for different truncature levels for the enrichment by POD modes on two points of the parametric space (strong and weak contact) outside the training set P tr .

The enrichment method by POD+FE modes can enable us to reach dual and primal relative errors lower

P 1 = 170MPa P 1 = 30MPa h 1 = 3.2mm h 1 = 5.25mm Rel. err.
Rel. err. nb POD than the enrichment method with only FE modes. Using a maximal enrichment κ max = 1, the enrichment method by extended POD naturally leads to higher cardinal of the primal extended RB (see Table 1).

P OD2 κ(B[A λ , A c ]V[A c , :]) l = card(V) mod.
Regarding the errors, we can see that for a parametric point already well treated by the FE enrichment method (see P 1 = 170 MPa and h 1 = 3.2 mm), few additional POD modes are enough to lower the dual relative error. However, when the dual relative errors remain high with the FE enrichment method, it is useful to add more POD modes to lower the dual relative error. For example, on the point P 1 = 30 MPa and h 1 = 5.25 mm, adding 9 POD modes (corresponding to P OD2 =3.5 10 -7 ) allows us to divide the dual relative error by four. We can outline from Table 1 that the choice of P OD2 remains not trivial as a too low values lead to phenomena associated with over-fitting that degrade the accuracy of the solutions.

Let us notice an enrichment method with POD modes only will not be efficient to lower the condition number of the projected contact rigidity matrix without impacting the precision of the reduced solutions.

In Figure 10, the output of the greedy algorithm on the primal RB extended with 9 supplementary POD modes (denoted 'l+9' in the legend) are reported. For an equivalent primal RB dimension, the POD+EF enrichment strategy becomes interesting only for low errors, see Figure 10a. However regarding to condition number of the projected contact rigidity matrix, the POD+FE enrichment strategy is quite always improving the dual errors for a given prescribed condition number, see Figure 10b. It is worth noting that the solvability condition is reached for a comparable condition number. This confirms the strong relation between solvability, precision and condition number.

The two enrichment methods are now compared on the whole parametric space. For that purpose, HHR solutions are calculated on 6400 points regularly distributed in the parametric space. We set the maximum condition number to κ max = 10 which seems to be a good compromise between precision and dimension of the extended primal RB.

The maps of the dual relative errors are reported in Figure 11 for the FE enrichment strategy (Fig. 11b) and POD+EF enrichment strategy with 9 suuplementary POD modes (Fig. 11a). It clearly shows an improvement due to the POD+EF enrichment method. The dimension of the final extended primal RB ( l) is in this case 7 modes larger (instead of 9), which implies that the improvement of the condition number is slightly facilitated by the addition of POD modes.

Whatever the enrichment strategy, it remains high dual relative errors on the left side of the maps, for low pressures P 1 . This highlights the complexity to efficiently reduced a problem on such parametric space, where the contact can be strong, weak or even absent, and on an large or small surface. This naturally open the way to study clustering strategies on the parametric space.

Before presenting numerical results for the clustering strategies, we first study here the impact of the size of the RID over the enrichment strategies. We hence consider a larger extended RID with more potential contact points. This RID is simply obtain considering a higher dual POD RB threshold of = 1 × 10 -4 for the DEIM point selection, see Section 5.2. In this case, the RID contains N C λ = 60 potential contacts with always the same primal RB of l = 36 POD modes, see Figure 12. As a consequence of the augmentation of N C λ , the solvability condition requires a larger enriched primal RB. For this RID, the FE enrichment method reaches the solvability condition for l = 61, that is to say 35 supplementary FE modes. But as previously shown, for this extended basis, the dual reduced solutions are not enough satisfying. Further enrichment based on the condition number of the projected contact rigidity matrix is required.

In Figure 13, we report dual relative error maps using both FE and POD+FE enrichment with a maximal condition number κ max = 10. This results in a primal RB with a cardinal of l = 83 for the FE enrichment method and of l = 91 for the POD+FE enrichment method with 9 supplementary POD modes (l + 9). Compared to Figure 11a, the first conclusion to be drawn is that the increase of N C λ does not improve the error maps for the FE enrichment method, see Figure 13a. In counterpart, we can see a significant improvement for the POD+FE enrichment method, cf. Figures 11b and13b . An important point here is to remember that the RID contains more Lagrange multipliers when N C λ = 60 than when N C λ = 34, meaning that the computation of the dual relative error as defined by (18) implies more forces in the case N C λ = 60. Contact forces for the parametric point P 1 = 170MPa and h 1 = 3.2mm for both configuration of N C λ with full FE enrichment (κ max = 1) are reported in Figure 14. In both situations, the dual forces are locally very accurate on the RID (2.6% for N C λ = 60, and 2.5% for N C λ = 34). However, having more potential contact points in the RID, and therefore more contact forces solved is obviously an advantage both for the physical interpretation of the simulation and the reconstruction of the dual solution, but the drawback is the dimension of the primal RB that grows with N C λ . This motivates furthermore to look forward clustering approaches to limit the increase of the primal RB cardinal while reaching the same accuracy.

Clustering strategies

From regular snapshots We focus here on the clustering strategy proposed in Section 4.1. The same regular training set P tr composed by 500 snapshots is used. The partition of parametric space is built through groups of 4 snapshots forming rectangles in the 2D parametric space representation, see Figure 15a. In Figure 15b we can see the results of the clustering algorithm for K = 3 clusters, each point representing a snapshot and the background color representing a cluster. The clusters contain respectively 170, 194 and 170 snapshots in their local snapshots matrices. Figure 15b shows that the clustering strategy on the parametric space naturally focuses on the position h 1 of the potential contact zone. It confirms that a large variation of the potential contact zone is the preponderant parameter of nonlinearity.

As the parametric space is locally reduced into the cluster, the local HHR models are obtained from a primal POD threshold P OD = 10 -8 . The choice of this threshold still enables a low-rank approximation as confirmed by Figure 16. The primal RB size in each cluster is l = {31, 39, 36}. A local RID is built for each local model. As previously a DEIM on the primal RB and a contribution from a DEIM on a dual POD RB (with a threshold = 10 -6 ) is used to obtain the set of nodes L. The resulting RIDs are reported in Figure 17. The number of potential contacts treated in each RID is N C λ = {36, 40, 35}. As a direct results of the clusters shapes, the RID of the local reduced order models are much more localized. Figure 17 which represents the zone of status change (contact/gap).

On each cluster, the primal RB extension is performed through the FE enrichment strategy only (since P OD is already low). The results obtained with a maximum enrichment (κ max = 1) for the two points studied in Table 1 are reported in Table 2. For the point P1 = 170MPa and h 1 = 3.2mm, already accurately treated by the enrichment strategies performed on the whole parametric space, similar dual errors are obtained (with a significant improvement on the primal error). For the point P1 = 30MPa and h 1 = 5.25mm the results are much better than the ones obtained on the whole parametric space whatever the enrichment strategy (see Table 1). Moreover, the size of the primal RB remains limited ( l = 67).

Cluster κ(B[A

λ , A c ]V[A c , :]) l = card(V) l = card(V)
Generally speaking, the proposed clustering strategy helps smoothening the dual relative errors. It particularly improves the accuracy on the parametric points where the relative errors are high without the clustering strategy, see Figure 18. Concerning the RB dimensions, this clustering strategy leads to comparable sizes than the strategies introduced in Section 5.2.1. However the dual errors are greatly improved by the clustering compared to Figure 11. The error maps of Figure 18 are even better than the ones obtained on an extended RID with the POD+EF enrichment strategy, see Figure 13b. In this latter case, the bases dimensions were almost 40% larger.

We can underline, that even with the clustering strategy, it remains difficult to accurately approximate very low contact pressure parametric zones. The number of clusters K can not be too much increased. Indeed, the more cluster the smaller the local parametric space becomes. But each cluster must have enough snapshots to build a representative POD RB and then an accurate HHR model. More snapshots could be initially used, but another interesting way to investigate would be to consider a recursive clustering strategy (with possibly adding snapshots) inside clusters having the least accurate HHR solutions.

From irregular snapshots

We focus here on an irregular distribution of the snapshots in the parametric space P obtained through a snapshot selection by the greedy algorithm introduced in [9]. The resulting training set P tr contains 200 snapshots, which can be visualized in Figure 20a. Let us notice that they are mainly located around small values of contact pressure P 1 .

In order to appreciate the solutions obtained by the clustering algorithm introduced in Section 4.2, we first briefly present some key results obtained by the HHR method with enrichment trained on irregular snapshots of P tr (without clustering). The HHR model is built using primal POD thresolhold P OD = 10 -6 , resulting in a primal RB of dimension l = 33. The RID is obtained as the one shown in Figure 12 (DEIM on the primal RB, and a contribution from a DEIM on a POD dual RB obtained with the treshold = 10 -4 ). It finally contains N C λ = 57 potential contacts. For the POD+FE enrichment method an error threshold P OD2 = 3.5 × 10 -7 is used (as in Figs. 13b and11b) resulting in a primal RB of dimension (l + 5) = 38 before enrichment. Both enrichment strategies are moreover applied setting κ max = 10 as a maximum condition number in Algorithm 1. The FE enrichment method results in a primal RB of cardinal l = 60, while the POD+FE method results in a primal RB of cardinal l = 82. We present in Figure 19 the dual relative errors on the parametric space for both reduced models. The errors are in good agreement with those obtained with 500 regular snapshots on a RID of similar size, see Figure 13. It is worth underlying that better results are globally reached (especially for low contact pressure) and with a more compact primal RB (lower dimension). This confirms the interest of greedy-selected snapshots [9].

In order to evaluate the clustering strategy based on irregular snapshots (cf. Section 4.2) on this training set, we build K = 3 clusters through the k-medoid algorithm (as for the previous clustering strategy based on regular snapshots). In Figure 20a, the cluster labels are reported. The Random Forest classifier is then trained to obtain a prediction function, represented in Figure 20b as a partition of the parametric space. Let us underline, that in this case the clusters are not uniquely defined thanks to the contact position h 1 but also according to the pressure values P 1 for small h 1 . It confirms that a small contact pressure applied to a reduced contact area is difficult to capture. The dual errors map is shown in Figure 22. The clustering strategy helps smoothening the dual relative errors, as observed previously. The errors are considerably reduced for RB of comparable sizes to those without clustering (cf. Figure 19), even for low contact pressures.

As improvement, we can see that an adaptive snapshots strategy into cluster could be interesting since big clusters (see cluster 1) may lack of snapshots and since the cluster delimitation could be improved where snapshots are sparsely distributed.

Conclusion

Efficient strategies to improve the dual solution precision of reduced order modeling of contact mechanics problems treated by Lagrange multipliers have been proposed. We have focused on the Hybrid Hyper-Reduced (HHR) method where the dual reduced basis returns to a restriction on a Reduced Integration Domain of the full order dual basis. First, the strong correlation between the condition number of the projected contact rigidity matrix and the accuracy of the dual solution has been highlighted. Two strategies of enrichment of the primal POD reduced basis have been introduced. These strategies allows us to verify the solvability condition of the reduced saddle-point problem (if not respected initially) but also to improve the accuracy of the method, especially the dual solution. Secondly for large parametric space, especially for large variation of the contact zone where the non-smoothness of the contact conditions makes difficult low rank dual approximations, clustering-based strategy have been proposed. Precise local enriched HHR models are then built leading to globally enhance the accuracy of the reduced model using limited sizes of reduced bases. Clustering strategies have been presented for regular or irregular snapshots. These different strategies have been applied on 2-dimensional test cases demonstrated the very good accuracy of the method on small parametric spaces, and the usefulness contribution of the enrichment and clustering strategies on large parametric spaces. The HHR method is able to accurately treat non-smooth cases where the contact forces can be strong or weak, and very localized or extended. Future works will focus on the extension of the HHR to more complex geometries where node-to-segment contact detection is required.

Figure 3 :

 3 Figure 3: Impact of the condition number of the projected contact rigidity matrix on the dual and primal relative errors.

Algorithm 1 :

 1 denotes the DoF number associated to the p th element of A C . Greedy algorithm based on the condition number of the projected contact rigidity matrix to enrich the primal RB. Data: maximum condition number κ max , matrices B and V, the set A C (primal FE DoFs on the potential contact zone of the RID) Result: The set of FE DoFs to activate L C ⊂ A C and the primal enriched RB

  Snapshots Partition of the parametric space POD on the partitionbased snapshots matrices Dissimilarity matrix based on the Grassmann distance Clustering Local ROMs

Figure 4 :

 4 Figure 4: Clustering framework

  (a) Initial state. (b) Deformation for µ = 0.3.

Figure 5 :

 5 Figure 5: Half-disks test case. FOM configurations

  Figure confirms the very good local accuracy on the Lagrange multipliers.

  (a) Reduced integration domain (RID).

  Lagrange multipliers for different values of µ / ∈ P tr .

Figure 6 :

 6 Figure 6: Half-disks test case. HHR model : RID and Lagrange multipliers. l = 9, N C λ = 8.

  Error on the minimum energy e ener (µ).

  Primal error on the displacement field e u (µ) on Ω and dual error of the Lagrange multipliers e λ (µ) on the RID.

Figure 7 :

 7 Figure 7: Half-disks test case. HHR errors evaluations of 101 points. l = 9, N C λ = 8

  Geometry of the 2D axisymmetric test case. Maillage initial (b) FOM mesh

Figure 8 :

 8 Figure 8: 2D axisymmetric test case. Sketch of the problem and FOM discretization.

  Figure 9: 2D axisymmetric test case. Reduced integration domain (RID)

  Dual relative error regarding the number of FE modes added to the primal POD RB enriched ('l+9') or not ('l').

  Relative dual errors regarding the condition number of the projected contact rigidity matrix for the two enrichment strategies.

Figure 10 :

 10 Figure 10: 2D axisymmetric test case. Study of the two enrichment strategies: by FE modes only ('l') and by POD+FE ('l+9', 9 POD modes added), with l = 36 and N C λ = 34. Point: P 1 = 170MPa and h 1 = 3.2mm

  FE enrichment. l = 55. POD+FE enrichment (l+9). l = 62.

Figure 11 :

 11 Figure 11: 2D axisymmetric test case. Dual relative errors maps (in %) for the HHR model with l = 36, N C λ = 34 on a test space of 6400 points. Enrichment strategies applied with κ max = 10.

Figure 12 :

 12 Figure 12: 2D axisymmetric test case. Larger extended RID. l = 36, N C λ = 60

  POD+FE enrichment (l+9). l = 91.

Figure 13 :

 13 Figure 13: 2D axisymmetric test case. Dual relative errors maps (in %) for the HHR model with l = 36, N C λ = 60 on a test space of 6400 points. Enrichment strategies applied with κ max = 10.

  (a) HHR with l = 36, N C λ = 34 and l = 70. (b) HHR with l = 36, N C λ = 60 and l = 96.

Figure 14 :

 14 Figure 14: 2D axisymmetric test case. Contact forces for P 1 = 170MPa and h 1 = 3.2mm. HHR with FE enrichment (red arrows) and FOM (black arrows). (forces superposed on the right side).

  (a) Partition of the parametric space. (b) Clusters for K = 3.

Figure 15 :

 15 Figure 15: 2D axisymmetric test case. Partition and clustering of the parametric space from regular snapshots.

Figure 16 :

 16 Figure 16: 2D axisymetric test case. Approximation error on the snapshot matrices for the case with (K=3) and without clustering.

35 )Figure 17 :

 3517 Figure 17: 2D axisymmetric test case.

Figure 18 :

 18 Figure 18: 2D axisymmetric test case. Dual relative errors map (in %) for the clustering strategy with K = 3 clusters. Test space of 6400 points. FE enrichment method with κ max = 10. HHR models with l = {31, 39, 36}, N C λ = {36, 40, 35} and l = {54, 67, 57}

  FE enrichment. l = 60.

  POD+FE enrichment (l+5). l = 82.

Figure 19 :

 19 Figure 19: 2D axisymmetric test case. Irregular snapshots in the training set. Dual relative errors maps (in %) for the HHR model with l = 33, N C λ = 57 on a test space of 6400 points. Enrichment strategies applied with κ max = 10.

  (a) Cluster labels on the irregular snapshots. (b) Prediction function obtained with the Random Forest classifier.(label 0:blue, label 1:grey, label 2:red).

Figure 20 :

 20 Figure 20: 2D axisymmetric test case. Clustering strategy for irregular snapshots. K = 3 clusters

Figure 21 :

 21 Figure 21: 2D axisymmetric test case. Irregular snapshots. RIDs associated to the local ROMs.

Figure 22 :

 22 Figure 22: 2D axisymmetric test case. Irregular snapshots. Dual relative errors map (in %) for the clustering strategy with K = 3 clusters. FE enrichment method with κ max = 10. HHR models with l = {27, 35, 23}, N C λ = {36, 56, 35} and l = {52, 77, 48}

Table 1 :

 1 2D axisymmetric test case: study of the enrichment strategy by POD modes and FE shape functions (l = 36 et N C λ = 34)

Table 2 :

 2 2D axisymmetric test case. Study of the clustering strategy with full local FE enrichment.

Acknowledgements

The authors are grateful to the PLEIADES project, financially supported by the CEA, EDF (Électricité de France) and FRAMATOME, that funded this research work.

A Approximation error

The approximation error is obtained by using the energy rate carried by the k first singular vectors of the SVD on the snapshot matrices: