V Keryvin

L Charleux

R Hin

J.-P Guin

J.-C Sanglebooeuf

Mechanical behaviour of fully densified silica glass under Vickers indentation

Keywords: Non-metallic glasses (silicates), Plasticity, Constitutive modelling, Densification, Nanoindentation

When subjected to constrained deformation loadings, such as during an indentation test, silica glass experiences complex deformation mechanisms including densification and volume-conservative shear plasticity. The densification mechanism may increase the density up to more than 21%. The question of the mechanical behaviour of an already fully densified glass sample naturally arises. This issue is one of the key points to address when one tries to propose a constitutive model of pristine silica glass. Indeed, during an indentation test, which is a popular test for this task, beneath the indenter tip, some regions might be fully densified. What is their behaviour, after saturation in densification and during loading, is therefore an issue to address.

Moreover, this is a crucial point for exploring the transition from plasticity to cracking, which is of paramount importance and a long-term objective for predicting the lifetime of glass products subjected to contact loadings such as impact or scratching.

In this paper, a quantitative identification of fully densified silica mechanical constitutive behaviour is made by using instrumented indentation testing and finite element analyses. The use of such an indirect method to assess the 1

Introduction

Silica glass (SiO 2 ) may be viewed as the archetype of a three dimensional network glass former. It is a highly brittle material that breaks in bending under few tens of MPa because of its extreme sensitivity to surface damage [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF]2]. Yet, under constrained mechanical loadings such as hydrostatic tests or indentation, it can deform permanently without breaking or even cracking [START_REF] Rouxel | Poisson's Ratio and the Glass Network Topology -Relevance to High Pressure Densification and Indentation Behavior[END_REF][START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF]. Under pure hydrostatic conditions, its density (respectively its volume) may increase permanently [START_REF] Bridgman | Effects of Very High Pressures on Glass[END_REF][START_REF] Cohen | Effects of Ultra high Pressures on Glass[END_REF][START_REF] Mackenzie | High-Pressure Effects on Oxide Glasses: II, Subsequent Heat Treatment[END_REF] by more than 20 % (respectively

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

decrease by more than 16 %) [START_REF] Rouxel | Poisson's Ratio and the Glass Network Topology -Relevance to High Pressure Densification and Indentation Behavior[END_REF][START_REF] Deschamps | Permanent densification of compressed silica glass: a Raman-density calibration curve[END_REF]. This has been related to some changes in the intermediate range order, among which the decrease of intertetrahedra angles [START_REF] Hehlen | Inter-tetrahedra bond angle of permanently densified silicas extracted from their Raman spectra[END_REF], and is referred to as pressure-induced densification (PID). It is now recognised that during indentation, densification occurs as well as volume conservative shear flow [START_REF] Rouxel | Indentation deformation mechanism in glass: Densification versus shear flow[END_REF]. These main mechanisms at stake have been observed separately and very recently properly modelled: PID as recalled previously [START_REF] Keryvin | Constitutive modeling of the densification process in silica glass under hydrostatic compression[END_REF] and shear flow by uniaxial compression experiments [START_REF] Lacroix | Plastic deformation and residual stresses in amorphous silica pillars under uniaxial loading[END_REF][START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF]. For the former mechanism (PID), the authors proposed a description of the densification process with a threshold pressure for the onset of densification, an increase of densification upon pressure, and another pressure threshold for the saturation in densification. They also pointed out the necessity to use a Finite Deformation framework and to account for the elastic stiffening with densification (for instance the bulk modulus will double) [START_REF] Keryvin | Constitutive modeling of the densification process in silica glass under hydrostatic compression[END_REF]. As for the latter mechanism, the authors of Ref. [START_REF] Lacroix | Plastic deformation and residual stresses in amorphous silica pillars under uniaxial loading[END_REF][START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF], using compression tests on micropilars, showed that silica glass behaves in compression as a perfectly elasto-plastic material (without any strain hardening in shear) with a compressive yield strength around 6-7 GPa. These two mechanisms take place together during an indentation test [START_REF] Ernsberger | Role of Densification in Deformation of Glasses Under Point Loading[END_REF][START_REF] Peter | Densification and flow phenomena of glass in indentation experiments[END_REF] so that their kinetics must be related to a combination of driving forces (shear and pressure). Different authors have proposed such constitutive modellings that are successful in describing the mechanical response of the instrumented indentation test (the force-displacement curve) [START_REF] Lambropoulos | Constitutive Law for the Densification of Fused Silica, with Applications in Polishing and Microgrinding[END_REF][START_REF] Xin | Densification of Fused Silica : Effects on Nanoindentation[END_REF][START_REF] Kermouche | Mechanical modelling of indentation-induced densification in amorphous silica[END_REF][START_REF] Gadelrab | Densification modeling of fused silica under nanoindentation[END_REF][START_REF] Keryvin | Densification as the Only Mechanism at Stake during Indentation of Silica Glass?[END_REF]. Some of them were also successful in observing densification fields beneath the residual imprint by Raman spectroscopy [START_REF] Perriot | Champagnon, Raman Microspectroscopic Characterization of Amorphous Silica Plastic Behavior[END_REF][START_REF] Kermouche | Mechanical modelling of indentation-induced densification in amorphous silica[END_REF] or by a chemical dissolution technique [START_REF] Niu | Locally enhanced dissolution rate as a probe for nanocontact-induced densification in oxide glasses[END_REF].

However, a piece of this jigsaw puzzle is still missing: what does the behaviour of the material become when the densification is saturated to ∼ 21%

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

? In other words, what is the mechanical behaviour of fully densified silica glass? This issue is of paramount importance not only from a condensedmatter physics point of view but also since the ultimate target of modelling the mechanical response of glasses to surface damage is to describe the transition from plasticity to cracking and it is highly liable that these densified or fully densified zones play a major role in it. This is the current objective of this paper. We will analyse instrumented indentation tests, in terms of force-displacement curves and residual imprints, to propose an adequate description of the mechanical behaviour of fully densified silica glass under Vickers indentation.. The paper is organised as follows. We will first describe how silica glass samples are densified up to more than 21%, how indentation tests are performed and how results are analysed.

