R J Ball

G C Allen

The measurement of water transport in porous materials using impedance spectroscopy

Keywords: mortar, dewatering, impedance spectroscopy ______________________________________________________________________ 1

This paper describes the application of electrical measurements to monitor the extraction (movement of water from the mortar) of water from calcium lime, natural hydraulic lime and Portland cement mortars placed on an adsorbent brick substrate.

Impedance measurements were used to identify the changes in bulk resistance of the mortar. A model has been developed combining sharp front theory and Boltzmann's distribution law of statistical thermodynamics to identify the point at which no further absorption of water in to the brick occurs. A linear relationship was found between the exponential of bulk resistance and square root of time during dewatering. A change in gradient was attributed to the end of dewatering.

Introduction

In their pioneering publication [START_REF] Hall | Water Transport in Brick, Stone and Concrete[END_REF] observed that "If there is a science of building the interplay of water and materials must be one of its themes". The water level in a normal stone wall is determined by the rate of absorption and evaporation which in turn depends on the porosity of the stone. Thus the movement of water through building materials, especially brick structures, may be a limiting factor not only in their durability but also in process of their construction. Climatic changes in the environment will lead to buildings experiencing periods of wetting and drying where the measurement and calculation of the transport of vapour and liquid water is essential for the understanding of their deterioration and restoration. This paper describes the application of electrical measurements to monitor the dewatering and drying of the mortar materials calcium lime, natural hydraulic lime (NHL) [START_REF] Allen | Hydraulic Lime Mortar for Stone, Brick and Block Masonry[END_REF] and Portland cement [START_REF] Bye | Portland Cement[END_REF]) mortars used in the bonding of absorbent brick substrates. Previously impedance spectroscopy has been applied to cement mortars to study rapid hydration reactions (McCarter et al 2000) but lime mortars stiffen in an entirely different way where the main reaction involves the absorption of CO 2 following hydration (Van Balen 2005). For this reason dewatering becomes especially important. Carbonation is encouraged in moist rather than saturated materials and is governed by the relative humidity (El-Turki et al 2007, El-Turki et al 2009) In this study a model has been developed combining sharp front theory (Collier et al 2007) and Boltzmann's distribution law of statistical thermodynamics to estimate the rate of dewatering from the recorded changes in bulk resistance. In the experiments a mortar of suitable thickness was chosen to allow the rate of desorption to be carefully studied in detail. This allowed the initial rapid dewatering stage and the subsequent slower rate of diffusion to be assessed when the remaining concentration of water had fallen to a level determined by the capilliary action of the substrate.

Experimental Method

Sample preparation

Single rather than hybrid mortars were manufactured from a pure calcium lime (CL90), natural hydraulic lime (NHL3.5) and Portland cement individually. These materials were supplied by Hanson Cement Ltd, Clitheroe, Lancashire, UK. The oxide composition, determined by X-ray fluorescence spectroscopy (XRF), of the calcium lime (CL90) and Portland Cement was reported by A El-Turki et al (2009) and that of the NHL3.5 by [START_REF] Ball | [END_REF]. The bulk densities of the CL90, NHL3.5 and OPC binders used were 650.6, 666.3 and 1119.3 kgm -3 respectively. The manufacturers recommended volume ratio for binder to aggregate of 0.5 and water to binder of 0.78 was used in all cases. A quarried Croxden sand composed of mainly silicon with traces of aluminium, potassium, titanium, iron and calcium of bulk density 1475.5 kgm -3 was used as aggregate . It's particle size distribution and energy dispersive X-ray analysis was reported by [START_REF] Ball | [END_REF]. The dry constituents of each mortar type were intimately mixed before the water was added, prior to further mixing. Due to the small quantities of mortar required for each experiment all mixing was done by hand for 5 minutes following the addition of water. After mixing the mortar was immediately applied to the brick and the experiment initiated within 1 minute.

Substrate sorptivity

Mortars were set upon square section brick prisms of width 25 mm and length 80 mm, sufficient to fully dewater the depth of mortar applied. Prisms were cut from golden purple bricks supplied by Ibstock Brick Ltd, Roughdales Factory. The sorptivity of each prism was determined using the method described by [START_REF] Hall | [END_REF]. Prior to measurement prisms were oven dried at 105 °C for at least 12 hours. The weight of each piece of brick was measured, to an accuracy of 0.1 g, at time intervals of 1, 2, 4, 9, 16 and 25 minutes when the bottom face was in contact with a shallow layer of water.

Brick pieces were dried again in the oven under the same conditions before the mortar was applied. The volume of water absorbed by the brick at each of the time intervals was calculated from the density of water, 1000 kgm -3 , and the weight change. Water volume was then plotted against the square root of time and a linear trend line fitted using a least squares algorithm. By dividing the gradient of the line by the cross sectional area of the prism a value of the brick desorptivity, R, was calculated.

Impedance Spectroscopy

The application of impedance spectroscopy to cementitious systems is well documented [START_REF] Christensen | [END_REF]. The impedance response of the wet mix was monitored using a Solartron 1260 impedance analyser over the frequency range from 10 MHz -100 Hz.

Each sweep between these frequencies contained 5 steps per decade and were carried out at a potential of 100 mV. 200 sweeps each lasting 70 seconds were repeated continuously over a time period of 225 minutes. The test cell shown in figure 1 (a) was constructed using two rectangular stainless steel electrodes 25 mm by 30 mm spaced 25 mm apart. These were positioned in a square section tube of the same dimensions as the brick prism to allow unidirectional flow of liquid between the wet mix and substrate.

The electrodes were positioned 15 mm above the brick surface.

Theory

The dewatering of mortar in contact with an absorbent substrate can be analysed using the sharp front model developed by [START_REF] Hall | Water Transport in Brick, Stone and Concrete[END_REF]. This model describes the relationship between the sorptivity, S, of the substrate, the desorptivity, R, of the slurry and the transfer sorptivity, A, between the slurry and the substrate. The withdrawal of water from a wet mix by an initially dry substrate occurs by capillary forces. The slurry can be divided into a 'filter cake' and 'wet mix' as shown in figure 1 (b). As dewatering progresses the thickness of filter cake increases until no wet mix remains or the brick substrate becomes completely saturated. The application of Darcy's law to the filter cake yields equation 1.

i K L K dt di i c c i c Ψ - = Ψ - = β α 1
Where i (m) is the cumulative volume of water desorbed per unit area of wet mix, t (s) is time, K c (ms -1 ) is the permeability of the cake, Ψi (MPa) is the capillary potential, L c (m) is the depth of the cake, and, α and β, are terms relating to the volume fraction of water. i can be defined by equation 2 in terms of the cumulative volume of water adsorbed, V H 2 0 (m 3 ) and the area of the wet mix in contact with the substrate,

A sub (m 2 ) . sub O H A V i 2 = 2
Integration of equation 1 gives equation 3.

2 1 2 1 2 1 2 t R t K i wm i c = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Ψ = β α 3
Alternatively if the area of the brick substrate, A sub , is represented by a capillary of radius, r, the flow of water can be described using the Washburn equation, equation 4 (Washburn 1921).

η ρ η θ γ 8 4 cos 2 g r h r dt dh v - = = 4
Where h is the depth of water adsorbed into the brick, γ cos θ is the surface tension and η is viscosity. Integration of this equation ignoring the gravity term yields equation 5. shown in equation 6. However in this case the flow is described by the properties of the brick substrate and not the wet mix.

2 1 2 1 4 cos t r h ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = η θ γ 6
If the initial volumetric water content of the mix at t = 0 is w 0 then the water content at time t, w, is given by equation 7. ( )

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - = kT E E N N j i j i exp 8
Where N i and N j are the numbers of constituent ions in states i and j and E i and E j are the respective potential energies, k is the universal constant equal to 1.38041 x 10 -23 JK -1 ;

and T is temperature, K.

When a solid with a charged surface is in contact with a liquid, containing ions, a double layer is formed. This comprises a tightly bound layer of fairly immobile ions adjacent to the solid and an oppositely charged ionic atmosphere (Dukhin and Dejaguin 1974). If stage i corresponds to ions within the double layer and stage j to those in the bulk water, then the difference in energy between the two states can be expressed in terms of the electrical Maxwell forces defined by equation 9.

( )

j i j i ze E E ϕ ϕ - = - 9
where ϕ is the electrical potential; e is the electronic charge; z is the summary valence of the mobile ions in the solution (Neilsen et al 1972). By combining equations 8 and 9

the fraction N i /N j representing the relative number of mobile electrical charges in the bulk water and double layer is obtained, equation 10:

( ) ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ - - = kT ze N N j i j i ϕ ϕ exp 10
If the number of ions in the double layer, N i , is assumed to be constant and those in the bulk water, N j , is proportional to the water content the ratio, N i /N j , is proportional to water content. The electrical potential difference φ i -φ j is proportional to the mortar bulk resistance according to Ohm's law. Thus, from equation 7 a relationship between mortar pore water content, w, and electrical resistance can be obtained, equation 11, where a ( )

2 1 0 exp t R A w bR a wm sub b - = - 12
Equation 12 can be re-arranged into the form, equation 13.

( )

a w t a R A bR wm sub b 0 2 1 exp + - = - 13
This equation predicts that if b, A sub , R wm , a and w o are constants a plot of exponential of (-b R b ) versus square root of time will be linear.

Results

Values of brick sorptivity for each of the prisms used are shown in table 1, ranging between 2.2 and 2.6 mm min 0.5 . Such variations are typical for the type of brick used and close enough to allow meaningful comparison between the mortars tested.

A typical complex plane plot for the CL90 mortar after a period of 30 minutes is shown in figure 2 Figure 3 shows plots of exponential of (-bR b ) versus square root of time for the calcium lime, natural hydraulic lime and Portland cement mortars respectively. A value of -0.001 J -1 A 2 s was used for b. A linear trend line was fitted to the beginning and end sections of data from each specimen. In each case a correlation coefficient of 0.99 was obtained.

The gradient of the initial section of the curve corresponding to dewatering is given in table 1. An increase in gradient was observed with mortar type from calcium lime, natural hydraulic lime and Portland cement mortars respectively.

Discussion

The mechanisms by which water is retained in mortars, plasters and renders in the wet mix will depend on their water content. This may take the form of free flowing water at complete saturation and physically bound water when surface tension effects become dominant. This approach has been reported previously by Zacharopoulou (2006) for lime putties.

The results presented here are those for mortars containing calcium lime, natural hydraulic lime and Portland cement binders. Figure 3 This may be expressed by sharp front theory to derive an equation for the time taken to dewater a wet mortar by an absorbent substrate. (C Ince et al 2009). As with Ficks law the rate at which the mix loses water depends on its thickness and assumes no restriction between the media involved in the transfer process giving a clear delineation at its endpoint.

Unfortunately though the inhomogeneous nature of the mortar brought about by the imperfect distribution of water, aggregate and binder and the variable composition of the adsorbate brick greatly complicates the dewatering process producing a less well defined endpoint. In the present study it is assumed that the change in gradient of the resistivity -time curve, shown in Figure 4, signifies a change in the dewatering process.