The numerical procedures, using both two-and three-dimensional finite element analyses of the Vickers indentation test, will be then described. After presenting and analysing both experimental and numerical results, we will eventually challenge our findings in terms of mechanical modelling and material properties to very recent experimental and simulation literature results.

Experimental and numerical methods

Material and experimental procedures

Silica glass (Vitreosil TM , Saint-Gobain, France) specimens (0.6 mm in diameter and 2 mm in length) were densified by means of an octahedral multianvil apparatus using a Walker cell. A typical run consisted in raising the load pressure of the main ram at a rate of 0.5 MPa/min of oil up to obtaining a pressure of 25 GPa on the sample. After reaching the target pressure, the specimens were maintained at high pressure for one hour and then slowly

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
unloaded. Density was measured with a better than 0.001 g/cm 3 accuracy by means of a density gradient method using two partially miscible liquids (iodobenzene and methylene iodide). The density of the pristine samples was 2.2 g/cm 3 . After unloading from the Walker cell, samples had a density increased by 21.6 %. They are further referred to as fully densified (FD).

Further details may be found elsewhere [START_REF] Rouxel | Poisson's Ratio and the Glass Network Topology -Relevance to High Pressure Densification and Indentation Behavior[END_REF].

Instrumented indentation tests (IIT) were carried out with a micro-indenter testing device (Fischerscope H100 XYp, Fischer, Germany) at ambient conditions (23˚C and 55% relative humidity). It has a load resolution of 0.02 mN and a depth resolution of 2 nm. The calibration of the instrument was done by using ISO-14577 standard on a reference block (BK7 TM borosilicate glass).

The indenter tip is a Vickers diamond pyramid. IIT tests were carried out both on a pristine silica sample (further referred to as PS, studied for sake of comparison) and on FD. A typical '4-5-4' loading sequence was used: 4 s to reach the maximum force P m , 5 s of holding time, and 4 s to unload the sample's surface. They were force-controlled and the P m values were 100 mN. The mechanical response of the indentation test is the force P vs. the displacement δ (counted positively). Due to the high reproducibility of the indentation test on the glass surface, five indents per chosen maximum load were performed.

Unloading curves were analysed to extract an indentation modulus, M , and an indentation hardness, H IT , by Oliver and Pharr's method [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF].

Atomic force microscopy made it possible to record the three-dimensional geometry of the residual indentation imprint after unloading. Images were captured by using the contact mode of an AFM (Bruker Nanoscope V, Dimension 3100, Santa Barbara, CA) equipped with silicon nitride tips (DNP) with a nominal tip radius of the order of 20 nm and apical angle of 40°. The AFM

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

tip being much sharper than the Vickers indenter it ensures that the residual imprint shape is captured without AFM probe shape convolution effect. Finally, the piezo scanner was calibrated in X-Y-Z directions with a grid having a 10 µm pitch (X-Y directions) of 200 nm deep squared holes (Z direction).

Regarding image post treatment, the size of the scanned areas containing an indentation imprint was large enough so that a sufficient area, unaffected by the indentation process, surrounds the imprint and may be used as a reference surface for post treatment. This surface of reference was extracted from the raw image by using a disk-shaped mask centered on the imprint. It was assumed that the surface of the tested samples, far from the indentation imprint, was flat (i.e. no tilt and no offset). Thus a linear fit was applied to this data set, which was subsequently subtracted from the raw image. It is important to note that removing a tilt on the whole data set does not modify the three-dimensional shape of the residual imprint. Meyer's hardness, H, was also measured as the ratio P m to the projected contact area, after unloading, on the residual imprints.

Numerical procedures

Finite element analyses (FEA) of the indentation process were performed using both two (2D) and three dimensional (3D) models constituted of a sample and an indenter. The 3D mesh is generated from a coarser 2D mesh. The 2D mesh is first described. The sample's mesh is divided into a core zone, beneath the indenter tip, where the mesh is fine, and a shell zone where the mesh is coarse far from the contact. The core zone is itself divided into a square zone with a 32x32 (2D) or 16x16 (3D) quadrangular structured mesh contained into an outer un-

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

structured zone made of quadrangular elements (32/16 again along the axis z = 0). The shell zone is decomposed into a transition zone where the element size is progressively increased and a outer zone, both with quadrangular elements. All elements are linear and use full integration 1 . The dimensions of the mesh are chosen in order to minimise the effect of the far-field boundary conditions. This is made by using a sufficient number of outer elements in the shell zone. The typical ratio of the maximum contact radius and the sample size is about 2 × 10 3 . The 3D model relies on 45°swept 2D meshes since the Vickers pyramid has symmetries so that only 1/8 th of it has to be meshed. The resulting elements are eight-node hexahedrons except along the axis where linear six-node prismatic elements are used. Detailed views of the meshes are shown in Figure 1.

As for the indenter, the same procedure used for the sample mesh is employed. For 2D simulations, a triangle with an half-angle of 70.23°(equivalent angle of the Vickers pyramid) is meshed. For 3D models, a triangle with an half-angle of 68°is first meshed before being swept.

The indenter material (diamond) is assumed to be isotropic, linear elastic (Poisson's ratio ν i of 0.07 and Young's modulus E i of 1140 GPa [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF]). The contact between the indenter and the sample's surface is strict (Signorini conditions) and taken as frictionless. The contact zone and potential piling-up will take place along the square elements of the core zone. The boundary conditions consist of a null radial displacement along the vertical axis for both the sample and the indenter and a null displacement on the outer nodes of the sample. The force on the indenter, P (taken as positive), is controlled and the displacement of the indenter far from its tip, δ (counted positively), is 1 No evidence of possible volumetric locking beneath the indenter tip has been observed.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