The extent of this inhomogenity is reflected by the time over which the change in gradient (transition) occurs.

When the slopes of the lines and the time at which the change in slope occurs are compared it is apparent that the type of binder used in the mortar has an important influence on behaviour. It is noteworthy that the steepest gradient was observed for the Portland cement mortar which contained the highest proportion of reactive hydraulic phases and where the chemical activity was greatest. Especially noteworthy was that in this measurement there was an indistinct change in gradient and hence little evidence of a two stage process. Moreover the conductance measurements reported by McCarter and Garvin (1989) also showed little change in gradient with changing water content in cement phases. Here dewatering appeared to be essentially continuous but the period of this process was longer than that for the natural hydraulic and calcium lime. This observation may be due to the fact that with Portland cement the rapid hydraulic reaction forms a silicate network through which the remaining free water must be progressively extracted with increased tortuosity.

In the case of CL90 the main setting reaction, carbonation, does not occur to any Knowledge of the dewatering time is of importance to the construction industry. This study provides an important method of measuring dewatering and subsequent water loss of a wet mix. Its application to mortars in-service would allow their subsequent mechanical properties, following stiffening, to be assessed. The bulk resistance of the mortar during and after dewatering will be influenced by factors such as the type and concentration of ions within the mix water. This is a function of the type of binder where cementitious binders, such as Portland cement, will contain a higher proportion of ions such as Cl -compared to calcium lime. The surface energy of the solid phases is also expected to influence the electrical characteristics of a mix. Consideration of these effects will become more significant particularly as more exotic additives derived from low embodied energy materials and industrial wastes become increasingly prevalent.

Conclusions

The results show a decrease in the exponential of (-bR b ) with square root of time for the mortars examined. This is consistent with a decrease in conductivity with loss of water and confirms that conduction is through the liquid as opposed to solid phase.

The time at which the change in slope occurs suggests the flow of water into the brick is controlled by the type of binder used to make the mortar.

The results demonstrate the ability of impedance spectroscopy to determine mortar bulk resistance and relate this to water flow into a brick.

7.

Tables Table 1. Values calculated from the experimental data relating to calcium lime, natural hydraulic lime (NHL) 3.5 lime and Portland cement mortars dewatered on Ibstock golden purple brick prisms. Square root of dewatering time, t 0.5 , min 0.5 5.7 6.4 7.7

Dewatering time (intercept of lines), min 32 42 60

List of Figure captions

  term, i, used in sharp front theory (equation 2) is comparable with h, rearranging this equation gives a 2 1 t relationship similar to that of equation 3, as

  mortar can be described as a system where solid particles of binder or aggregate are suspended in a liquid. The potential energy of ions dissolved in the liquid can be related to the number of ions by applying Boltzmann's distribution law of statistical thermodynamics, equation 8(Besanson 1974).

  and b are empirical parameters and R b is the resistivity. It should be noted that parameter b is temperature dependent.

  (a). The plot comprises two arcs which correspond to the mortar (left arc) and electrodes (right arc). The intersection of the arcs with the x-axis (real impedance) occurred at a frequency of 63 kHz and corresponds to the bulk resistance, R b . The intersection also occurred at this frequency for the natural hydraulic and cement mortars. An idealised plot and equivalent circuit diagram is shown in figure 2 (b).

  illustrates plots of the exponential of the recorded resistivity versus the square root of time. These may be analysed in terms of two regions relating to the mode of extraction of unbound or free flowing water. The initial slope describes the absorption of liquid from the wet mortar mixture into the brick by capillary forces and is governed by the chemical changes taking place.

  significant extent within the duration of the experiment. The absence of a structural network does not provide an increasingly restricted path thereby allowing a more rapid extraction of water. By contrast, for the natural hydraulic lime a very significant change in gradient was observed giving a very clear indication of the processes involved.Initially the hydraulic reaction from the incorporated dicalcium silicate phase produced a loose network of a gelatinous hydrated silicate limiting the process of free water extraction. At the intercept the remaining incorporated water is more tightly bound and the rate of the dewatering process is dramatically altered. It is possible that the loosely bound silicate gelatinous network may be drawn towards the mortar brick interface.Experiments with NHL3.5 mortars dewatered on a brick substrate have identified a crystalline silicate layer at the surface of the substrate. Clearly impedance measurements may give an indication of the time required to dewater the sample. Previous research (El-Turki et al 2009) has highlighted the importance of dewatering on the subsequent mechanical properties of a mortar.

Figure 1 .

 1 Figure 1. Experimental setup (a) Schematic diagram showing test cell, electrodes and brick prism. (b) Sharp front model of dewatering mortar.

Figure 2 .

 2 Figure 2. (a) Complex plane plot obtained after 30 minutes of dewatering for a calcium lime mortar showing the bulk resistance corresponds to a frequency of 63 kHz. (b) idealised complex plane plot and equivalent circuit model.

Figure 3 .

 3 Figure 3. Plots of exponential of (-b R b ) versus square root of time for Portland cement, natural hydraulic lime and calcium lime mortars showing two distinct linear regions.

Figure 4 .

 4 Figure 4. Typical exponential of (-b R b ) versus square root of time plot showing idealised structures for free and physically bound water.

Acknowledgements

The authors would like to thank the EPSRC for funding this work and Castle Cement and Ibstock Brick for supply of materials. Thanks are also due to Margaret Carter and Moria Wilson for helpful discussions.

8.


S Tardon

R Brüggemann
email: rudi.brueggemann@uni-oldenburg.de

Characterisation of the interface properties in a-Si:H/c-Si heterostructures by photoluminescence

Keywords: 78.55.-m, 72.20.Jv, 73.20.-r, 88.40.H-, photoluminescence, recombination, a-Si:H/c-Si, heterostructure, solar cell

published or not. The documents may come

Characterization of the interface properties in a-Si:H/c-Si heterostructures by photoluminescence

S Tardon, R

Introduction

Heterojunction solar cells based on hydrogenated amorphous silicon (a-Si:H) layers on crystalline silicon (c-Si) wafers have been developed over the last twenty years [1] and efficiencies in excess of 20% have been achieved [2]. A prerequisite to a high efficiency a-Si:H/c-Si heterojunction is a high quality amorphous/crystalline interface, as otherwise the solar cell properties are limited by recombination at the a-Si:H/c-Si interface [2][3][4][START_REF] Froitzheim | 3 rd World Conference in Photovoltaic Energy Conversion[END_REF][START_REF] Maydell | 31 st Photovoltaic Specialists Conference[END_REF][START_REF] Tardon | [END_REF]. Characterisation techniques, such as photoluminesence (PL) measurements, may be applied to study the influence of the interface defects on the minority carrier density in the volume of the cell [START_REF] Maydell | 31 st Photovoltaic Specialists Conference[END_REF][START_REF] Tardon | [END_REF].

The description of the photoluminescence (PL) radiation emitted from a semiconductor by Planck's generalised law [8,[START_REF] Würfel | Physik der Solarzellen[END_REF] gives access to the splitting of the quasi-Fermi levels (E Fn -E Fp ), where E Fn and E Fp are the electron and hole quasi-Fermi levels that describe the densities of free electrons and holes in the conduction and valence band, respectively. The energetic separation of the Fermi energies, (E Fn -E Fp ), in the volume of a wafer or solar cell is reduced via non-radiative transitions, such as bulk or surface/interface recombination via defects. The influence of surface recombination is noticeable in typical crystalline silicon (c-Si) solar cells, especially if the minority-carrier diffusion lengths are longer than the absorber thicknesses. It is thus possible to take the volume radiation recombination rate monitored by the PL measurement as an indicator of the interface quality that is of interest in a-Si:H/c-Si structures. Thus, we have applied such PL measurements for the characterisation of the interfaces for different wafer cleaning schemes and a-Si:H passivation layers [START_REF] Tardon | [END_REF]. The PL technique can also be applied to "precursor" structures after the different steps of solar cell processing, respectively, as has been shown on chalcopyrite cell structures [START_REF] Unold | Proceedings 16th European Photovoltaic Solar Energy Conference[END_REF]. These measurements have a strong relevance to the optimisation of solar cell performance since the magnitude of the quasi-Fermi level splitting gives an indication of the achievable open circuit voltage V OC in a final device [START_REF] Würfel | Physik der Solarzellen[END_REF].

Two aspects concerning the exploitation of PL measurements are dealt with in this paper: i) We analyse different steps in a standard production process of a-Si:H/c-Si solar cells with final AM1.5efficiencies > 14 % via quantitative photoluminescence measurements. ii) We apply PL measurements to obtain information on the interface defect distribution by studying p-type c-Si wafers in which the Fermi level positions at the amorphous / crystalline interface are adjusted by a variation in a-Si:H layer sequences. Numerical modelling assists in the interpretation of the results.

Theoretical Background

The spectrally dependent absorptivity A( ¢¡ ) of a material slab can be calculated as

£ ( ¤ ) £ ( ¤ ) 2 £ ( ¤ ) (1 ( ¥ ))(1 ( ¥ )e )(1 e ) ( ¥ ) 1 ( ¥ )R ( ¥ )e d d F B d F B R R A R - - - - + - = - ¦ ¦ ¦ § § § § § (1)
where R F ( ¨¡ ) and R B ( ¨¡ ) are the reflection coefficients at the front and back surface of the semiconductor, © ( ¡ ) is the absorption coefficient, and d is the thickness of the slab. For R F = R B , Eq.

(1) corresponds with Eq. ( 24) of [8].

Assuming that the carrier distributions of electrons and holes, n and p, can be described by quasi-Fermi level statistics and that the quasi-Fermi level splitting for electrons and holes, ∆E F = E Fn -E Fp , is independent of position, the photoluminescent photon flux dj , emitted from a semiconductor surface in an energy interval d( ¡ ) into one hemisphere is described by a generalised Kirchhoff law given by

[8,9] 2 2 2 ( ) ( ) d( ) d ( ) 4 ³c exp 1 Fn Fp A j E E kT = - - - ! ! ! ! ! ( 2
)
where k is the Boltzmann constant and c is the speed of light in vacuum. The absorptivity/emissivity term A( ¢"

) is that of Eq. (1).

Under typical excitation conditions, Eq. ( 2) reduces to

2 # 2 2 ( $ ) $ d ( $ ) exp exp d( $ ) 4% ³c Fn Fp E E A j kT kT ( ) ( ) 0 1 0 1 2 2 2 2 2 . ( 3
)
The exponential term with E Fn -E Fp reflects the product np/n i 2 for radiative recombination which shifts the spectral distribution without changing its shape, provided that E Fn -E Fp is spatially constant.

Experimental Details

The photoluminescence radiation was measured using a lock-in technique at room temperature using a laser diode with a wavelength λ = 782 nm as an excitation source, with a photon flux of 1.6 × 10 18 photons cm -2 s -1 (~ 4 × AM1.5). The photoluminescence was detected by a LN 2 -cooled

InGaAs detector positioned at the exit slit of a 0.25 m grating monochromator.

The calibration of the setup was performed using a tungsten filament lamp source placed at the sample position, allowing a quantitative measurement of the emitted photon flux in the photoluminescence experiment, and therefore the determination of the quasi-Fermi level splitting according to Eq. ( 2).

The samples were deposited at the Institut für Physikalische Elektronik in Stuttgart. Figure 1 sketches the sample series which was deposited on a double-sided polished p-type c-Si wafer (thickness 250 µm and 0.8-1 cm). It includes a symmetrical (i)a-Si:H/(p)c-Si/(i)a-Si:H structure (figure 1a) passivated with two 8 nm thin intrinsic a-Si:H layers and an (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(n)a-Si:H (figure 1b) structure. Technically, the structure of figure 1a received the additional nlayers so that the (i)a-Si:H/(p)c-Si interfaces are the same. These samples were designed in order to obtain symmetrical optical and structural properties at the back and front of the wafer with different Fermi levels at the interfaces. Figure 1c shows the precursor of the solar cell with different doped a-Si:H layers with an (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(p)a-Si:H sequence. Figure 1d shows the complete cell with the constituent layers deposited on the same remaining pieces of the wafer after respective surface cleaning and pre-treatment. Following these procedures, all pieces of the original wafer should have identical properties of the (p)c-Si/(i)a-Si:H interface.

Numerical modelling

Numerical modelling has been applied with the device modelling program SC-Simul [START_REF] Rösch | Proceedings of the 2 nd World Conference on Photovoltaic Solar Energy Conversion[END_REF][START_REF] Rösch | [END_REF]. It solves the full set of device and transport equations and has been successfully applied to a-Si:H/c-Si heterostructures [START_REF] Unold | [END_REF]14]. For the purpose of extracting information on the interface defect distribution we have simulated photoluminescence for each of the previously described types of samples and determined the electron and hole distributions from which E Fn -E Fp was calculated.

The defects at the interface N if were represented by a Gaussian profile of energetic width σ and peak energy E p, both of which are varied in the simulation study. The nature of the interface defect was taken as that of a dangling bond with its three possible charge states and corresponding densities of positive (D+), neutral (D0) and negative (D-) defects. The statistics in terms of a modified Shockley-Read-Hall approach to take account of the double-occupancy and correlation energy between the D+/0 and D0/-transitions was applied [15,[START_REF] Main | New Physical Problems in Electronic Materials[END_REF]. No approximations limiting the transitions to capture processes as in [START_REF] Olibet | [END_REF] have been introduced. The correlation energy considered here was U corr = 0.25 eV, similar to the a-Si:H structure [START_REF]Street Hydrogenated Amorphous Silicon[END_REF]. Band tails were not included in the representation of the interface recombination-centre distribution. Capture coefficient values for the electron and hole capture into charged defects were: c n D+ = 1 × 10 -7 cm 3 s -1 , c n D0 = 1 × 10 -8 cm 3 s -1 , c p D-= 1 × 10 -7 cm 3 s -1 and c p D0 = 1 × 10 -8 cm 3 s -1 , where the subscript denotes the carrier type and the superscript denotes the type of dangling bond into which the free carrier is captured. In the simulation, the a-Si:H/c-Si interface layer was represented by a 2 nm thin c-Si type defect layer with a chosen defect density, which is then translated to an interface defect density, N if . Good agreement between experimental and simulation results is achieved with N if = 6.4 × 10 11 cm -3 , which is the value in the employed simulations. A similar value has been determined in [START_REF] Froitzheim | [END_REF] from surface photovoltage experiments on a-Si:H passivated c-Si wafer with HF pre-treatment before the a-Si:H deposition.

Results

Photoluminescence spectra of a-Si:H/c-Si heterostructures

For the analysis of the PL spectra of samples with different front and back reflection properties it is necessary to properly take into account the absorptivity term in order to track the different outcoupling properties for the photoluminescence emitted from these structures. The absorptivity according to Eq.

(1) changes for the different layer sequences, especially after deposition of the metal back contact and the front TCO. Fig. 2 shows the calculated A( ¨" ) for three different cases: a c-Si wafer, passivated with a-Si:H layers on both sides and with identical front and back reflection coefficients, the wafer with an additional Al back contact and the final device with a ZnO front contact. The reflection coefficients were assumed to be wavelength independent. The connection to the experimental value is based on the measured reflection in the high-energy range of Fig. 2. For the solely a-Si:H covered surfaces, we have determined R = 0.308 neglecting the effect of the thin amorphous silicon layers. In the shown spectral range of Fig. 2, R F = 0.2 for the ZnO covered surface.

Figure 2 shows that at high photon energy A( ¨" ) according to Eq. ( 1) saturates at 1 -R F . Here, the product of © d is much larger than 1 so that the exponential terms can be neglected and only 1-R f remains in the expression for A( ¨" ) in Eq. ( 1) [8]. At lower photon energies, re-absorption is less important so that the back reflection properties determine the output. As we consider the sample with the Al back contact having mirror properties and R B = 1, A( ¢" ) increases for energies ¢" < 1.25 eV, compared with the symmetric structure. If the front reflection is less, as in the case of the ZnO layer, the output is even higher. It is noted that for the same quasi-Fermi level splitting, the measured PL intensity will be different solely due to the different absorptivity/emissivity spectra..

The front reflection of the structure also influences the photon flux that enters through the front interface. This variation must be taken into account when the PL spectra of Fig. 3 for the different structures are compared and discussed. For the symmetric structures there is a distinct increase in the quasi-Fermi level splitting of about 25 meV which will be discussed in more detail in the next section. For the processing steps of the solar cell denoted by the full symbols it can be seen that E Fn -E Fp increases after each step. For the ZnO layer sample, the excitation photon flux is increased as the front reflection of the ZnO is 0.1 at the photon energy of the excitation laser, in contrast to the 0.3 for the a-Si:H front layer sequences. We can thus calculate the increase in the quasi-Fermi level splitting, which results from the increased excitation flux.

Translating this increment of the photon flux into the splitting of the quasi-Fermi levels, assuming that each absorbed photon generates one electron -hole pair, yields ,ZnO 0 0 0.9 ( ) ln ln 7meV, 0.7

Fn Fp E E kT kT γ γ ¡ Φ ¡ Φ ∆ - = = = ¢ £ ¢ £ ¢ £ Φ Φ ¤ ¥ ¤ ¥
where Φ ,ZnO (Φ γ ) is the photon flux entering the silicon wafer with (without) the ZnO layer and Φ 0 is the photon flux impinging on the wafer. Fig. 3 shows that, in good agreement with these considerations, the experimental increase is 8 meV between the structure with and without the ZnO, i.e., the increase of the quasi-Fermi level splitting and of the PL output is solely due to the higher excitation level.

Figure 3 shows that the deposition of the Al layer leads to a further increase in PL yield over the whole spectrum. From the consideration for A from Fig. 2, we would expect an increase in the PL yield but we would not expect an increase in the high-energy range for photon energies > 1.25 eV if the quasi-Fermi level splitting remains unchanged. However, the experimental result for our structure is an increase in E Fn -E Fp of 6 meV from 661 meV to 667 meV after the Al deposition, which indicates an increased minority carrier density. We note, that despite the error of 4 meV for each value of E Fn -E Fp , the hierarchy in PL yield and the relative increase, respectively, from the different structures has always been maintained.

Simulation results for the symmetrical structures

The band diagrams of the two symmetric structures are depicted in Figs. 4(a) and 4(b). For the structure with the additional doped layer, i.e., the layer sequence (n)a-Si:H/(i)a-Si, the band bending is dominated by the doped n-type amorphous silicon because the intrinsic a-Si:H layer is very thin with a low charge density and its effect can be neglected. By comparison of both figures it can be seen that the occupancy and charge of the defects must change because the relative position of the interface Fermi level changes.

Figure 5 shows simulation results of the interface recombination rate R if (Eq. (A1)) for the (i)a-Si:H/(p)c-Si/(i)a-Si:H and (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(n)a-Si:H structures when the peak position E P of the defect distribution of width = 0.2 eV varies with respect to the valence band edge,

(E P -E V ).
The values of R if are higher for the structure with the additional n-layers, apart from the values at which E P is close to the c-Si conduction band edge.

There is a strong anti-correlation between the interface recombination rate and the splitting of quasi-Fermi levels, taken here as an average value from the c-Si volume, as depicted in Fig. 6. As expected, high interface recombination results in a decrease in minority carrier density and thus reduces the quasi-Fermi level splitting in the volume of the wafer (and vice versa). For the passivation with (i)a-Si:H, E Fn -E Fp reaches a minimum at mid-gap while adding (n)a-Si:H increases E Fn -E Fp when defects are moved towards the conduction band.

From the experimental results of Fig. 3, a difference of (25 ± 8) meV between the splitting of quasi-

Fermi levels E Fn -E Fp is observed. A value of (600 ± 4) meV is measured when the wafer is passivated with (i)a-Si:H, which increases to (625 ± 4) meV with the addition of an (n)a-Si:H layer on both sides.

Inspection of Fig. 6 shows that these experimental findings are met when E P -E V lies in the range 0.7 to 0.8 eV.

To take account of the nature of the defects and their energetic distributions, Fig. 7 = 8 × 10 -9 cm 3 s -1 , adopted by Froitzheim in the simulation part of his work [START_REF] Froitzheim | [END_REF]. The full symbols show a similar variation in E Fn -E Fp with a slightly larger E Fn -E Fp difference between the two structures. The circle in Fig. 7 indicates the range of E P -E V for which a variation in the parameters leads to the experimentally determined difference in E Fn -E Fp of about 25 meV.

Complementary results for asymmetrical structures

Discussion

Our analysis takes into account the different aspects of the variations in the coupling of light into and out of the wafer structures, and we have adjusted the analytical spectra with different A( ¨" ) to obtain a correct description of the absorptivity/emissivity, as shown in Fig. 3. This gives an insight into why changes in the quasi-Fermi levels occur after different processing steps, notably after the ZnO and Al deposition. When the ZnO layer is added, there is an increase in the excitation and emission flux due to variations in the reflectivity and absorptivity of the wafer which causes a higher quasi-Fermi level splitting. However, the analysis has shown that this increase is solely attributed to the higher excitation flux after ZnO deposition and not to any intrinsic changes in internal field.

Similarly, reflection of PL radiation at the Al back contact increases the PL flux, which is coupled out. Because of re-absorption, there is no additional flux in the high energy range of the PL spectrum in which E Fn -E Fp is usually determined. This has been demonstrated by placing a separate mirror at the back side of a symmetric structure [14]. Thus, after Al deposition, no increase in E Fn -E Fp is expected.

However, a higher E Fn -E Fp has been detected for the structure with the Al back contact. This can be accounted for by a small increase in the doping level at the back interface by diffusion of Al into the (p)a-Si:H layer, and consequent enhancement of the passivation quality at the back of the wafer.

Further systematic investigations (e.g. by SIMS profiling) will be needed to fully understand this effect.

Concerning the symmetric structures in Fig. 3, the quantitative PL yields an increased output when comparing the (i)a-Si:H/(p)c-Si/(i)a-Si:H and (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(n)a-Si:H structures.

The appeal of this comparison is that the (i)a-Si:H/(p)c-Si interface defect distribution is the same in both structures. The only difference would be the changed occupancies under illumination which are influenced by the highly-doped (n)-a-Si:H layer, and are included in the numerical solutions for the whole system.

The simulation results have shown that there is a variation in the quasi-Fermi level splittings of the two structures when the simulation parameters are varied slightly. This is observed e.g. for different widths of the distribution in Fig. 7. The simulations show also that a flat distribution, with Gaussian width > 0.5 eV, does not predict the experimental result but in fact to the opposite effect, with a higher PL yield for the (i)a-Si:H structure. For very narrow distributions with width < 0.01 eV the difference in quasi-Fermi level splitting increases to much higher values than observed experimentally.

It is not the aim of the simulation study to track in detail all the features of the experimental results.

Given the scatter in the experimental values, which depend on the different measurement positions on the sample, and the range of plausible simulation parameter values one can say that the experimental trend is very well reproduced, especially since the simulation results agree with the experimental trend across the three different structures in Fig. 8. Also, the good agreement in the open-circuit voltage between the experimental and simulated solar cell (d of Fig. 1) adds confidence in the conclusion that, within the model we have used, the interface defect peak lies in the range 0.7 eV -0.8 eV above the c-Si valence band edge, with a Gaussian width around 0.1 eV.

Conclusions

An advantage of photoluminescence measurements, as a contactless method, is to enable the analysis of consecutively processed layers during device deposition, and thereby the evaluation of each of the process steps. By taking into account the different contributions to the PL yield, from the excitation flux to the two terms in Planck's generalised law -quasi-Fermi level splitting and absorptivity -we have distinguished between cases in which a higher PL flux is related to a larger quasi-Fermi level splitting, i.e., minority carrier density, and those where other factors are involved. While the higher PL yield after ZnO deposition is solely related to a higher excitation flux, the deposition of Al was shown to result in a higher quasi-Fermi level splitting with improved back-interface properties.

From the differences in relative and absolute terms in the PL yield of the symmetric (i)a-Si:H/(p)c-Si/(i)a-Si:H and (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(p)a-Si:H structures with the same (i)a-Si:H/(p)c-Si interface properties, we deduce the energetic position of the interface defect distribution. Although the energetic value is not unambiguous, simulations using a range of plausible parameter sets indicate a position of 0.7 to 0.8 eV above the c-Si valence band edge and a Gaussian width in the range of 0.1 eV. The analysis is supported by electrical measurements of the solar cell open-circuit voltage, which could also be reproduced with the same simulation parameter set.

Fig. 3

 3 shows the spectral quantitative PL-yield measured from the different structures. The black dotted lines are analytically derived from Eq. (2). Adjustments of the analytical spectra yield an error of 4 meV for the deduced values of E Fn -E Fp . Full symbols show the results from asymmetric structures with differently doped layers and Al and ZnO additions. The open symbols represent measurements with the symmetrical structures.

Fig. 6 .

 6 Fig. 6. The full symbols are from additional simulations with σ = 0.1 eV and slightly different sets of capture coefficients c n D+ = 8 × 10 -8 cm 3 s -1 , c n D0 = 4 × 10 -8 cm 3 s -1 , c p D-= 4 × 10 -8 cm 3 s -1 and c p D0

Figure 8

 8 Figure 8 summarises experimental and simulation results for the quasi-Fermi leve splitting in the symmetric structures, where the range of experimental data points results from measurements at different positions on the sample. The asymmetric structure with a "back surface field" layer (structure c of Fig. 1) shows a further increase in E Fn -E Fp . The simulation results, determined with the same parameter set, reproduce the correct trends. In order to substantiate the agreement between the simulated and the experimental results, we have also measured electrically the open-circuit voltage of the solar cell (structure d of Fig. 1) under AM1.5 illumination, and obtain a value of 590 mV. Simulation of the solar cell with identical parameters to those used in the wafer modelling results in V OC = 591 mV, under the same AM1.5 illumination conditions.

Figure 1 :

 1 Figure 1: Illustration of the samples investigated: a) c-Si absorber passivated symmetrically with two 8 nm thin (i)a-Si:H layers, b) the structure of a) with additional two 8 nm (n)a-Si:H layers, c) asymmetric (n)a-Si:H/(i)a-Si:H/(p)c-Si/(i)a-Si:H/(p)a-Si:H and d) (n)a-Si:H/(p)c-Si solar cell with symmetric 8 nm (i)a-Si:H buffer layers and a 20 nm (p)a-Si:H "back surface field" layer and an Al back contact.

Figure 2 :

 2 Figure 2: Calculated absorptivity of the (p)c-Si wafer for different reflection coefficients at the front and back. For the samples according to Fig.1 (a-c), the same reflectivity at the front and back of the wafer has been assumed (R F = R B = 0.308) (dotted line). The addition of a ZnO layer reduces the reflection at the front in the spectral range shown (R F = 0.2), which results in an increase of A at higher energies from 0.692 to 0.8 (dashed line). The addition of a back contact as in Fig.1(d) produces a reflector (R B = 1) and A increases for energies ¢" < 1.25 eV (full line).

Figure 3 :

 3 Figure 3: Quantitative photoluminescence after different steps of c-Si wafer and a-Si:H/c-Si solar cell processing (full symbols) and of symmetric structures (open symbols) with the corresponding quasi-Fermi level splitting. The dotted lines are calculated with the analytical expression of Kirchhoff´s generalised law (Eq. 2).

Figure 4 :Fig. 1 .

 41 Figure 4: Calculated band diagram in thermodynamic equilibrium of the 2 symmetric structures of Fig. 1. In the diagram, the energetic distribution of defects at the interface is shown schematically.

Figure 5 :

 5 Figure 5: Total interface recombination rate versus energetic position of the defect peak relative to the valence band edge.

Figure 6 :Figure 7 :Figure 8 :

 678 Figure 6: Corresponding splitting of quasi-Fermi levels for different energetic positions of the defect peak at the interface for the (i)a-Si:H (squares) and (n)a-Si:H/(i)a-Si:H (circles) passivated wafers. The simulation has been performed for two different parameter inputs: = 0.2 eV and (open squares and circles) and = 0.1 eV (crossed squares and circles). The capture coefficients were also different (see text). The experimental results can be reproduced in the numerical calculation by an energy peak around 0.7 eV for interface defects

Figure A1 :

 A1 Figure A1: Transition diagram for the recombination via dangling bonds showing the possible transitions of electrons and holes for the three charge conditions. Emission is denoted by dashed arrows, free-carrier capture by full arrows.

  Figure 1

  Figure 8

& ' & ' = -

Acknowledgements

The authors thank G.H. Bauer for valuable discussions and suggestions. They also thank their colleagues at IPE, Stuttgart, for sample preparation and M. Scherff at FernUniversität Hagen for complementary measurements. Financial support by BMBF contract 01SF0115 and EWE Stiftung is gratefully acknowledged.

Appendix

The statistics of the interface defects is taken to follow the well-established statistics developed for dangling bonds in a-Si:H which takes into account the three charge states positive, neutral and negative [15,[START_REF] Main | New Physical Problems in Electronic Materials[END_REF]. Thus, the defect density at an energy (E) is the sum of the densities in the different charge states D + (E), D 0 (E) and D -(E). The electron interface-defect recombination rate R if is given by the net loss via two paths as illustrated in Fig. A1. Taking into account the capture and emission rates, R if densities can be derived [15,[START_REF] Main | New Physical Problems in Electronic Materials[END_REF].


N Bonfils

N Chevaugeon

N Moës

Treating volumetric inequality constraint in a continuum media with a coupled X-FEM/Level-Set strategy

Some mechanical problems involve inequality kinematic constraint. This study deals with an original approach to handle those difficult problems. The main issue is the treatment of the variational inequalities due to the fact that the constrained area is a priori unknown. The method, introduced here, is to find the exact constrained area iteratively starting from an intial trial one. Thanks to numerical tools such as level-set and X-FEM we turn the constrained minimization problem into a shape equilibrium problem.

Introduction

The goal of this study is to propose a new approach to deal with continuum media subjected to a kinematic inequality constraint. A kinematic constraint is a restriction on the deformations a body can undergo. In this document, we will treat the problem of a media subjected to a volumetric inequality constraint. In fact, this paper is first step toward handling an extensibility limit constraint in hyperelastic materials.

Classical methods used to enforce inequality constraint within a media are the Lagrange multiplier method [START_REF]Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF], the penalty method [START_REF]Computational Contact Mechanics[END_REF], and the augmented Lagrangian method [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. The penalty method does not permit to enforce exactly the kinematic constraint and can lead to a bad conditioning of the tangent stiffness matrix. Both the Lagrangian and augmented Lagrangian method assure the respect of the constraint. However, their numerical treatment imply the use of specific algorithms such as the Uzawa method [START_REF]Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] or the Generalized Newton Method [START_REF] Curnier | A generalized Newton method for contact problems with friction[END_REF], and parameters have to be specified. These parameters play a major role in the convergence of those algorithms and it is often very difficult to choose them appropriately. Note that in most cases, those algorithms are used to treat inequality constraint on the boundary of a media (unilateral contact condition). The constraint studied here is a volumetric constraint.

After a brief statement of the constraint minimization problem studied here, we will describe in the first section an already existing method [START_REF] Ben Dhia | Problèmes de contact frottant en grandes transformations : du continu au discret[END_REF] (close to augmented lagrangian method) used for contact problems and that can be applied to our problem. This is an interesting method since it does not need the use of a Generalized Newton Method and can be solved using a kind of Active Set algorithm [START_REF]La méthode des contraintes actives appliquée au contact unilatéral[END_REF]. However, we will show that in our particular case, this method is not really efficient. First, it offers poor rate of convergence and the presence of user-set parameters for the active set algorithm may deteriorate its stability.

In the next section, we will present our new approach. The strategy consists in turning the minimization problem under inequality constraint into a shape equilibrium one. The inequality constraint induces variational inequalities because we do not know a priori the shape of the active constraint and unactive constraint zone in the media. The strategy proposed here is to replace the inequality constraint by an equality constraint over an unknown domain. This domain is called the active constraint domain. Hence, we assimilate the media to a two elastic phases media. On one phase, no constraints are imposed and on the second one an equality constraint is imposed. Each phase has a mechanical behaviour that depends on the state of the constraint (active or inactive). With the help of the configurational mechanic and the notion of configurational force, the exact shape of each phase is found iteratively. Finally, the notion of level-set [START_REF]Level Set Methods and Fast Marching Methods : Evolving Interfaces in Computational Geometry[END_REF] and the X-FEM [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] will be used to give an implicit description of the interface separating the two phases. Several numerical tests will be shown in the final section so as to compare the first method to the new one.

Statement of the problem

A linear elasticity two-dimensional problem is considered, with an isotropic homogeneous material in a plane strain state. Consider the domain Ω as shown in Figure 1. On the boundary, a Dirichlet condition u = u d is applied on Γ u and a Neumann condition σ.n = t d is prescribed on Γ t . Finally, the whole domain is submited to a body force f .

In this study we will adopt a particular form of the constraint in the small strain setting. The constraint reads: γ(ε) = Trε -α ≤ 0 on Ω with ε, the classical strain tensor and Trε the trace of ε.

V = {v ∈ H 1 (Ω)|v = u d on Γ U } (1)
V 0 = {v ∈ H 1 (Ω)|v = 0 on Γ U } (2)
And the total potential energy of the mechanical system:

Π(u) = Ω ψ(ε)dΩ - Ω f .udΩ - Γ t t d .udΓ (3)
where ψ(ε) = 1 2 ε : C : ε and C is the fourth-order elasticity tensor. We are looking for the infimum of Π(u)

inf v∈S Π(v), S = {v ∈ V | γ(ε(v)) ≤ 0} (4)
The space S is the space of kinematically admissible displacements. The problem ( 4) is called the primal problem. The introduction of the inequality constraint in the search of the infimum yields to a variational inequality [START_REF]Les inéquations en mécaniques et en physique[END_REF].

If u is an infimum, then

Ω ε(u) : C : ε(v -u)dΩ ≥ Ω f .(v -u)dΩ + Γ T t d .(v -u)dΓ ∀v ∈ S (5)
This result is a variational inequality characterizing the solution u of the problem (4). This is also the primal variational formulation of this problem. Unlike classical variational principles encountered in solid mechanics, we obtain an inequality because of the constraint u ∈ S.

A classical approach

The variational formulation

Constrained optimization problems like (4) are difficult to formulate and to solve. Indeed the inequality constraint prevents us from solving it directly. The goal of this section is to transform a constrained problem into an unconstrained one. In this section, we will recall how our constrained minimization problem can be transformed into saddle point problem using the method of Lagrange multiplier [START_REF]Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] (duality principles). First, to take into account the constraint γ(ε), we add to the elastic strain energy function ψ(ε) the following characterisic function

ψ(ε) = ψ r (ε) + I - γ(ε) (6)
where ψ r , is the parent or regular potential (ψ r (ε) = 1 2 ε : C : ε) and I -(x), is the characteristic function of R -such that

I - x = 0 if x ≤ 0 +∞ otherwise (7)
The constitutive law then reads:

σ = C : ε + pI (8) γ(ε) ≤ 0 , p ≥ 0 , γ(ε)p = 0 (9)
Now, let the Lagrange multiplier trial functions be q. These functions reside in the following space

Q + = {q ∈ Q|q ≥ 0}, Q = {q ∈ L 2 (Ω)} (10)
Since

I - γ(ε) = sup q∈Q + {qγ(ε)} (11)
the above infimum problem becomes a saddle-point problem with the following functional :

L(v, q) = Ω ψ r (ε) + qγ(ε) dΩ - Ω f .vdΩ - Γ t t d .vdΓ (12)
We are looking for the saddle-point of this functional

inf v∈V sup q∈Q + L(v, q) (13)
Unfortunately, there is still an inequality constraint. Indeed, the lagrange multiplier q has to respect the positivity condition q ≥ 0. We still have to solve a constrained problem.

If (u, p) is the saddle point of the functional L(v, q), then

DL(u, p)[v] = 0 ∀v ∈ V 0 (14) DL(u, p)[q -p] ≤ 0 ∀q ∈ Q + ( 15
)
where DL(u, p)[qp] denotes the directional derivative of L(u, p) with respect to p in the arbitrary direction (qp) and DL(u, p) [v] the directional derivative of L(u, p) with respect to the motion in the arbitrary direction v.

Rewriting ( 14) and ( 15), we obtain once again a variational inequality of the form [START_REF]Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]:

Ω ε(u) : C : ε(v)dΩ + Ω pTrε(v)dΩ = Ω f .vdΩ + Γ T t d .vdΓ ∀v ∈ V 0 (16) Ω γ(ε(u))(q -p)dΩ ≤ 0 ∀ q ∈ Q + ( 17
)
An alternative is to use an Augmented Lagrangian formulation [START_REF]Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF] which is a composition of an ordinary lagrangian method and a penalty method. This alternative method permits to regularize the saddle point problem [START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented Lagrangian treatment[END_REF] and to give an unconstrained problem. Then, we augment once again the functional with a term that ensures the positivity condition [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF].

The functional becomes :

L + (v, q) = Ω ψ r (ε)dΩ + 1 2ρ Ω q + ργ(ε) 2 + dΩ - 1 2ρ Ω q 2 dΩ - - Ω f .vdΩ - Γ t t d .vdΓ (18)
where ... + denotes the positive part of the considered quantity, and ρ is a postive parameter, called augmentation parameter.

The saddle point problem becomes

inf v∈V sup q∈Q L + (v, q) ( 19
)
The constraint on q has been removed and the stationnary condition of the above functional reads:

If (u, p) is the saddle point of the functional L + (v, q), then

DL + (u, p)[v] = 0 ∀v ∈ V 0 (20) DL + (u, p)[q] = 0 ∀q ∈ Q (21)
The stationary condition leads to the following formulation

Ω ε(u) : C : ε(v) + ρχ(g)Trε(u)Trε(v) dΩ + Ω χ(g)pTrε(v)dΩ = Ω f .v + χ(g)ραTrε(v) dΩ + + Γ T t d .vdΓ ∀v ∈ V 0 (22) -1 ρ Ω (1 -χ(g))pqdΩ + Ω χ(g)Trε(u)qdΩ = Ω χ(g)αqdΩ ∀q ∈ Q (23)
where

χ(x) = 0 if x < 0 1 if x ≥ 0 ( 24
)
and g designates the augmented multiplier [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]

g = p + ργ(ε) (25)
The parameter ρ does not affect the correctness of the theorical formulation, but from a numerical point of view its choice does affect the convergence of the scheme and the condition number of the tangent stiffness matrix.

The strong form associated to this problem is depicted in Table 1.

Ω mechanical equilibrium div(σ) + f = 0 σ.n = t d constitutive relation σ = C : ε + χ(g)gI p -χ(g)g = 0 g = p + ρ(Trε -α) kinematic conditions u = u d ε = 1 2 (∇(u) + ∇(u) t )
Table 1: The strong form

Numerical strategy 3.2.1. The active-set algorithm

Starting from the formulation stated previously, we have to solve it numerically. Firstly, this formulation represents a set of non-linear equations even if our setting is the linear elasticity. Indeed, the function χ presented above introduces a non-linearity. This non-linearity due to the constraint condition could be solved with a Generalized Newton Method [START_REF] Curnier | A generalized Newton method for contact problems with friction[END_REF], but we will use here another method. A kind of fixed point algorithm described in [START_REF]Une formulation pour traiter le contact frottement en 3d dans Code Aster[END_REF][START_REF] Ben Dhia | Problèmes de contact frottant en grandes transformations : du continu au discret[END_REF][START_REF] Ben Dhia | Hybrid frictional contact particles-in elements[END_REF][START_REF]Eléments mixtes de contact frottant en grandes transformations et applications[END_REF][START_REF]Contribution à la modélisation mécanique et numérique des problèmes de contact-impact[END_REF][START_REF]Approche X-FEM pour la fissuration sous contact des structures industrielles[END_REF] is here used. This algorithm is close to the active set method [START_REF]La méthode des contraintes actives appliquée au contact unilatéral[END_REF][START_REF]Algorithme des contraintes actives et contact unilatéral sans frottement, G. Dumont[END_REF] and permits to update iteratively the values of χ. This active-set algorithm is explained below

Algorithm 1 Active-set algorithm Initialization of χ (0) = 0 loop Computing of u i et p i Updating χ i = χ i (u i , p i ) if χ i = χ i-1 (convergence criterion) then break end if end loop
At each step of the loop, the field χ is completly known over the whole domain. We can then solve a linear problem, and then update the state of χ. In this algorithm, all the states of the χ function are changed at the same time. (this function is evaluated at the Gauss points since we use a Gauss quadrature integration). Unfortunaletly, the convergence of this algorithm is proved only if just one state is changed at each iteration. We prefer here to change all the states at the same time [START_REF]Eléments mixtes de contact frottant en grandes transformations et applications[END_REF] so as to obtain a significant gain in speed.

Convergence criterion We evaluate the χ function at the Gauss points during the numerical integration procedure. The active-set algorithm has converged when the values of χ at the Gauss points has not changed. Let χ i the vector containing the values of χ of all Gauss points at iteration i of the algorithm. The algorithm has converged when χ iχ i-1 = [0]. Unfortunately, in this algorithm, cases of looping between several states of χ may occur. These oscillations prevent the algorithm from converging and may deteriorate the accuracy of the results.

The discretization

Once the non-linearity due to the presence of the inequality constraint is solved, we have to establish the spatial discretization. We denote p h et u h the discrete Lagrange multiplier field and discrete displacement field. The discretized weak form at iteration i of the active set algortithm yields Find (u h , p h ) such that

a(u h , v h ) + b(p h , v h ) = l(v h ) ∀ v h ∈ V h (26) b(u h , q h ) + c(p h , q h ) = p(q h ) ∀q h ∈ Q h ( 27
)
with

a(u i h , v h ) = Ω h ε(u i h ) : C : ε(v h )dΩ + Ω h χ i-1 ρTrε(u i h )Trε(v h )dΩ (28) b(p i h , v h ) = Ω h χ i-1 p i h Trε(v h )dΩ (29) l(v h ) = Ω h f .v h dΩ + Γ T t d .v h dΓ + Ω h χ i-1 ραTrε(v h )dΩ (30) c(p i h , q h ) = -1 ρ Ω h (1 -χ i-1 )p i h q h dΩ (31) b(u i h , q h ) = Ω h χ i-1 Trε(u i h )q h dΩ (32) p(q h ) = Ω h χ i-1 αq h dΩ ( 33
)
where Q h and V h are the discrete finite dimensions subspaces of Q and V .

The discretization of the weak form involves two fields : the displacement u and the Lagrange multiplier p. The shape functions for the Lagrange multiplier differ from those for the displacements. Hence, we will use different symbols for the two approximations. The classical FEM discretization on an element Ω e reads

u h | Ω e = N ∑ i=1 u h i N u i p h | Ω e = N ′ ∑ i=1 p h i N p i ( 34
)
where N and N ′ is the number of coefficients describing the approximation, and N i is the shape function associated to the coefficient i. If the use of a Lagrange multiplier field ensures the respect of the constraint, its discretization involves a major problem. Indeed, the interpolation used has to respect the well-known Ladyzhenskaya-Babuska-Brezzi compatibility condition [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] which guarantees the stability of the mixed method. Denoting

b ′ (v h , q h ) = Ω Trε(v h )q h dΩ (35)
This condition reads inf

q h ∈Q h sup v h ∈V h b ′ (v h , q h ) v h 1 q h 0 = β h ≥ β > 0 (36)
Where . 1 and . 0 designate respectivlely the H 1 and L 2 (Ω) norms, and β is a positive real independent of the mesh size. Among all mixed elements, only a few respect this condition. Moreover, we have used mixed elements with discontinuous pressure because we have noticed that this kind of elements permits the stability of the active-set algorithm. These elements can be used because we do not need the derivative of the pressure field in the weak formulation. One basic discontinuous pressure element has been selected. The mixed element chosen is the P2 + /P1 disc (see Figure 2). This element is composed of a quadratic displacement enriched by a cubic bubble function and a linear discontinuous pressure. The upper script + designates the enrichment by the bubble function. The displacement is enriched with this special function so as to increase the degree of the displacement interpolation and to ensure the stability of the formulation.

More preciselly, the shape functions for each component of the displacement field are the regular second order lagrange shape function (P2), plus a bubble function. The bubble function that we used is defined in each triangle as : b

(x) = λ 1 (x)λ 2 (x)λ 3 (x)
where the λ i are the barycentric coordinates relative to the three nodes of the triangle. The pressure field is defined independently in each element by a linear combination of the three order one lagrange shape function on the triangle. This finite element formulation therefore give a total number of degree of freedom ndo f = 2 * nv + 2 * ne + (2 + 3) * nt, where nv is the number of vertices, ne is the number of edges and nt is the number of triangle.

A new approach

The strategy

The essential difficulty of the problem described in section 2 is that the region where the constraint is active is not known a priori. This observation has motivated a new strategy of resolution. We consider here the same domain Ω and the same boundary conditions as before. But the strategy proposed here is to split the whole domain Ω into two distinct subdomains separated by Γ. The first domain denoted by Ω c is the constrained zone or active zone where an equality constraint is imposed. The second domain denoted by Ω f is the unconstrained or unactive zone, where no constraints are imposed. On each subdomain, the same linear elastic isotropic homogeneous material is considered. But in the domain Ω c we impose the following constraint

Ω c ∪ Ω f = Ω (37) Ω c ∩ Ω f = / 0 (38) ∂Ω c ∩ ∂Ω f = Γ (39)
γ(ε) = Trε -α = 0 on Ω c ( 40
)
The elastic energy in domains Ω f and Ω c are given by

ψ c = 1 2 ε : C : ε + I 0 γ(ε) = µε D : ε D + κ 2 α 2 + I 0 γ(ε) (41)
ψ f = 1 2 ε : C : ε = µε D : ε D + κ 2 (Trε) 2 (42)
where κ = λ + 2 3 µ is the bulk modulus, with λ and µ the Lamé coefficients and I 0 (x), is the characteristic function of 0 such that

I 0 (x) = 0 if x = 0 +∞ otherwise (43)
The corresponding constitutive laws then read

σ c = 2µε D + καI + pI (44) σ l = 2µε D + κTrεI (45)
Note that p is not the classical pressure (Trσ), it is what is needed for σ to be in equilibrium We can notice that the behaviour in domain Ω f and Ω c no longer involve inequalities and differ from the original behaviour [START_REF] Curnier | A generalized Newton method for contact problems with friction[END_REF].

The total potential energy Π ′ (u) of the system reads

Π ′ (u) = Ω c ψ c (ε)dΩ + Ω f ψ f (ε)dΩ - Ω f .udΩ - Γ T t d .udΓ (46)
The interface between the two domains can be seen as a surface of discontinuity and the strategy is to propagate correctly the front separating the two elastic phases Ω c and Ω f so as to find the exact location of Γ. Then, we have replaced the initial problem of a media subjected to a kinematic inequality constraint by a problem of a media subjected to an equality constraint over an unknown restricted domain Ω c .

The interface Γ is a coherent interface, thus the jump of the displacement field u has to vanish across the interface Γ, preventing separation, penetration or slipping of the two domains. We do have the continuity of the normal stress vector.

Let Γ be fixed for the moment. We are looking for the infimum of Π ′ (u)

inf v∈V Π ′ (v) (47)
Rewriting (47) as the following saddle point problem

inf v∈V sup q∈Q L ′ (v, q) (48)
where L ′ (v, q) designates the following Lagrangian functional

L ′ (v, q) = Ω c ψ c (ε) + qγ(ε) dΩ + Ω f ψ f (ε)dΩ - Ω f .udΩ - Γ T t d .udΓ (49)
The stationary condition of the above functional leads to the variational form Find (u, p) such that:

a(u, v) + b(p, v) = l(v) ∀v ∈ V 0 b(u, q) = p(q) ∀q ∈ Q (50) a(u, v) = Ω ε(u) : C : ε(v)dΩ (51) b(p, v) = Ω c pTrε(v)dΩ (52) l(u, v) = Ω f .vdΩ + Γ T t d .vdΓ (53) b(u, q) = Ω c
Trε(u)qdΩ (54)

p(q) = Ω c αqdΩ (55)
The strong formulation associated to this new problem reads: Unlike the previous approach, we obtain a fully linear set of equations. Indeed, the nonlinearity introduced by the function χ has been removed. The intrinsec non-linearity due to the presence of the inequality constraint will be solved by changing iteratively the shape of each phase. Note that this formulation does no involve any parameter since the equality constraint on Ω c is solved using an ordinary Lagrangian method.

Ω c Γ Ω f mechanical equilibrium div(σ) + f = 0 [[σ.n]] = 0 div(σ) + f = 0 σ.n = t d on ∂Ω c ∩ Γ T σ.n = t d on ∂Ω f ∩ Γ T constitutive relation σ = C : ε + pI σ = C : ε Trε -α = 0 kinematic conditions u = u d on ∂Ω c ∩ Γ U [[u]] = 0 u = u d on ∂Ω f ∩ Γ U ε = 1 2 (∇(u) + ∇(u) t ) ε = 1 2 (∇(u) + ∇(u) t ) ε = 1 2 (∇(u) + ∇(u) t )
We denote

ψ(ε(x), x) = ψ c (ε c (x)) if x ∈ Ω c ψ f (ε f (x)) if x ∈ Ω f ( 56
)
Let us consider now the directional derivative of Π ′ (u) for a configurational change θ of Γ. Here, u is the solution of the infimum problem (47) for a given location Γ of the interface. We could then write Π ′ (u) = Π ′ (u(Γ), Γ).

Following [START_REF] Pradeilles-Duval | Mechanical Transformations and discontinuities along a moving surface[END_REF], we get:

DΠ ′ (u(Γ), Γ)[θ] = Ω f Dψ f [θ]dΩ + Ω c Dψ c [θ]dΩ + Ω f ψ f D(dΩ)[θ] + Ω c ψ c D(dΩ)[θ] - - Γ t t d .Du[θ]dΓ - Ω D f [θ].udΩ - Ω f .Du[θ]dΩ - Ω f .uD(dΩ)[θ] (57) DΠ ′ (u(Γ), Γ)[θ] = Ω f σ : ε(Du[θ])dΩ + Ω c σ : ε(Du[θ])dΩ - Γ t t d .Du[θ]dΓ - Ω f .Du[θ]dΩ + + Ω f (ψ f -f .u)∇.θ -σ : (∇u.∇θ) -D f [θ].udΩ + + Ω c (ψ c -f .u)∇.θ -σ : (∇u.∇θ) -D f [θ].udΩ (58)
The first four term above in the right hand side cancel out because of the equilibrium. Then noting the expression of the Eshelby tensor [START_REF] Eshelby | The elastic energy-momentum tensor[END_REF]

P, P = (ψ(ε) -f .u)I -∇ u T .σ (59)
we get

DΠ ′ (u(Γ), Γ)[θ] = Γ θ.[[P]].ndΓ (60)
After calculus (Appendix A), we get

DΠ ′ (u(Γ), Γ)[θ] = Γ τθdΓ (61)
where

τ = (λ + 2µ) -1 p 2 c 2 = (λ + 2µ) 2 (Trε f -α) 2 (62)
Remark The term τ is the driving force or the configurational force acting on the interface Γ. One can notice that the configurational force is always positive whatever the location of the interface.

The exact location of the interface is the one that verify the following optimality criteria. Indeed, if the configurational equilibrium is satsified, then p is continous accross Γ and p ≥ 0. Moreover, Trε is also continous accross Γ and Trε ≤ α. In this case, the couple (u, p) verify the equilibrium conditions, the boundary conditions, and the Karush-Kuhn-Tucker conditions of optimality [START_REF]La méthode des contraintes actives appliquée au contact unilatéral[END_REF]. Hence, (u, p) is the solution of the problem (4). The final algorithm, therefore will consist in finding the position of Γ such as the configurational equilibrium defined in Table 3 is ensured . Then, we verify the conditions of optimality exposed in Table 3. Around the points where these two conditions are violated, a new trial boundary needs to be introduced . We then evolve once again the whole new interface so as to verify the configurational equilibrium.

Ω c Γ Ω f configurational equilibrium div P + ∇ f .u = 0 τ = 0 div P + ∇ f .u = 0 optimality conditions p c ≥ 0 Trε f -α ≤ 0

Analytical solution on a one-dimensional example

Let us consider the simple one dimensional problem in the small strain elastostatic setting with uniaxial state of strain as shown in Figure 4. We wish to limit the uniaxial strain of a truss by a value α under a linear prescribed load f (x) = f x.

We split the bar into two domains. The constrained domain [0, y] and the unconstrained domain [y, L]. The exact displacement, strain, and pressure field expression along the bar for any given location y of the interface are given in Appendix B.

We choose α such that the exact location of the interface y e is at the middle of the bar.

α = 3 f 8(λ + 2µ) L 2 (63) Denoting ε = ε(λ+2µ) f L 2
and σ = σ f L 2 , the exact strain and stress fields ε and σ are depicted in Figure 6.

The evolution of the potential energy and the configurational force with the location of the interface are depicted in Figure 5 and the analytical expression are given in Appendix B. When the interface is not located at the exact position, the strain field is discontinuous at the interface. On can notice that the potential energy of the system is stationary, when the interface is located at ξ = 1 2 . It means that at the exact location of the interface, the configurational force is zero, as one can see in Figure 5. Finally, one can see that the configurational force is always positive whatever the location of the interface and is minimum at the exact location. Hence, the total potential energy is always increasing with the location of the interface ξ. The exact location makes the potential energy stationary but not minimum as one can see in Figure 5.

Φ (λ+2µ) f 2 L 5 ξ 1 0.5 0 -0.06 -0.065 ∂Φ ∂ξ (λ+2µ) f 2 L 5

Shape sensitivity

The goal of this section is to establish the change in (u, p) as a result of the migration of the interface. Consider the perturbed constrained zone Ω τ c and the perturbed unconstrained zone Ω τ f . Let us introduce a velocity field θ responsible for the migration of the interface and the following transformation:

Ω c → Ω τ c x τ = x + τθ(x) ∀x ∈ Ω c Ω f → Ω τ f x τ = x + τθ(x) ∀x ∈ Ω f
The shape change of each phase is expressed using a sufficiently smooth velocity field θ and a parameter τ. The trick introduced in [START_REF]Shape sensitivity analysis in linear elastic fracture mechanics[END_REF] is to use the classical continuum terminology and assimilate the change of shape to a motion of a body. We rewrite the state equation (50) with the parameter τ.

a(τ, u, v) + b(τ, p, v) = l(τ, v) ∀v ∈ V 0 (64) b(τ, u, q) = p(τ, q) ∀q ∈ Q (65)
Taking the derivative of the above equation with respect to τ

a(τ, ∂u ∂τ , v) + b(τ, ∂p ∂τ , v) = ∂l(τ, v) ∂τ - ∂a(τ, u, v) ∂τ - ∂b(τ, u, q) ∂τ ∀v ∈ V 0 b(τ, ∂u ∂τ , q) = ∂p(τ, q) ∂τ - ∂b(τ, u, q) ∂τ ∀q ∈ Q (66)
Defining ∂u ∂τ = u and ∂p ∂τ = ṗ, the above formulation becomes Find ( u, ṗ) such that:

Ω ε( u) : C : ε(v)dΩ + Ω c ṗTrε(v)dΩ = Ω c p ∇θ T : ∇vdΩ - Ω c p div(θ)Tr(∇v)dΩ + + Γ T t d (-n.∇θ T n + div(θ))vdΓ + + Ω f (x) div(θ).vdΩ + Ω (∇ f θ).vdΩ + + Ω T : ∇vdΩ ∀v ∈ V 0 (67) with T = [C : (∇u∇θ) s +C : ∇ s u∇θ T -C : ∇ s u div(θ)] Ω c Trε( u)qdΩ = Ω c
Tr((∇u∇θ) s )qdΩ -

Ω c Tr(∇ s u) div(θ)qdΩ + + Ω c α div(θ)qdΩ ∀q ∈ Q (68)
We can notice that the two bilinear forms of (66) are similar to the bilinear forms of (50). Only the linear forms on the right hand side are different and they imply the knowledge of the fields (u, p) and the velocity field θ. This weak form can be discretized as the previous one.

If we use the same approximation for the two couple of fields ( u, ṗ) and (u, p), we will obtain two linear systems to solve with the same tangent stiffness matrix. From a numerical point of view, this remark has a huge importance, since we will not have to compute twice the tangent matrix. If we use a LU solver for the two linear systems, we then make the LU decomposition just once.

Numerical Strategy

The goal of this section is to present the numerical tools and the numercial procedure to move the interface Γ towards its exact location.

Level-set/X-FEM

To represent the surface of discontinuity separating the domain Ω c et Ω f we used the notion of level-set [START_REF]Level Set Methods and Fast Marching Methods : Evolving Interfaces in Computational Geometry[END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed : algorithms based on hamilton-jacobi formulations[END_REF]. The level-set method permits an implicit representation of surfaces through a higher dimensional function φ(x). A level-set is a function that evaluates at each point of the domain the signed distance to the surface of discontinuity. The surface Γ is then defined by

Γ = {x | φ(x) = 0} (69)
By convention, the sign of the level-set is negative if the point is located in Ω c and positive if the point is located in

Ω f . φ(x) < 0 if x ∈ Ω c (70) φ(x) > 0 if x ∈ Ω f (71)
The iso-zero of the level-set locates the true location of the interface. The level-set is evaluated at the nodes of the mesh and we interpolate these values with the classical linear interpolation functions N i .

φ(x) = ∑ i N i (x)φ i (72)
Because of the presence of the surface of discontinuity Γ, the strain tensor field presents a discontinuity. To avoid poor rates of convergence we enrich the classical approximation through a partition of unity technique for strain jump [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]. The classical discretization (34) becomes

u h | Ω e = N ∑ i=1 u h i N u i + N ′ ∑ j=1 N u j F(x)b j (73)
Additional degrees of freedom b j have been added at the nodes for which the support is split by the interface. The enrichment function F, is chosen such that its gradient is discontinuous along Γ. Thus, we can represent the discontinuity in strain. The function F is defined as [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF]:

F(x) = ∑ i |φ i |N i (x) -| ∑ i φ i N i (x)| (74)
The function F is depicted in Figure 7, for a one dimensional case.

In our case, we enrich the linear part of the approximation. It means that the N u j functions are the classical linear interpolation function even if the N u i are quadratic. The pressure field p is only defined over the domain Ω c . To take into account the presence of Γ, the weak form involving p is integrated only on Ω c . Thanks to the level-set notion, the interface Γ do not need to be meshed. We just have to take care of the integration step as desribed in [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite element method[END_REF].

Front evolution

We have shown in section 4.1 that the condition in Table 3 ensures that (u, p) is the solution of ( 4). Numerically, we will first find the location of the interface that ensures the condition τ = 0. When this location is found, we verify the two other one.

Thanks to (62), the condition τ = 0 on Γ is equivalent to the condition p = 0 on Γ. Rewriting this condition as a variational form, we get

Γ pq * dΓ = 0 ∀q * ∈ L 2 (Γ) (75)
We now take the directional derivative of (75) for a configurational change v, where v is a normal velocity responsible for the migration of the interface Γ.

D Γ pq * dΓ [v] = Γ Dp[v]q * + Γ pq * DdΓ[v] ( 76
)
Since

DdΓ[v] = ∇ sur .v = ϑv.n (77)
where ϑ denotes the curvature of the interface and ∇ sur .v is the surfacic divergence of v We thus get

D Γ pq * dΓ [v] = Γ Dp[v]q * dΓ + Γ ϑpv.nq * dΓ (78)
The algorithmic scheme is the following. We first initialize Γ = Γ 0 . And then at iteration k, we calculate the velocity v k using the following Newton-Raphson procedure

Γ k pq * dΓ + D Γ k pq * dΓ [v k ] = 0 (79) Γ k pq * dΓ + Γ k Dp[v k ]q * dΓ + Γ k ϑpv k .nq * dΓ = 0 (80)
If we use a modal representation for the velocity field v k and the test functions

q * v k = ∑ i F i ṽi n (81) q * = ∑ i F i q * i ( 82
)
The modes F i are functions defined over the interface Γ. To do so, Γ is parametrized using a curvilinear abscissa s ∈ [0, 2π]. Two types of modes are here introduced to describe the field.

• cosinus modes

F i = cos(is) (83)
• sinus modes

F i = sin(is) (84)
If we neglect the curvature term Γ k ϑpv k .ndΓ that becomes small as p tends to zero, we get

Γ k Dp[v k ]q * dΓ = - Γ k pq * dΓ (85)
that leads to the following matricial system

[S][ ṽ] = [P] (86)
where [S] is the matrix containing the sensitivity of p to the migration of the interface for each modes F i n, and [P] is the vector containing the coefficients of the modal projection of the pressure field p on the interface Γ.

We now present the procedure to calculate the matrix [S]. For example, we want to determine the sensitivity of p to the migration of the interface for the mode i . This field is denoted Dp[F i n].

To calculate this field using the formulation developed in section 4.3, we first need to extend the velocity F i n to the whole domain Ω.

Let us define the extended velocity θ i as

θ i = F i (s(x))∇φ(x)r(x) (87)
Where the gradient of the level-set ∇φ is used to make the velocity normal to the interface. Since the level-set function is a signed distance function, it must be of norm 1. The last function r(x) is a function which permits to extend the velocity Ω. This function is only non zero in a thick band, where its value is 1. This band can be defined as a layer of elements around the interface. The Figure 8 shows the band around a circular interface.

The size of the band determines the size on which the right side of ( 67) and (68) need to be computed. The optimal size of the band with regard to numerical accuracy will be the subject of further work. But, in all numerical experiments shown in this study, the band has a thickness of approximately 4 elements. As defined above, s is a curvilinear coordinate along the interface Γ. This curvilinear coordinate needs to be extended in Ω for equation 87 to make sense. One possibility is to compute s at each node v of the mesh by taking the value of s on the closest node to v on the interface. This algorithm is very similar to some implementation of the velocity extension algorithm very common in level-set method, where a field only known on the iso-zero of the level-set need to be extended every where else [START_REF]Level Set Methods and Fast Marching Methods : Evolving Interfaces in Computational Geometry[END_REF].

Once the field Dp[θ i ] is known, we determine the modal projection of this field on the interface Γ. The coefficients found constitutes the vector i of the matrix [S]. Once the velocity field v k is found, the final step is to propagate the level-set by solving the following Hamilton-Jacobi equation. [START_REF]Non planar 3d crack growth by the extended finite element and level-sets. part ii : level set update[END_REF] To solve it, we use the triangulated versions of the levelset methods as presented in [START_REF]Level Set Methods and Fast Marching Methods : Evolving Interfaces in Computational Geometry[END_REF], more preciselly, the so called explicit positive coefficient scheme.

∂φ(x,t * ) ∂t * + (v k .n) ∇φ(x,t * ) = 0 (88)
Where t * designates a fictious time.

Once the level-set is propagated, we have to renormalize the level-set using the same type of equation.

Convergence criterion After each propagation step, and once the level-set is normalized, the equilibrium (50) is solved. Then, the following criteria is computed

c p = Γ p 2 dΓ 1 2 ≤ ε ( 89
)
where ε is a tolerance that can be set arbitrarily.

If the convergence criterion is not verfied, a new interface velocity is computed, and the levelset is once again propagated.

The whole algorithm, composed of two overlaped loops, is detailed below

Discretization

To enforce the constraint imposed on Ω c , we have used a Lagrange multiplier approach. We have to discretize both the displacement and the Lagrange multiplier fields. As explained in the previous part, we have to use specific mixed elements that satisfy the LBB condition. Unfortunately, the ridge function introduced in the X-FEM enrichment modify the finite dimensions subspace V h . However, it is shown in [START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF] that the stability of enriched mixed element is preserved through the numerical inf-sup test. Where the pressure is defined (all elements that have at least one node inside Ω c ), the element used is the P2/P1 cont (see Figure 9). This element is composed of a quadratic displacement and a linear continous pressure. Here, the use of continuous pressure element is possible unlike with the classical approach where this type of elements causes problems of convergence of the active-set algorithm. Other element are simply classical quadratic lagrange element for the displacement field. Each component of the displacement field is further enriched, at all the node ve that have there support cut by Γ, by the function F defined previously. This formulation therefore give a total number of degree of freedom ndo f = 2 * nv + 1 * nvp + 2 * ne + 2 * nve where nv is the total number of nodes, ne the total number of edges, nvp the total number of node where a pressure is defined, and nve, the total number of node where the displacement field is enriched.

Numerical examples

Convergence studies

Two cases will be studied here. A one-dimensional case and a cylinder case for which the analytical solution is known. The analytical solution will be used to determine the h-convergence of each method. To compare the exact and the numerical solution, we use the L 2 norm for the pressure field (90) and the energy norm for the displacement (91).

e p = Ω (p h -p ex ) 2 dΩ 1 2
( Ω p 2 ex dΩ)

1 2
(90)

e e = Ω (ε D h -ε D ex ) : C : (ε D h -ε D ex )dΩ 1 2 Ω ε D ex : C : ε D ex dΩ 1 2 (91)
Algorithm 2 Algorithm of front evolution 1: Compute (u, p) using (50) 2: while Conditions of optimality in Table 3 are not verified do 3:

Detect points where these conditions are violated 4:

Initialize an interface around these points Compute Dp[q i ] for each mode 11:

Compute the front velocity using (

Update Level-Set 13:

end loop 14: end while

One-dimensional case

We apply now the new approach to solve the simple one-dimensional problem in the small strain elastostatic setting with uniaxial state of strain described in section 4.2. We wish to limit the uniaxial strain of a truss by an arbitrary value α under a prescribed load f (x) = f x. We set α = 3 8 , (λ + 2µ) = 1, f = 1 and L = 1.

The analytical solution of this problem is given in Appendix B. The results of convergence are presented in Figure 10. With the classical approach, the L2 norm error on the pressure and the energy norm error on the displacement exhibits a rate of convergence in O(h ). The rate of convergence is degrated because of the presence of a singularity in the exact solution as we can see in Figure 6. With the new approach, thanks to the X-FEM enrichment, we can catch the singularity and we recover the optimal rate of convergence in O(h 2 ).

Two-dimensional case

The second case studied here is an axisymmetric problem. It is presented in Figure 11. We consider a cylinder subjected to a volumetric force f = u r , where u r is the unit radial vector. The cylinder is held fixed at radius a = 0.5 and the external radius is b = 1. We choose λ = 0 and µ = 1 2 for the Lamé coefficients. We want to limit the value of Tr(ε) by α = 0.5. The analytical solution is given in Appendix C.

The results of convergence are presented in Figure 12. With the classical approach, the L 2 norm error on the pressure and the energy norm error on the displacement exhibits a rate of convergence in O(h presence of a low regularity in the exact solution. With the new approach, thanks to the X-FEM enrichment, we can catch the low regularity and we recover the optimal rate of convergence in O(h 2 ) for the L 2 norm error on the pressure. Unfortunately, we do not recover the optimal rate of convergence in O(h 2 ) for the energy norm error on the displacement. It comes from the fact that the interface is discretized with straight edges while the displacement is quadratic. We could improve the convergence using the method described in [START_REF] Dréau | Studied X-FEM enrichment to handle material interfaces with higher order finite element[END_REF].

A more complex case and related numerical experiments.

We want here to apply the new approach to a more complex two-dimensional case as shown in Figure 13. The example studied here is a zero centered 2D rectangular plate subjected to a body force f (x) and held fixed at the bottom. The height b of the plate is set to 1 and the width a is set to 0.5. We choose λ = 0.9 and µ = 0.6 for the Lamé coefficients.

The following body force is chosen :

f (x, y, z) = xu x + yu y ( 92
)
Figure 13: A more complex case This kind of body force permits us to have a non homogeneous repartition of Trε. Without any constraint, the equilibrium solution can be computed and in this case the iso-values of Trε are ellipses. With the constraint on, we do not know of an analytical solution for this example. It would be too complex to calculate the exact shape of Ω c , but this example exhibits an interesting repartition of Trε such that the numerical research of the shape of Ω c involves several modes.

The starting position of the interface Γ is chosen circular and zero-centered. The initial radius of the cylinder is fixed arbitrarily at r = 0.1. Only four modes are used to describe the shape changes. After four iterations, the convergence criteria is reached and the interface is close to an ellipse. It can be compared to the constrained zone as computed by the classical algorithm, by displaying the values of χ. On figure 14, left, we therefore show a representation of χ, only evaluated at gauss point during the computation, as a linear by element field. On top of this field we display the iso-zero value of the level-set representing the interface as obtained with the new approach. As can be seen on the figure, both method give similar results in term of constrained zone. The new approach has the advantage of defining the zone in a more explicit sense.

The problem we are solving here are for now too small to do a proper efficency comparaison of the two methods: the solving time is below one minute on regular desktop for both method. The new method is probably, for equal meshes, more expensive at each iteration than the classical one. Indeed The sensibilities for each velocity mode need to be computed, and then an Hamilton Jacoby equation need to be solved. This last comment can still be nuanced a bit since the pressure field needed to enforce the equality constrain is only defined inside the constrained zone, delimited by the zero value of the level-set, and the Hamilton Jacoby equation could be solved in a narrow band only. Over all, the most expensive computations are the assembly and the factorisation of the linear system needed to solve the elasticity problem. For the mesh we used for the present test case, this cost is equivalent in both method, since the total number of degree of freedom is 11576 for the classical method and reach 10689 at convergence for the new method.

The advantages of the new method is that it is able to give more precise results on a given mesh than the classical one, as shown by the convergence rate in previous section. It provides a clear answer to the question of the position of the interface, which enable the enrichment its enrichment to better capture the irregularities of the fields and eventually the description of more involved physics at this interface.

In order to stress the presented algorithm a little bit more, we tried to solve the same problem with different starting position of the interface. Figure 14, middle and right, show the evolution of the constrained zone during the iterations. The middle picture show the evolution when the starting interface is an off center circle. The right picture show the evolution when the starting interface is made of two circles. In both cases the algorithm converge to the same final shape as the one obtain with the previously defeined starting interface. To obtain the correct results, more modes are needed in this two cases than for the first one. In fact this tests permitted to discovers that the numbers of mode used are optimally set when it is proportional to l, the total lenght of the interface divided by h, a characteristic size of the finite elements. Indeed, when the interface lenght is very small, there is no point to compute a large number of modes, since those modes are going to be very badly represented on the mesh. When the relative lenght of the interface become larger, more mode can be used. In the experiment shown here, we used a number of mode n such as n is the closest integer to l 2h . Our numerical experiments also shown that, when during the iteration the local curvature of the interface become to high compared to the element size, the algorithm may take a very long time or even fail to converge. The case when we started with the interface made of two circle is one such a case. Indeed, when the two circles merge, there is a point where the local curvature could even become singular. This problem might be due to the fact that we neglicted the curvature term in equation 80. To avoid computing the curvature, we used the following fix : as long as the residual in the newton method is big enough, the level-set function is smoothed prior to compute the new advance. We used a very simple smoothing the level-set : at each node v, the value of the level-set is replaced by the result of a least-square fit of the level-set values at each nodes connected to the current node v. In black, χ = 1, in grey χ = 0, for the classical method. The interface between constrained and unconstrained zone for the new method is the white line. Middle and Right : evolution of the interface during iterations for the new approach. Middle: the starting interface is one off centered circle. Right: the starting interface is made of two circle.

Conclusion

We have presented in this study a new approach to deal with continuum media subjected to an inequality constraint. The main difficulty of inequality constrained minimization problems is that the shape of the constrained zone is a priori unknown. Hence, the strategy proposed in this paper is to fix arbitrarily an initial constrained zone and to propagate it correctly until the correct shape is found following a specific criteria.

After a theoretical presentation of the problem and the classical approach, we have developed an algorithmic scheme that ensures to find the exact location of the interface. This new approach, relies on the notion of level-set , which is a powerful tool to represent moving interfaces in a domain, coupled with a X-FEM enrichment. Several numerical examples have been performed and exhibit a significant gain in convergence compared to the classical approach. Indeed the exact representation of the constrained zone avoids the loss of convergence due to the lack of regularity of the exact solution.

It must be made clear that proposed approach has some limitations related to the discretization. If the active zone is too small (or too thin), for the mesh size, it will not be detected. The detection will only be possible if the loading increases and the zone enlarges. If the knowledge of the loading for the active zone initiation is important, adaptivity could be introduced.

Finally, the new use of level sets to handle inequality contraints (or ILS in short for Inequality Level Set) introduced in this paper could be applied to other type of constraints such as inextensibility constraints in hyperelastic materials.

p(x) = p f (x) if x ∈ Ω f p c (x) if x ∈ Ω c (94)
In the same way, we denote by ε c the strain field on the domain Ω c and by ε f the strain field on the domain Ω f and we have

ε(x) = ε f (x) if x ∈ Ω f ε c (x) if x ∈ Ω c ( 95
)
ψ f -ψ c = µ(ε c : ε c -ε f : ε f ) + λ 2 (α 2 -(Trε f ) 2 ) (96) = µ(ε c + ε f ) : (ε c -ε 1 ) + λ 2 (α + Trε f )(α -Trε f ) (97) = (µ(ε c + ε f ) + λ 2 (α + Trε f )I) : (ε c -ε f ) (98) = ( σ l + σ c 2 - p c 2 I) : (ε c -ε f ) (99)
We thus get

τ = n[[P]]n = ( σ l + σ c 2 ) : (ε c -ε f ) -(σ.n).((∇u c -∇u f ).n) (100)

Appendix B

In this appendix, we detail the analytical solution for the 1D problem

u(x) = αx if 0 ≤ x ≤ y f 6(λ+2µ) (y 3 -x 3 ) + f 2(λ+2µ) L 2 (x -y) + αy if L ≥ x ≥ y (105) ε(x) = α if 0 ≤ x ≤ y f 2(λ+2µ) (L 2 -x 2 ) if L ≥ x ≥ y (106) p(x) = f 2 (L 2 -x 2 ) -(λ + 2µ)α 0 ≤ x ≤ y 0 if L ≥ x ≥ y ( 107
)
Let us now compute the potential energy Φ in the system for any given location y of the interface.

C 2 = - 2 3 + C 1 2 (120)
The constant r e is the zero of the equation g(y) = 0, with g g(y) = -

2 3 y + C 1 2 - C 2 y 2 - α 2 + C y 2 (121)
Bibliography

Figure 1 :

 1 Figure 1: Problem statement Let use define the space of kinematically admissible displacements

Figure 2 :

 2 Figure 2: P2 + /P1 disc mixed element

Figure 3 :

 3 Figure 3: Problem statement 2

Figure 4 :

 4 Figure 4: A one-dimensional example

Figure 5 :Figure 6 :

 56 Figure 5: Evolution of the potentiel energy (left) and its derivative (right) with ξ

Figure 7 :

 7 Figure 7: The enrichment function

Figure 8 :

 8 Figure 8: The band of elements around the level-set on the left, and the normal velocity on the right

Figure 9 :

 9 Figure 9: P2/P1 cont mixed element

Figure 10 :

 10 Figure 10: Convergence of the L 2 norm error on the pressure (left), and the energy norm error (right)

Figure 12 :

 12 Figure 12: Convergence of the L 2 norm error on the pressure (left), and the energy norm error (right)

Figure 14 :

 14 Figure 14: Left : Comparaison between the classical and new the method. In black, χ = 1, in grey χ = 0, for the classical method. The interface between constrained and unconstrained zone for the new method is the white line. Middle and Right : evolution of the interface during iterations for the new approach. Middle: the starting interface is one off centered circle. Right: the starting interface is made of two circle.

2 (

 2 Trε f -α) 2 (103)On the intertface, we have due to the continuity of the normal stress vectorp = (κ + 3 4 µ)(Trε f -α) = (λ + 2µ)(Trε f -α)

Table 2 :

 2 The Γ problem

Table 3 :

 3 The Γ optimality

). As in the one-dimensional case, the rate of convergence is degraded because of the Figure 11: Axisymmetric problem

Appendix A

In this appendix, we demonstrate the equation (62). We consider the velocity θ normal to the interface. If we take into account the continuity of the displacement field and normal stress accross the interface Γ and assume the volumetric force f to be continous accross Γ, we get

We will next denote by p c the pressure field on the domain Ω c and by p f the pressure field on the domain Ω f . We then have

Denoting ξ = y L , and after calculus, we obtain

where τ(y) = [[P]] y is the configurational force at the interface y.

Appendix C

In this appendix, we detail the analytical solution for the axisymmetric problem

where r designates the radial coordinate. And D is the a constant defined by