recorded. The experimental maximum force was P m = 10 mN 2 . The problem is solved using the commercial software ABAQUS TM (version 6.15). The preand post-processing tasks were made with Abapy toolbox [24]. Details and views of the meshes are also visible in the documentation associated with Ref. [24].

The mechanical behaviour of the densified silica glass is assumed to be elasto-plastic (rate-independent) with a linear isotropic elasticity (Young's modulus E and Poisson's ratio ν). The plastic part, topic of this investigation, will be addressed hereinafter The finite strains framework is employed.

Results

Experimental results

The curves Pδ of both PS and FD specimens are shown in Figure 2.

FD is obviously much stiffer than PS. At a given load, the penetration depth is lower during the loading stage. Because of geometrical similarity of sharp indentation [START_REF] Tabor | The physical meaning of indentation and scratch hardness[END_REF][START_REF] Cheng | Scaling, dimensional analysis, and indentation measurements[END_REF], the loading part of the curves is a parabola whose curvature C is the indentation loading stiffness 3 . Values for both specimens are reported in Table 1, and it can be seen that FD is 44 % stiffer than PS. After unloading, the residual penetration depth, δ f , is consistently lower for FD vis-à-vis PS. The ratio of the reversible work released during unloading, W r , to the total work involved during loading, W t , is also indicative of the 2 Actually only P m /8 is imposed because of the symmetries reported hereinbefore, for the 3D mesh. 3 Actually, it is not a stiffness, and this why it is called the loading pre-factor usually. Indeed the stiffness is the slope of the Pδ curve, which is then C × δ, so a linear stiffness. Still, taking aside the kinematics δ, the parameter of interest describing the stiffness is C.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

indentation response. This ratio is 18% higher for FD than for PS. This might qualify it as "more elastic" than pristine silica.

The indentation modulus, M , is also reported in Table 2. The contribution of the indenter stiffness was subtracted and, for linear isotropic cases, M is E/(1 -ν 2 ). Elastic parameters of densified samples have been characterised in literature [START_REF] Zha | Acoustic velocities and refractive index of SiO2 glass to 57.5 GPa by Brillouin scattering[END_REF][START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF]. Taking the average reported value of 0.22 for ν and the M value from this work, E is found to be 106 GPa,which is in the range 110 ± 6 GPa reported in literature [START_REF] Deschamps | Elastic moduli of permanently densified silica glasses[END_REF]. These elastic parameters are reported in Table 2.

The residual imprints are shown in Figure 3 along with the profiles extracted along the four edges and the four faces. All imprints on FD were free of cracks, while edge cracks (quasi circles around the imprint, on the faces, see Ref. [29? ] for a detailed description of cracking systems in silica glass) were observed for PS. PS exhibits a sinking-in phenomenon, since all residual depths are below the reference level of the initial surface of the sample.

The profiles extracted (where the contribution of a edge crack is noticeable) confirm this. As for FD, the opposite situation is observed. An unusual long range piling-up is visible where the matter is above the initial reference level at a distance up to ∼ 3 µm from the imprint center. This is the case along the faces but also along the edges, which is rather singular for an inorganic glass. Usually, we can have piling-up along the faces on some silicate glasses [START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF] but not along the edges. For metallic glasses, piling-up may be even more pronounced but only along the faces [START_REF] Keryvin | Indentation and scratching mechanisms of a ZrCuAlNi bulk metallic glass[END_REF][START_REF] Keryvin | Hardness, toughness, brittleness and cracking systems of an iron-based bulk metallic glass by indentation[END_REF], not along the edges. This is therefore quite a rather unique situation for this fully densified silica glass.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Numerical results

Since the fully densified silica glass may not densify anymore, it sounds logical to assume first that the only mechanism for permanent deformation is volume-conservative plastic flow. The subsequent absence of any strain hardening mechanism in shear [START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF][START_REF] Molnár | Densification dependent yield criteria for sodium silicate glasses -An atomistic simulation approach[END_REF] made us use a von Mises material for describing the behaviour of FD, i.e. a rate-independent elasto-plastic model without strain hardening with the second invariant of the deviatoric part of the stress tensor (the von Mises equivalent stress) acting as the sole driving force. No contribution of the two other invariants is considered. In this simple model, the only plastic parameter is the yield strength σ y . We then performed a straightforward identification procedure by minimising the gap between the simulated force-displacement curve and the experimental one (from Figure 2) with this plastic parameter as the only free variable 4 . A remarkable agreement is found for σ y = 6.5 GPa, as observed in Figure 4.

Material parameters are presented in Table 2. Due to elastic stiffening, the compressive yield strain is strongly decreased, by 40 %, when comparing FD to PS.

Discussion

Molecular dynamic simulations (MD) have proven to be relevant tools for examining the multiaxial response of silica glass to mechanical testing, which is usually uneasy experimentally. Analysing the results in terms of driving forces (shear, pressure) versus mechanisms (densification, shear flow,. . . )

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

might give some insight on the choice of a given constitutive modelling a the continuum level, as well as the modifications in the short-to-medium range order. Yuan and Huang [START_REF] Yuan | Brittle to ductile transition in densified silica glass[END_REF] found that shear flow in densified silica can be attributed to a less rigid local environment associated with 5-fold Si coordination defects, more amenable to local shear deformation as compared with the 4-fold coordination state. Very recently, Molnar et al. [START_REF] Molnár | Densification dependent yield criteria for sodium silicate glasses -An atomistic simulation approach[END_REF] conducted MD simulations on binary NaO-SiO 2 glasses and observed that the shear yield strength was unaffected by the densification process. In both cases, the mechanical behaviour of fully densified glasses was found to be adequately described by pressure-dependent yield criteria such as Drucker-Prager (when pressure is an additional driving force) or extended Drucker-Prager (with the third stress invariant acting as another additional driving force). This situation ressembles that of some metallic glasses or polymeric glasses whose behaviour has been described by such models [START_REF] Donovan | Compressive deformation of amorphous Pd40Ni40P20[END_REF][START_REF] Patnaik | Spherical indentation response of metallic glasses[END_REF][START_REF] Keryvin | Indentation as a probe for pressure sensitivity of metallic glasses[END_REF]37]. These models involve a higher number of plastic parameters, in terms of yield criterion and plastic flow rule, than for the von Mises model we have used and for which only σ y is required. What is remarkable is that this latter simple model suffices to describe quite perfectly the mechanical response of the indentation test on FD. MD simulations (again very recent) in Ref.

[38], on pure silica glass, indicated that the shear yield strength was also unaffected during the densification process. Our present work as well as experimental results published in the literature [START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF] are very consistent with this statement. Indeed, experimental values reported in Table 2 show that the yield strength of PS and FD samples is the same.

Usually, the classical method for analysing the unloading stage of the force-displacement curve [START_REF] Oliver | An improved technique for determining hardness and elastic modulus using load and displacement sensing indentation experiments[END_REF] fails in extracting the correct indentation mod-

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

ulus (or Young's modulus provided Poisson's ratio is known) when piling-up occurs. Indeed, piling-up makes this method underestimate the contact area (∼ 30 % for a Zr-base metallic glass [START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF]) and therefore overestimates the indentation modulus. Yet our indentation results are very consistent with literature ones in terms of elastic parameters. We examine the reasons for this possible discrepancy by FEA. In Figure 5, are compared the numerical profiles along the edges and the faces to the experimental profiles. The comparison is again very satisfactory, excepted at the indentation imprint apex where a discrepancy is visible. The later may result from the finite size of the experimental indenter (tip roundness and truncated length) as AFM probe tip shape effect may be ruled out due to its much sharper geometry. Apart from this very localised zone, the profiles are well described and the residual piling-up along both the faces and the edges are correctly assessed. Thanks to FEA, we have also plotted the profile at maximum load, therefore before unloading. We can observe that a massive sinking-in develops during loading so that the residual piling-up is not a direct consequence of shear flow but rather of the elastic spring-back of the loaded surface during the unloading stage. This means that it is relevant to use the classical method for extracting the indentation modulus and this is why our values are in total agreement with literature results.

A more precise comparison between experimental data and numerical simulations is shown in Figure 6. It is a topography view of the residual imprint.

Apart from the red part (right), the comparison is very relevant. As for the red part, it shows that the numerical simulation underestimates the pile-up height, maybe due to a small local tilt of the sample since this red part is not present along the four faces of the imprint. We believe that this discrepancy can be ruled out by performing enhanced simulations with a much finer mesh

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

in the core zone, but this is not the topic of the present work.

To go a step further on this issue, we must now refer to results from a numerical benchmark proposed by Charleux et al. [START_REF] Charleux | A method for measuring the contact area in instrumented indentation testing by tip scanning probe microscopy imaging[END_REF] with two different classical constitutive equations (Hollomon and Drucker-Prager) to extend this observation. Using the same procedure, we can get some insight on possible piling-up and sinking-in both during loading and loading. This is made (in 2D) for the sole Hollomon model with a strain hardening coefficient n = 0 corresponding to the von Mises model, which was shown earlier in this work to be a relevant candidate to quantitatively describe the mechanical response of FD to the indentation test (cf. Figures 4 and5). Results are presented in Figure 7. It is observed, that up to a yield strain ε y of ∼ 1% piling-up occurs both during loading and unloading. After this value, sinking-in prevails during loading while piling-up is always observed after unloading, even if it decreases to very low values as ε y increases. The case for FD corresponds to an experimental value of ε y = 6.1%; it is superimposed and matches very well the numerical results.

The classical relationship between hardness and compressive yield strength [START_REF] Tabor | The physical meaning of indentation and scratch hardness[END_REF] is H ≈ c Y (Tabor's rule-of-thumb), where Y is the yield strength of the material for a given representative strain. For non strain hardening materials, it is σ y . For crystalline metals, c, the constraint factor, is taken as 3 and is frequently used as a convenient rule-of-thumb. Actually, this is strictly speaking between Vickers's hardness H v and σ y [START_REF] Tabor | The physical meaning of indentation and scratch hardness[END_REF]. There is a ∼ 7% difference between H and H v , so c ∼ 2.8. As reported in Table 2, this ratio is ∼ 1.25 for PS and 1.9 for FD. While it can be be treacherous to use Tabor's rule for a material like PS exhibiting an additional densification mechanism under indentation, this should be quite straightforward for FD since it is adequately modelled by

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

a von Mises material. Yet, for FD, the ratio hardness-to-yield strength is much lower than Tabor's rule. If this number of 2.8 has been revisited since the original work of Tabor, the recognised values nowadays range from 2.6 to 2.7 [39].

It has been shown that for pressure-sensitive glasses, like metallic glasses or polymeric glasses, this ratio can be higher [40, 41, 37]. However, for FD this ratio is much lower than the value for pressure-insensitive materials like metals. Using FEA (2D), we have plotted in Figure 8, the ratio hardness-to-yield strength (H/σ y ) for von Mises materials (for ν = 0.22) as a function of the yield strain ε y in the particular situation where the indenter is of Vickers type.

For this indenter, the so-called indentation strain ε i = cotan θ , where θ is the half angle of the equivalent cone (here 70.27°), is 0.358. Classically we have two regimes described for example in Ref. [START_REF] Cheng | Scaling, dimensional analysis, and indentation measurements[END_REF]: a fully plastic regime, for low values of ε y , corresponding to a plateau value of ∼ 2.65, the so-called constraint factor c; and an elasto-plastic regime where H/σ y is lower than c.

The experimental case of ε y = 6.1% for FD is superimposed on Fig. 8, and matches very well with the numerical results. This indicates, that when using a Vickers indenter, we are still in the elasto-plastic regime of indentation for FD, explaining the low value of the ratio hardness-to-yield strength, which is not in this case the constraint factor defined for the sole fully plastic regime of indentation. The Vickers indenter is not sharp enough to reach this latter regime.

Concluding remarks

We have characterised the mechanical properties (elastic stiffness, hardness) of a fully densified silica glass samples (21.6 % increase in density) by using instrumented indentation testing. We have also shown that the ma-

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

terial's response to Vickers indentation (force-displacement curve, residual imprint) is correctly assessed by a constitutive equation usually relevant for some crystalline metals, i.e. a rate-independent pressure insensitive elastoplastic model without strain-hardening, usually referred to as a von Mises material. This singular comparison with metallic alloys ends rapidly since the yield strength and yield strain values (respectively 6.5 GPa and 6.1 %) are far higher than those for crystalline alloys. These values are in agreement with very recent experimental and simulation studies on densified silica or silicate glasses. Apparent disagreement with indentation features (a perfect assessment by the Oliver and Pharr's method of Young's modulus in the presence of a massive piling-up around the residual imprints ; a ratio hardness-to-yield strength of 1.9 much lower than the constraint factor of ∼ 2.7 for crystalline alloys) have been explained again by the very high value of the yield strain. Indeed, finite element analyses have shown that, during the loading stage of indentation, a massive sinking-in occurs and it is because of the very high elastic spring-back that unusual (as compared to classical alloys or metallic glasses) piling-up is observed around the residual imprints. The same analyses have highlighted that the material is in the so-called elastoplastic indentation regime and not in the fully plastic one, for this Vickers indenter, where a hardness-to-yield strength value of ∼ 2.7 is expected. This mechanical description of fully densified silica glass mechanical behaviour under Vickers indentation, to our opinion, paves the way for new insights into the constitutive modelling of pristine silica glass during indentation loading, which is a key to understand the transition from plasticity to cracking in these oxide glasses to increase their use in service lifetime.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT M A N U S C R I P T A C C E P T E D M A N U S C R I P T A C C E P T E D M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT C [GPa] W r /W t [%] δ f [nm] M [GPa] H [GPa] H IT [GPa]

Figure 1 .

 1 Figure 1. View of the 2D and 3D meshes used in the FEM simulations for the material close to the contact location. The core and shell zones are represented respectively in red and blue. The 3D mesh is a 45 o swept version of a coarse 2D mesh. l c stands for the horizontal length of the core zone.

Figure 2 .

 2 Figure 2. Nano-indentation force-displacements curves for pristine silica glass and fully densified silica glass. Five tests are carried out for each type, and experimental reproducibility is confirmed.

Figure 3 .

 3 Figure 3. Topography AFM images (5×5 µm 2 , top) presenting the residual imprints of aVickers indentation (100 mN) for the pristine silica glass (left) and the fully densified silica (right). For sake of presentation the altitudes below -10 nm are masked (white part). For pristine silica (used only for comparison), edge cracks are visible. For fully densified silica, colors different from blue are indicative of a residual pile-up (z > 0). A unusual long range residual pile-up all around the imprint is observed. Altitude profiles are then extracted (bottom) along the four edges (green) and the four faces (magenta) for both specimens. No piling-up is visible along the faces and the edges for pristine silica. Piling-up is present along the faces and the edges for fully densified silica silica.

Figure 4 .

 4 Figure 4. Force-displacement curves of the nano-indentation test in fully densified silica at 100 mN. Comparison between experimental results and numerical simulation (3D FEA).

Figure 5 .

 5 Figure 5. Comparison between experimental (AFM) results and numerical simulations (3D FEA) of the residual imprint profile in fully densified silica indented at 100 mN by a Vickers pyramid. The four experimental profiles are plotted together for sake of presentation either along the edges (right) or the faces (left). The profiles at maximum load, as predicted by the simulations (blue color), indicate that there is no piling-up during loading.

Table 1 .

 1 Mechanical properties of pristine silica and fully densified silica from instrumented indentation testing. C is the loading pre-factor. W r is the reversible work during unloading and W t the total work involved in the test. δ f is the residual penetration depth, and M is the indentation modulus of the material. H is Meyer's hardness measured from the residual imprint. H IT is the IIT hardness following Oliver and Pharr's method[42].

Pristine 106 ± 1 31.7 ± 0.1 486 ± 1 69 ± 1 8.40 ± 0.30 8.55 ± 0.05 Fully densified 153 ± 1 37.3 ± 1.2 412 ± 3 112 ± 2 11.5 ± 0.5 11.5 ± 1.0

We must mention that this procedure suffers from non uniqueness when two or more parameters are to be identified. This is not our case where only one parameter is determined, as pointed out in [26, Fig.7].

Acknowledgements

We would like to thank H. Ji (Univ. Rennes 1, France), T. Hammouda (Univ. Blaise Pascal, France) for experimental help and T. Rouxel (Univ. Rennes 1, France) for stimulating discussion.

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT

Condens. Matter 20 (2008) 114119.

[37] K. E. Prasad, V. Keryvin, U. Ramamurty, Pressure sensitive flow and constraint factor in amorphous materials below glass transition, J. Mater. Res. 24 (2009) 890-897. [40] V. Keryvin, Indentation of bulk metallic glasses: Relationships between shear bands observed around the prints and hardness, Acta Mater. 55 (2007) 2565-2578.

[41] V. Keryvin [42] W. Oliver, G. Pharr, Measurement of hardness and elastic modulus by instrumented indentation: Advances in understanding and refinements to methodology, J. Mater.

Res. 19 (2004)

Pristine (from [START_REF] Lacroix | Plastic deformation and residual stresses in amorphous silica pillars under uniaxial loading[END_REF][START_REF] Kermouche | Perfectly plastic flow in silica glass[END_REF]) 72 0.15 ∼ 6.5-7.0 ∼ 9.5 1.2-1.3

Fully densified (this work) 106 ± 1 0.22 6.5 6.1 1.89 ± 0.08

Position, (x, y)

[µm] . Normalised indentation stress (H is Meyer's hardness and σ y is the material's yield strength) versus normalised indentation strain (cotan θ , where θ is the half angle of the equivalent cone, here 70.27°; ε y is the material's yield strain). Numerical simulations (2D, ν = 0.22) allow to distinguish between the fully plastic regime (FP), for which H/σ y is a constant c = 2.65 and the elasto-plastic regime (EP), for which H/σ y is an increasing function of 1/ε y . The experimental value for fully densified silica (FD) is superimposed.


Sylvain Chatel
email: sylvain.chatel@supelec.fr

Pierre Haessig

Romain Bourdais
email: romain.bourdais@centralesupelec.fr

Security of the Distributed Model Predictive Control

In this paper, we explore the security of a distributed model predictive control (DMPC) scheme. A global system is decomposed into several local subsystems which share a common resource. For instance a heating problem among several rooms with a limited amount of power. Each subsystem owns a model predictive control (MPC) controller which cooperates and works iteratively with all the other controllers in order to reach the system's objective. Here, we focus on the security of such a system. Indeed, we report that a user could disturb the DMPC through several attacks in order to destabilize the resource distribution. Hence, examples of cheating protocols and counter-cheat ideas are presented and discussed.

I. INTRODUCTION

Model predictive control (MPC), also known as receding horizon control, is becoming more and more popular in the industrial process control [START_REF] Venkat | Distributed mpc strategies with application to power system automatic generation control[END_REF], [START_REF] Camponogara | Distributed model predictive control[END_REF]. Among other benefits, the attractiveness of the MPC lies in its ability to represent clearly the constraints of the optimization problem [START_REF] Jia | Distributed model predictive control[END_REF]. Over the past decades, the MPC has become widespread and many successful applications have been developed in the process industry.

In an MPC, the objective is to get the optimal control input over a given horizon by solving a discrete-time optimal control problem. Traditionally, this kind of control is done by a controller in a centralized way. This controller have full knowledge of the system over the horizon. Although efficient, this schemes might sometimes be inadequate. For instance, for large-scale interconnected systems such as power and water distribution or even traffic systems (with plug and play), a centralized MPC might not be ideal or even technically feasible [START_REF] Jia | Distributed model predictive control[END_REF]. Moreover, the fact that the controller is omniscient and fully knowledgeable might be an obstacle because some agents might not be able to divulge information about their local subsystem. This might be the case as [START_REF] Camponogara | Distributed model predictive control[END_REF] points out for the newly regulated power markets in the United States.

Hence the idea of developing a distributed system. The global system is decomposed into several subsystems. Each subsystem has a MPC controller which communicates with the others. Local control inputs are then locally computed and determined with the few information shared. For instance, when a decentralized control scheme (where local control input are determined based on local measurements for instance) is needed, a distributed MPC is a better choice. For example, for water distribution systems, a centralized MPC is decomposed into a DMPC using a coordinator (e.g. augmented Lagrangian [START_REF] Camponogara | Distributed model predictive control[END_REF] or Uzawa method [START_REF] Cohen | Optimisation des grands systèmes[END_REF] depending on the type of dependence between the subsystems).

Several distributed MPC examples are available in the literature. DMPC frameworks were developed in [START_REF] Acar | Some examples for the decentralized receding horizon control[END_REF], [START_REF] Sawadogo | Decentralized predictive controller for delivery canals[END_REF] and [START_REF] Dunbar | Distributed receding horizon control for multi-vehicle formation stabilization[END_REF]. In the latest, a MPC framework is developed for a system with several independent subsystem dynamics but nonetheless linked through their cost functions.

In a DMPC, each agent communicates with the others via a simple coordinator. This latest, returns a simple information (e.g. a Lagrangian multiplier for an augmented Lagrangian technique). All the information required for the system to work efficiently is communicated through this small amount of data. As pointed out in [START_REF] Brooks | Demand dispatch[END_REF], though efficient the DMPC still suffers from some risks. Indeed, DMPC might sometimes be harmful to the physical system. For instance in a complex system such as a power grid, human error, malicious and misleading agents could easily disturb the grid. According to Brooks in [START_REF] Brooks | Demand dispatch[END_REF], one of the best ways to mitigate those security threats is to involve a combination of counter measures such as using fault detection softwares, limiting the size of each entity or even certifying aggregators. However, to the best of our knowledge, there has been no comprehensive study on security threats and counter measures at the level of DMPC algorithms (as opposed to the lower level layers like communication protocols). This article is a first attempt to bridge this gap.

The paper is organised as follow. Section II introduces the model of DMPC for power distribution we used in our study. In Section III, a study of the base cases is conducted. The security threats are implemented and studied in section IV. In Section V, the main contribution of this study are summarized.

II. MODEL

In this paper, we consider a linear time-invariant (LTI) system composed of m interconnected subsystems. Each subsystem represents a room with a thermal resistance R thi and a thermal conductivity C thi . Taking into account the exterior temperature T ext , the objective is for the temperature of each room T i to be as close as possible to the reference temperature T i id at each time. All rooms are getting their power from a limited global amount of power U max , and each user has a maximum admissible power u i max due to its physical configuration. We note u i the power consumed by user i. Since the consumed power cannot exceed the maximum admissible power, and since the global power resource is limited, we have the first two constraints :

∀i ∈ 0, m -1 , 0 ≤ u i ≤ u i max (1) m-1 i=0 u i ≤ U max (2)
The main objective is to minimize the consumed energy represented by the sum. In order to take into account the will of each user to reach its ideal temperature, we introduced a comfort factor α. Hence the optimization problem is to minimize the following cost function subjected to 1 and 2.

J u = m-1 i=0 u i + m-1 i=0 α i (T i -T i id ) 2 (3)
In order to determine the J u , we used the thermal model presented in Figure 1.

With this model, it is possible to use the state-space representation of the problem. We can set X[t k ] as the temperature vector of each room. U [t k ] represents the consumed energy of each room and U ext [t k ] the external temperature. With those notations, we can write :

X[t k+1 ] = A.X[t k ] + B.U [t k ] + B .U ext [t k ] (4)
Given equation ( 4), and the model, we can replace the temperature in equation ( 3). This enables us to have a quadratic problem in U :

J u = U T P U + q T U + cst (5)
In order to proceed to a centralized MPC, we just have to minimize 5 subjected to 1 and 2. However, in order to have a DMPC, we need to decompose the computation to all users. To this end, we used an Uzawa method. The idea is to unconstrain the optimization problem and relax the constraint in 2. Then, after solving the optimization problem, we iterate the Lagrangian multiplier λ.

λ k+1 = λ k + p * ( m i=0 u * i [t k ] -U max ) (6)
In 6, p represents the step of the Uzawa iteration and u * i [t k ] the optimal consumption of user i at t k . This iterative procedure is done until the difference between

m-1 i=0 u * i [t k ] FIGURE 2. Structure of the distributed problem
and U max is lower than a threshold. We note λ opt the value of the multiplier fulfilling this condition.

To summarize, the DMPC requires to solve the following quadratic problem over a prediction horizon N for a simulation lasting over N sim :

   min J u [t k ] = U [t k ] T P U [t k ] + q T U [t k ] + cst[t k ] + λ[k opt ] subjected to ∀t k ∈ 0, N sim -1 , ∀i ∈ 0, m -1 , 0 ≤ u i ≤ u i max (7)
In our study, we simulate 24 hours, which is N sim = 240 points with a time step ∆t = 0.1h.

III. BASE CASES

In this section, we present the results for the nominal experiment of power distribution in several cases : the static case, the centralized dynamic case and the distributed dynamic case. All those simulations were made via the python module we developed. This package relies on M. Andersen, J.Dahl and L. Vandenberghe cvxopt package for convex optimization [START_REF] Boyd | Convex Optimization[END_REF].

A. Static case 1) Centralized case: Let us consider a problem with three different users. We consider, at first, that all users are identical. We set the system so that U max = 3kW is inferior to the sum of all admissible temperature. Power allocation = 2.50000/2.5

u * U max /n <u * > u max u id FIGURE 3. Optimal power distribution in a static symmetrical case
Now, if we change the comfort factor of user 2 in order for him to be more comfortable. To do so, let us take α 2 = 10 × α 0 . When comparing Figure 3 and 4, we notice that the third user (n • 2 is indeed preferred in the power distribution and then its comfort becomes better to the disadvantage of the others (i.e. the square deviation to the ideal temperature of the third user (n • 2 is the smallest of all users).

2) Distributed case: When we decompose the centralized computation with the Uzawa method, we are able to model a distributed optimization. At first we consider a symmetrical situation. Now, if we change the comfort factor of user 2 in order for him to more comfortable. To do so, let us take α 2 = 10 × α 0 . When comparing Figure 5 and 6, we also notice that user 2 is more satisfied to the disadvantage of user 0 and 1. We notice as well that the total amount of consumed energy is superior to U max . This phenomenon is due to the constraint relaxation executed in the Uzawa method.

Through this nominal static study, we were able to determine the influence of all the different parameters on the optimal solution among which the Uzawa step determination to have a convergence of the Uzawa iteration. We then used those results in the study and set p = 1.5.

B. Dynamic case

Now, we consider a dynamic simulation on a horizon of 24 hours. Three different users are using a limited amount of energy U max = 3kW . The ideal temperature is determined with a time based profile : between 6:30-8:30 am and 18:00-22:00 pm the user wishes to have a temperature T pres = 22 • C (people are present in the room), the rest of the time the temperature is set to T abs = 18 • . As before, we now consider a situation where user 2 is preferred and has a comfort factor ten times superior to others. FIGURE 9. Optimal power distribution in a dynamic case for user 2 when it is preferred, α 2 /α 0 = 10, focus between 6:30-8:30 am We notice comparing Figures 8 and9 that in the second case, user 2 is given mush more energy when needed and hence can reach its ideal temperature much faster.

From this point on, we are able to use this model to begin the study of cheating and misleading scenarios.

IV. SECURITY ISSUES : MISLEADING AND CHEATING

USERS

In this section, we will focus on the security issues of the DMPC. First, we will analyse the different ways of cheating before putting it into application on static and then dynamic situations.

A. Security concepts

First of all, we must clarify the concept of security breach in a DMPC. Through our research, we were able to detect two kinds of threatening users : greedy users (whose objective is to maximize their comfort) and nihilist users (whose objective is to destroy the system).

In Figure 10, we can see an overview of the different threat we detected. The reader can note that even if other strategies might exist, we realized that all the other strategies can be somehow reduced to one of those we studied.

Moreover, for our study we considered that the transmission canal is ideal and that the communication cannot be altered. Indeed with the current knowledge, security protocols enable us to safely communicate and to change ones identity. This implies that in our model cheating by usurping ones identity or by altering the communication is not considered. Moreover, unless the demand profile of the users are complementary, the strategy of regrouping is ineffective.

To summarize, the base attack cases are as follow :

• Change its own comfort factor α i .

• Do not listen to the DMPC talks ie the Uzawa iterations results to satisfy the global power limitation. • Listen partially to the DMPC talks. This means the user only takes a fraction of the multiplier λ. • Duplicate itself. In other words, a user broadcast to the others that he is two or more users in order to get more power.

B. Security in static situations

Now, we implement those strategies into the static model. First, we implemented the influence of the comfort factor α. The idea was to determine to which extent this parameters changes the power distribution. In Figure 11, we represented the optimal power distribution and the temperature deviation for two users when the comfort factor is a parameter. We noticed that, as expected, when the one user comfort factor is superior to the one of the others, it receives much more energy and as therefore a better comfort. Then, we determine the influence of the ideal temperature on the distribution. In Figure 13 is represented the power distribution and the temperature deviation as functions of the broadcasted temperature. We notice that the higher this temperature is, the more energy the user receives and the less its deviation is. Please note that in our model if a user receives energy, this user has to consume it (hence the positive deviation in Figure 13). Finally, we created a situation where a user could listen or not to the talks of the DMPC. This means the user can decide if it takes into account the Lagrangian multiplier returned by the coordinator. We introduced this notion as the deafness β. If the user behaves and listen, β = 0, but if he is deaf, and does not take the talk into account, β = 1.

   min J u [t k ] = U [t k ] T P U [t k ] + q T U [t k ] + cst[t k ] + (1 -β)λ[k opt ] subjected to ∀t k ∈ 0, N sim -1 , ∀i ∈ 0, m -1 , 0 ≤ u i ≤ u i max (8)
In Figure 14, the nominal situation is presented and in Figure 15 a situation were user 1 is deaf is presented. We do notice that when the user is deaf, its comfort is much better than when it listens to the talk. Of course this deafness factor β can be set between 0 and 1. Figure 16 shows the influence of the deafness factor of user 1 on the temperature deviation of both users. With no surprise, when a user is not listening to the talks all the others are penalized.

C. Security in dynamic situations

Now, let us consider the dynamic model. Our developed python module is able to deal with centralized MPC and distributed MPC. However, due to the complexity of the DMPC, its calculation requires more computation power. So for the sake of the simplicity the results of this part were obtained with a classic MPC but they can all be extended to the DMPC.

1) First, we consider a cheating user: We want to observe the influence of the comfort factor. In Figures 17 and18 we can observe the situation with two users. In blue we can see the the distribution in the nominal case were α 1 = α 2 = 10. Then we change the value of the comfort factor for the second user. We notice that its comfort is drastically better (follows the reference temperature more than the green and the blue lines) to the detriment of the first user (in green). FIGURE 18. Power distribution through time for two users. In blue is the nominal situation with identical comfort factor. Focus between 6:00-7:00 am Our study pointed out that the more α 2 > α 1 , the more user 2 has a great comfort and the opposite for user 1. Hence, by modifying α, one could easily cheat the DMPC in order to gain access to more energy.

2) Then, we consider a nihilist user.: As previously stated, a nihilist user is a user whose objective is to disrupt the system and even destroy it. In our first subsection, we presented several strategies in order to disable the DMPC system. During this phase, we evaluated the influence of different parameters to detect if they could be a liability in the DMPC. We notice that by augmenting the comfort factor of all users, we could reach a threshold and make the Uzawa iteration not converging (or at least reaching the maximum admissible iteration k max ). Hence if a user succeeds to change the value of all the α, it might be able to break the DMPC.

Secondly, we tried the other method : to use a particular sequence of u * i [t k ] and not the result of the optimization.

We created a sequence such as u 2 [t k ] = u i max if k is an odd number and u 2 [t k ] = 0 else. In Figure 19, we show the distribution through time for the non-nihilist user (user 1). In Figure 20, we show the same but for the nihilist user (user 2). When summing u * 1 and u * 2 we notice that the max is superior to the amount of energy available U max . This proves that this sequence disable the DMPC.

D. Protection

Finally, we studied security counter-measures against those threat. We noticed through our experiment that when we were altering the parameters, we modified as well the following functions :

•

k U zawa → (u * (k U zawa )| α ) • λ → (u * (λ)| α ) • k U zawa → (λ(k U zawa )| α )
In other words, the first (resp. second) function represents the power distribution (of one user for a given α) as a function of the number of iteration realized in the Uzawa method (resp. the value of the Lagrangian multiplier). Similarly, the third one returns the value of the Lagrangian multiplier as a function of the number of iteration realized in the Uzawa method.

In Figure 21, we can observe the nominal situation. Then in Figure 22 we can see a cluster of lines created for different values of α.

We then decided to realize a large cluster of lines for alpha linearly distributed. The result is presented in Figure 22 for α in range 0 to 150.

Thanks to Figure 22, we can set a threshold (for instance on the line α 12 ), this will divide the space in two : above the line (u * (k U zawa )| α12 ), user will be considered as cheaters, else the users are cleared. Respectively below the line (u * (λ)| α12 ), user will be considered as cheaters, else the users are cleared. This is a way to detect efficiently cheating and nihilist users using a comfort variation approach.

V. CONCLUSION

In this paper, we presented new results on security breach and measures for distributed model predictive control. By developing a thermal model we were able to experiment various approaches on, firstly, how to destroy or cheat, and secondly on how to detect such behaviours. By analysing all the different method, we managed to reduce those situation to a few base cases much easier to analyse and export to other situations. All those results are encouraging and call for further investigation on the subject.

FIGURE 1 .

 1 FIGURE 1. Thermal model of a room. R th i represents the thermal resistance and C th i the thermal conductivity of the room. U is the consumed power, T the temperature of the room and Text the external temperature.

FIGURE 5 .

 5 FIGURE 5. Optimal power distribution in a static symmetrical distributed case

FIGURE 6 .

 6 FIGURE 6. Optimal power distribution in a static distributed case with different priorities, α 2 /α 0 = 10

FIGURE 7 .FIGURE 8 .

 78 FIGURE 7. Optimal power distribution in a dynamic case for user 2 in a symmetrical problem

FIGURE 11 .FIGURE 10 .

 1110 FIGURE 11. Parametric study of the comfort factor for two users

FIGURE 13 .

 13 FIGURE 13. Parametric study of the broadcasted ideal temperature when different from the real ideal temperature for three users

FIGURE 14 .FIGURE 15 .

 1415 FIGURE 14. Power distribution for three user in the nominal case

FIGURE 16 .

 16 FIGURE 16. Parametric study of the deafness factor on the deviation to ideal temperature

  17

FIGURE 17 .

 17 FIGURE 17. Power distribution through time for two users. In blue is the nominal situation with identical comfort factor.

  17

  17

FIGURE 19 .FIGURE 20 .

 1920 FIGURE 19. Power distribution through time for user 1 when applying a destroying sequence of consumed power for user 1

FIGURE 22 .

 22 FIGURE 21. u * (k), u * (λ) and λ(k) in nominal situation