Marie Rebouah

Guilherme Machado

Grégory Chagnon
email: gregory.chagnon@grenoble-inp.fr

Denis Favier

Anisotropic Mullins stress softening of a deformed silicone holey plate

Keywords: Mullins effect, stress-softening, finite element implementation, constitutive equation

Introduction

Both natural and synthetic elastomers are widely used in industrial design. Numerical simulations are often used to develop new parts. These simulations are fundamental in the optimization process of parts design. One of the main difficulties of the engineers is the choice of an adapted constitutive equation able to represent the behavior of the rubber like materials. The choice is often oriented towards an hyperelastic constitutive equation requiring an energy density. Finite element codes present a large choice of strain energy densities forms, even if the [START_REF] Mooney | A theory of large elastic deformation[END_REF] constitutive equation stays one of the most used.

However rubber like materials present a lot of other phenomena than pure hyperelasticity to take into account in their behavior, as for example: the Mullins effect, the hysteretis and time dependent behavior. The Mullins effect can be very important as the mechanical behavior of the material can totally change after a first loading. For some rubber parts, it is very important to take into account the Mullins effect, but very few constitutive equations are implemented in industrial finite element codes. Eventhough, many researchers have developed isotropic constitutive equations for the Mullins effect and proposed a finite element implementation (see for example [START_REF] Cantournet | Mullins effect and cyclic stress softening of filled elastomers by internal sliding and friction thermodynamics model[END_REF]; [START_REF] Chagnon | Development of new constitutive equations for mullins effect in rubber using the network alteration theory[END_REF]; [START_REF] Gracia | A comparison between pseudo-elastic and damage models for modelling the Mullins effect in industrial rubber components[END_REF]; [START_REF] Guo | A composites-based hyperelastic constitutive model for soft tissue with application to the human annulus fibrosus[END_REF]).

Another important point is that, many new experimental data are proposed in the literature to emphasize that the Mullins effect is strongly anisotropic [START_REF] Muhr | Experimental determination of model for liquid silicone rubber: Hyperelasticity and Mullins effect[END_REF][START_REF] Park | Anisotropy in gum and black filled SBR and NR vulcanizates due to large deformation[END_REF][START_REF] Pawelski | Softening behavior of elastomeric media after loading in changing directions. Constitutive models for rubber[END_REF][START_REF] Besdo | Properties of rubber like materials under large deformations explained by self-organizing linkage patterns[END_REF][START_REF] Laraba-Abbes | A new Taylor-made methodology for the mechanical behavior analysisof rubber like materials: II. Application of the hyperelastic behavior characterization of a carbon-black filled natural rubber vulcanizate[END_REF]Diani et al., 2006a;[START_REF] Hanson | Stress softening experiments in silica-filled polydimethylsiloxane provide insight into a mechanism for the Mullins effect[END_REF][START_REF] Itskov | Experimental observation of the deformation induced anisotropy of the Mullins effect in rubber[END_REF][START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF][START_REF] Dorfmann | A constitutive model for the Mullins effect with changes in material symmetry[END_REF]. Different constitutive equations have been proposed but they are not adapted to finite element implementation. The only formulation implemented in a finite element code was proposed by [START_REF] Göktepe | A micro-macro approach to rubber-like materials. Part III: The micro-sphere model of anisotropic Mullins-type damage[END_REF] who extended the approach of [START_REF] Miehe | A micro-macro approach to rubber-like materials -part i: the nonaffine-micro-sphere model of rubber elasticity[END_REF]. This model presents good results, but its formulation makes difficult its numerical use into finite element codes, as it needs complex algorithm implementation and also some linearizations (see for example [START_REF] Miehe | Superimposed finite elastic-viscoelastic-plastoelastic stress response with damage in filled rubbery polymers. experiments, modelling and algorithmicimplementation[END_REF]; [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF]). A formulation written with strain invariants permits to avoid such technics and ensure a better convergence. Some constitutive equations has been developed for living tissues [START_REF] Peña | On the Mullins effect and hysteresis of fibered biological materials: A comparison between continuous and discontinuous damage models[END_REF][START_REF] Bose | Computational aspects of a pseudo-elastic constitutive model for muscle properties in a soft-bodied arthropod[END_REF][START_REF] Kroon | A new constitutive model for multilayered collagenous tissues[END_REF] but they are limited to materials presenting two reinforced directions. In this way, in this paper, a new anisotropic model based on strain invariant formulation, is proposed and implemented in a finite element code.

Recently, [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF][START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF]) developed a large database for a filled silicone rubber including on one hand cyclic classical experimental tests and on the other hand uniaxial tests realized after different uniaxial and biaxial tension tests. This database is, here, used to build a new constitutive equation easily implementable in finite element codes. In Section 2, the constitutive equation is detailled and the ability of the model to describe recent experimental data is discussed. In Section 3, the subject of the numerical implementation of the model is treated. In Section 4, the ability of the model to describe complex structures is tackled by means of a specific test on a rectangular plate containing five holes. The global and local estimations of the model are compared to experimental measures. Finally, Section 5 contains some concluding remarks of the modeling.

Anisotropic modeling of the Mullins effect

General formulation in strain invariants

Different anisotropic approaches to model Mullins effect were proposed in literature, but none of them was only expressed in term of strain invariants. [START_REF] Shariff | An anisotropic model of the Mullins effect[END_REF] and [START_REF] Itskov | A thermodynamically consistent phenomenological model of the anisotropic Mullins effect[END_REF] took into account three principal damage directions to reproduce a special behavior in the direction orthogonal to loading. In a more general way, the spatial repartition of [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF] was used by many authors to create an anisotropic model. [START_REF] Diani | Observation and modeling of the anisotropic visco-hyperelastic behavior of a rubberlike material[END_REF] and [START_REF] Dargazany | A network evolution model for the anisotropic Mullins effect in carbon black filled rubbers[END_REF] generalized the network evolution proposed by [START_REF] Marckmann | A theory of network alteration for the Mullins effect[END_REF] to an anisotropic approach by taking into account the maximum elongation in each spatial direction. Later [START_REF] Merckel | A simple framework for full-network hyperelasticity and anisotropic damage[END_REF][START_REF] Merckel | A Mullins softening criterion for general loading conditions[END_REF] introduced a new framework and proposed a softening anisotropic criterion adapted to complex loading states.

The stress softening phenomenon has been often associated to the presence of fillers in the rubber, but [START_REF] Harwood | Stress softening in natural rubber vulcanizates. Part 2. stress softeningeffects in pure gum and filler loaded rubbers[END_REF] showed that stress softening can also occur in unfilled rubber, even if it is reduced compared to filled rubber. For silicone rubbers, [START_REF] Meunier | Mechanical experimental characterisation and numerical modelling of an unfilled silicone rubber[END_REF] observed no Mullins effect for an unfilled one, whereas [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF][START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF] observed stress softening for a filled one. As a consequence, it can be considered that fillers in silicone rubbers are mainly responsible of the Mullins effect. Thus, as proposed by [START_REF] Govindjee | Mullins' effect and the strain amplitude dependence of the storage modulus[END_REF] the strain energy density W is additively decomposed into two parts: one that represents the energy density of the chains linked to other chains W cc and an other part that represents the energy density of the chains linked to filler W cf , the total strain energy density is W = W cc + W cf . It is considered that only (W cf ) can evolve with the Mullins effect. As a consequence W cc is represented by a classical hyperelastic isotropic energy density and W cf must be represented by an anisotropic strain energy that can evolve with the deformation history of the material. The ideal representation would be to propose a full integration of all spatial directions as proposed by [START_REF] Wu | On improved network models for rubber elasticity and their applications to orientation hardening in glassy polymers[END_REF] in hyperleasticity, but it is not adapted to finite element implementation. A spatial discretization is needed. Forty-two initial spatial directions, noted A (i) , are introduced, these directions are those proposed by [START_REF] Bazant | Efficient numerical integration on the surface of a sphere[END_REF]. Therefore, the strain energy density is written as:

W = W cc (I 1 , I 2 ) + n ∑ i=1 ω (i) F (i) W (i) cf (I (i) 4 ) (1)
where I 1 , I 2 are the first and second strain invariants of the right Cauchy-Green strain tensor C. The strain in each direction is defined by means of i) .CA (i) . ω (i) represents the weight of each direction and F (i) is the Mullins effect evolution function. The initial direction A (i) are transformed in a (i) by a (i) = FA (i) , where F is the deformation gradient.

I (i) 4 = A (
Classicaly, in an isotropic approach, the evolution function F (i) would be written through the strain energy density, but [START_REF] Chagnon | On the relevance of continuum damage mechanics as applied to the Mullins effect: theory, experiments and numerical implementation[END_REF] showed that the first invariant can also be used. In an anisotropic approach, the elongation in each direction is used [START_REF] Diani | Observation and modeling of the anisotropic visco-hyperelastic behavior of a rubberlike material[END_REF], knowing that the elongation is the square root of the invariant I 4 . According to the conclusions of [START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF], it is chosen to describe the stress-softening function according to I 1 and I 4 (i) . For each direction (i), an evolution function which depends on three terms F 1 , F 2 , F 3 is proposed:

F = 1 -F 1 (I 1 max -I 1 )F 2 ( I 4 max(i) -I 4 (i) ) F 3 ( I 4 max(i) I 4 max ) (2)
Where I 1 max and I 4 max(i) represent the maximum values taken during the material history by I 1 and I 4 (i) respectively. I 4 max = max i (I 4 (i) ) is the maximum dilatation in space and time. As proposed by [START_REF] Zuñiga | A new phenomenological model for stress-softening in elastomers[END_REF], a function that is constant during first loading and that evolves with the maximum and current deformations is imposed for the evolution function.

A particular form for evolution function

A [START_REF] Mooney | A theory of large elastic deformation[END_REF] constitutive equation is chosen for (W cc ), and a [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF] quadratic equation K(I

(i) 4 -1) 2 is chosen for W ( cf i).
A first particular form is proposed for the stress-softening function, considering that a minimum of parameters should be introduced:

F (i) = 1 -η √ I 1max -I 1 I 1max -3 ( I (i) 4 max -I (i) 4 I (i) 4 max -1 ) ( I (i) 4 max I 4 max ) 4
(3)

In this way, the evolution function depends only on one parameter: η. The large experimental database proposed by [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF][START_REF] Machado | Induced anisotropy by the mullins effect in filled silicone rubber[END_REF]) on a filled silicone rubber is used to fit the model. These experimental results are decomposed into three parts: first the classical uniaxial tension, planar tension and biaxial tests realized by means of a bulge test (Machado et al., 2012a), second the complex tensile tests with change of directions after the first loading; and third biaxial tensile tests followed by uniaxial tensile tests. The three hyperelastic parameters C 1 , C 2 and K (i) are obtained by fitting the different first loading curves. The followings values are obtained:

C 1 = 0.05MPa, C 2 = 0.03MPa and ∀i K (i) = 0.10MPa. The last parameter is fitted to represent every second loading curves: η = 4. The simulations of the cyclic uniaxial tensile, pure shear and equibiaxial tensile tests are presented in Fig. 1. It appears that the model describes well the stress softening for all these tests. The model predicts well uniaxial and pure shear tests whereas equibiaxial first loading curve is underestimated. This phenomenon is expected since first loading depends only on the hyperelastic equation. As explained by [START_REF] Marckmann | Comparison of hyperelastic models for rubber-like materials[END_REF], it is difficult to fit all the different tests with the same energy density. The proposed model is now compared with the experimental data of the two complex preconditioning methods. First, Fig. 2, presents the results for tensile tests with a change of loading direction between the first and second loadings. The results from the model do not superimpose exactly experimental data, but all trends are quite well described for the different directions. Second, the model predictions are compared with biaxial pre-stretching tests results. The biaxial loading is characterized by the biaxiality ratio defined as µ = ln(λ min )/ln(λ max ) (where λ min and λ max are the minimum and maximum in-plane principal elongation). Tests with different biaxiality ratios were used for the simulation. The comparison of the second loading curves is presented in Fig. 3. It appears that the stress softening is slightly overestimated but the return point on the first loading curve is very well described.

Finite element implementation

Considering the model proposed in the last section, the strain energy density depends on the first, second and forty-two fourth invariants, i.e., W = W(I 1 , I 2 , I (i) 4 ). As the numerical imple- Curve a presents simulation of the model for the second load after an equibiaxial test; curve g presents simulation of the model for the second load after a biaxial test of biaxiality rate µ = 0.7; curve h presents simulation of the model for the second load after a biaxial test of biaxiality rate µ = 0.5 mentation needs a quasi-incompressible formulation [START_REF] Federico | An energetic approach to the analysis of anisotropic hyperelastic materials[END_REF][START_REF] Holzapfel | A new constitutive framework for arterial wall mechanics and a comparative study of material models[END_REF][START_REF] Holzapfel | A viscoelastic model for fiber-reinforced composites at finite strains: continuum basis, computational aspects and applications[END_REF][START_REF] Ogden | Nonlinear elasticity, anisotropy, material stability and residual stresses in soft tissue[END_REF][START_REF] Ogden | Non-linear Elastic Deformations[END_REF][START_REF] Sun | Finite element implementation of a generalized Fung-elastic constitutive model for planar soft tissues[END_REF], the equation must be rewritten in terms of incompressible invariants as follows:

W = W( Ī1 , Ī2 , Ī4 (i) ) + U(J) (4)
with J 2 = I 3 = det(C) is the volume variation, and Ī1 = I

-1/3 3 I 1 , Ī2 = I -2/3 3 I 2 , Ī4 (i) = I -1/3 3 I (i) 4
are the incompressible invariants. The problem of the implementation of anisotropic hyperelastic energy densities has ever been tackled in literature, see for example [START_REF] Peña | On the Mullins effect and hysteresis of fibered biological materials: A comparison between continuous and discontinuous damage models[END_REF][START_REF] Bose | Computational aspects of a pseudo-elastic constitutive model for muscle properties in a soft-bodied arthropod[END_REF][START_REF] Kroon | A new constitutive model for multilayered collagenous tissues[END_REF]. After the decomposition given in Eq. 4, the second Piola-Kirchhoff stress tensor S is obtained by the strain energy derivation, as

S = S + S vol = 2 { ∂ W ∂ Ī1 ∂ Ī1 ∂C + ∂ W ∂ Ī2 ∂ Ī2 ∂C + n ∑ i=1 ∂ W ∂ Ī4 (i) ∂ Ī4 (i) ∂C } + J ∂U ∂J C -1 (5)
Where S and S vol are respectively the isochoric and volumetric parts of the second Piola-Kirchhoff stress tensor S. The lagrangian deviatoric tangent modulus, is finally given by

C L = 4 { ∂ ∂C ( W, 1 ∂ Ī1 ∂C ) + ∂ ∂C ( W, 2 ∂ Ī2 ∂C ) + n ∑ i=1 ∂ ∂C ( W, (i) 4 ∂ Ī4 (i) ∂C )} (6a) (6b) Where W, k notes ∂ W ∂ Īk
. The eulerian deviatoric modulus is deduced by means of the push-forward operation Φ * as C e = 1 J Φ * ( C L ). Finally the stress and the tangent modulus are given by the first and second derivates of the strain energy function with respect to the first, second and fourth invariants. The general form of C L and C L vol are presented in [START_REF] Kaliske | A formulation of elasticity and viscoelasticity for fibre reinforced material at small and finite strains[END_REF]. For the present model, the particular forms of the derivatives of the energy are expressed as: W,

(i) 4 = 2k w (i) F (i) ( Ī(i) 4 -1) (7) W, (i) 44 = 2k w (i) { ∂F (i) ∂ Ī(i) 4 ( Ī(i) 4 -1) + F (i) } (8) W, (i) 14 = 2k w (i) ∂F (i) ∂ Ī1 ( Ī(i) 4 -1) (9)
with the following conditions

if Ī1max = I 1 ∂F (i) ∂ Ī1 = 0 elseif ∂F (i) ∂ Ī1 = η 1 2 ( Ī1max -Ī1 Ī1max -3 ) -1 2 ( 1 Ī1max -3 ) ( Ī(i) 4 max - Ī(i) 4 Ī(i) 4 max -1 ) ( Ī(i) 4 max Ī4 max ) 4 (10) if Ī(i) 4 max = Ī(i) 4 ∂F (i) ∂ Ī(i) 4 = 0 elseif ∂F (i) ∂ Ī(i) 4 = η √ Ī1max -Ī1 Ī1max -3 ( 1 Ī(i) 4 max -1 ) ( Ī(i) 4 max Ī4 max ) 4 (11)
Finally, the model is implemented via a UMAT in the Finite Element code Abaqus and simulations can be realized.

Validation of the model on a structure calculus

Presentation of an experimental holey plate tension test

In order to test the ability of the model to simulate adequately a complex structure, an experimental test was designed. A RTV-3428 silicone rectangular plate sample with 2mm in thickness was molded using the same protocol presented in [START_REF] Machado | Analysis of the isotropic models of the Mullins effect based on filled silicone rubber experimental results[END_REF]. After molding, five holes were generated. The final dimensions of the holey plate are described in Fig. 4(a). The plate was put into a tensile device which allows to record the force undergone by the plate. The force measure is synchronised with a Stereo Digital Image Correlation (SDIC) system. The SDIC technique allows the determination of three dimensional surface displacements and strain fields of the plate (Machado et al., 2012a). Thus, under the incompressibility assumption, in each point of the plate, strain tensor is evaluated. The aim of this paper is to observe the Mullins effect, so a cyclic tensile test is designed. The strain field of the plate is studied for a macroscopic engineering strain of 50% after different maximal engineering strain: 50 %, 100% and 150%. It is chosen to study particularly three paths in the

Comparaison with model predictions

Finite element simulations of the experiments were performed with the model in plane stress hypothesis with linear elements. The objective of this section is to compare the experimental and numerical results for an engineering strain of 50%. The following notations are used in this part, 1C EF 50 and 1C exp 50 represent the first load for the numerical model and experimental results, respectively. In the same way, 2C EF 100 , 2C exp 100 , and 2C EF 150 , 2C exp 150 represent the deformation at 50% after a first load of 100% and 150%, respectively.

Global behavior of the plate

The global displacement and the force undergone by the plate for the experimental and numerical tests are compared in Fig 5 . It appears that the first and second loadings are very close. A small difference appears for the first loading but, it should be due to the underestimation of biaxial loading by the model. This small difference is only due to the hyperelastic constitutive equation and not to the stress softening function.

Contours of the plate and strain fields measurements

Fig 6 presents the maximal principal logarithmic strain field obtained for the plate by finite element simulation and experimental measurement for 2C 100 configurations. Considering the strain field distribution along the plate surface, it appears that the FE simulation represents quite well the experimental data. This is also valid for 1C and 2C 150 configurations but these results are not presented in the paper. The contours of the plate are compared for 1C 50 , 2C 100 and 2C 150 configuration and the results are presented in Fig 7 for a deformation of 50%. The experimental tests are the gray picture and the FE simulations the contour lines in red. It is possible to observe that experimental results and numerical simulations are almost superposed.

Observation of specific paths

To focus on the local behavior of the plate, the maximal principal logarithmic strain along the different paths defined in Fig. 4 is presented in Fig. 8. A comparaison between the experimental and numerical tests is presented. According to the Fig. 8(a) the maximal principal logarithmic strain is greater for the second loadings than for the first loadings, but Fig. 8(b) and (c) show that the maximal principal logarithmic strain is more important for the first loadings that for the second loadings. This difference is explained by the local state of deformation reached along each path. The maximal deformation is reached along the path 1 so the stress-softening too. Since the stress-softening is maximal along the path1 it generates more important deformation for the along this same path seconds loadings . Along the paths 2 and 3 the stress-softening is smaller that along 8 the path1 so the deformations generated along these paths decrease for the second loadings since the global deformation of the holey is still the same as the first loading. Fig. 8 also shows that the results are quite similar even if there is not a perfect spatial superposition between experimental and FE simulations. This difference can be also explained by the out-ofplane deformation, induced by the geometry of the plate, during the experimental test. Indeed for the experimental tests an out-of-plane displacement, close to the holes of the silicone plate, was observed. This out-of-plane displacement does not occur in numerical simulation since it is defined in plane stress. Thus, the superposition of the results is better along the path 2 (Fig. 8(b)) than along the paths 1 and 3 (cf. Fig. 8(a) and (c)) since the out-of-plane displacement observed is minimal near the hole of the path2 and maximal near the holes of paths 1 and 3.

Conclusion

A new model was presented to take into account the anisotropy induced by stress softening during deformation of a rubber like material. The proposed form depends both of global maximal equivalent strain deformation and directinal maximal strains. The model was written using only strain invariants which permits to have a formulation well adapted to finite element implementation. Details about tangent modulus calculation were given in the paper.

The ability of the proposed model to simulate experimental data was tested first on classical tests (uniaxial, planar and equibiaxial tension) and second on tests where the loading directions were changing between first and second loadings. The analysis of a holey plate was realized by means of digital image correlation which allows to validate the global force-displacement response and local deformation of the plate following three different paths (with diffrent type of local strain). Finally, all the results highligt the good efficiency of the anisotropic stress-softening function which depends on only one parameter. This model can be used in finite element softwares. It allows to take account very particular loading paths.

Figure 1 :

 1 Figure 1: Comparison of the model (solid lines) with experimental data (dotted lines) from Machado et al. (2010) (a) cyclic uniaxial tensile test, (b) cyclic pure shear test (c) cyclic equibiaxial test.

Figure 2 :

 2 Figure 2: Comparison of the model with uniaxial pre-stretching experimental data (dotted lines) from Machado et al. (2012b). (a) Second tensile test oriented at 90 • (b) Second tensile test oriented at 45 • (e) Second tensile test oriented at 30 • (f) Second tensile test oriented at 0 •

Figure 3 :

 3 Figure3: Comparison of the model (solid lines) with biaxial pre-stretching experimental data (dotted lines). Curve a presents simulation of the model for the second load after an equibiaxial test; curve g presents simulation of the model for the second load after a biaxial test of biaxiality rate µ = 0.7; curve h presents simulation of the model for the second load after a biaxial test of biaxiality rate µ = 0.5

Figure 4 :

 4 Figure 4: (a) Silicone plate sample containing five holes and paths of interest, (b) Biaxiality ratio along the three paths

Figure 5 :

 5 Figure 5: Comparison of global force response between the proposed model and the experimental holey plate.

Figure 6 :Figure 7 :

 67 Figure 6: Maximal strain logaritmic fields of the silicone plate for the first loading : (a) numerical results and (b) experimental results.

Figure 8 :

 8 Figure 8: Maximal principal deformation for a global strain of 50% along the path1 (a), the path2(b) and the path3 (c).


Nathann Cohen
email: nathann.cohen@gmail.com

Several Graph Problems and their LP formulation

is meant as an explanation of several graph theoretical functions defined in Sage's Graph Library (http://www.sagemath.org/), which use Linear Programming to solve optimization of existence problems.

Maximum average degree

The average degree of a graph G is defined as ad(G) = 2|E(G)| |V (G)| . The maximum average degree of G is meant to represent its densest part, and is formally defined as :

mad(G) = max H⊆G ad(H)
Even though such a formulation does not show it, this quantity can be computed in polynomial time through Linear Programming. Indeed, we can think of this as a simple flow problem defined on a bipartite graph. Let D be a directed graph whose vertex set we first define as the disjoint union of E(G) and V (G). We add in D an edge between (e, v) ∈ E(G) × V (G) if and only if v is one of e's endpoints. Each edge will then have a flow of 2 (through the addition in D of a source and the necessary edges) to distribute among its two endpoints. We then write in our linear program the constraint that each vertex can absorb a flow of at most z (add to D the necessary sink and the edges with capacity z).

Clearly, if H ⊆ G is the densest subgraph in G, its |E(H)| edges will send a flow of 2|E(H)| to their |V (H)| vertices, such a flow being feasible only if z ≥ 2|E(H)|
|V (H)| . An elementary application of the max-flow/min-cut theorem, or of Hall's bipartite matching theorem shows that such a value for z is also sufficient. This LP can thus let us compute the Maximum Average Degree of the graph.

LP Formulation :

• Mimimize : z • Such that : -a vertex can absorb at most z ∀v ∈ V (G), e∈E(G) e∼v
x e,v ≤ z each edge sends a flow of 2 ∀e = uv ∈ E(G), x e,u + x e,v = 2

• x e,v real positive variable REMARK : In many if not all the other LP formulations, this Linear Program is used as a constraint. In those problems, we are always at some point looking for a subgraph H of G such that H does not contain any cycle. The edges of G are in this case variables, whose value can be equal to 0 or 1 depending on whether they belong to such a graph H. Based on the observation that the Maximum Average Degree of a tree on n vertices is exactly its average degree (= 2 -2/n < 1), and that any cycles in a graph ensures its average degree is larger than 2, we can then set the constraint that z ≤ 2 -2 |V (G)| . This is a handy way to write in LP the constraint that "the set of edges belonging to H is acyclic". For this to work, though, we need to ensure that the variables corresponding to our edges are binary variables. corresponding patch : http://trac.sagemath.org/sage_trac/ticket/7529 Given a graph G whose edges are weighted by a function w : E(G) → R, a solution to the T SP is a hamiltonian (spanning) cycle whose weight (the sum of the weight of its edges) is minimal. It is easy to define both the objective and the constraint that each vertex must have exactly two neighbors, but this could produce solutions such that the set of edges define the disjoint union of several cycles. One way to formulate this linear program is hence to add the constraint that, given an arbitrary vertex v, the set S of edges in the solution must contain no cycle in G -v, which amounts to checking that the set of edges in S no adjacent to v is of maximal average degree strictly less than 2, using the remark from section 1.

We will then, in this case, define variables representing the edges included in the solution, along with variables representing the weight that each of these edges will send to their endpoints.

LP Formulation :

• Mimimize e∈E(G) w(e)b e
• Such that : This problem is polynomial by a result from Edmonds. Obviously, nothing ensures the following formulation is a polynomial algorithm as it contains many integer variables, but it is still a short practical way to solve it. This problem amounts to finding, given a graph G and an integer k, edge-disjoint spanning trees T 1 , . . . , T k which are subgraphs of G. In this case, we will chose to define a spanning tree as an acyclic set of |V (G)|-1 edges.

-Each vertex is of degree 2 ∀v ∈ V (G), e∈E(G) e∼v b e = 2 -No cycle disjoint from a special vertex v * * Each edge sends a flow of 2 if it is taken ∀e = uv ∈ E(G -v * ), x e,u + x e,v = 2b e * Vertices receive strictly less than 2 ∀v ∈ V (G -v * ), e∈E(G) e∼v x e,v ≤ 2 - 2 |V (G)| • Variables -x e,v

LP Formulation :

• Maximize : nothing

• Such that :

-An edge belongs to at most one set

∀e ∈ E(G), i∈[1,...,k] b e,k ≤ 1 -Each set contains |V (G)| -1 edges ∀i ∈ [1, . . . , k], e∈E(G) b e,k = |V (G)| -1 -No cycles * In each set, each edge sends a flow of 2 if it is taken ∀i ∈ [1, . . . , k], ∀e = uv ∈ E(G),
∀v ∈ V (G), e∈E(G) e∼v x e,v ≤ 2 - 2 |V (G)| • Variables :
b e binary (is e in the tree ?) c v binary (does the tree contain v ?) x e,v real positive variable (flow sent by the edge) corresponding patch : http://trac.sagemath.org/sage_trac/ticket/8403

The linear arboricity of a graph G is the least number k such that the edges of G can be partitionned into k classes, each of them being a forest of paths (the disjoints union of paths -trees of maximal degree 2). The corresponding LP is very similar to the one giving edge-disjoint spanning trees LP Formulation :

• Maximize : nothing

• Such that :

-An edge belongs to exactly one set

∀e ∈ E(G), i∈[1,...,k] b e,k = 1 -Each class has maximal degree 2 ∀v ∈ V (G), ∀i ∈ [1, . . . , k], e∈E(G) e∼v b e,k ≤ 2 -No cycles * In each set, each edge sends a flow of 2 if it is taken ∀i ∈ [1, . . . , k], ∀e = uv ∈ E(G), x e,k,u + x e,k,v = 2b e,k * Vertices receive strictly less than 2 ∀i ∈ [1, . . . , k], ∀v ∈ V (G), e∈E(G) e∼v x e,k,v ≤ 2 - 2 |V (G)| • Variables
b e,k binary (is edge e in set k ?) x e,k,u positive real (flow sent by edge e to vertex u in set k)

Acyclic edge coloring

An edge coloring with k colors is said to be acyclic if it is proper (each color class is a matching -maximal degree 1), and if the union of the edges of any two color classes is acyclic. The corresponding LP is almost a copy of the previous one, except that we need to ensure that k 2 different classes are acyclic.

corresponding patch : http://trac.sagemath.org/sage_trac/ticket/8405

For more information on minor theory, please see http://en.wikipedia.org/wiki/Minor_%28graph_theory%29 It is a wonderful subject, and I do not want to begin talking about it when I know I couldn't freely fill pages with remarks :-) For our purposes, we will just say that finding a minor H in a graph G, consists in :

1. Associating to each vertex h ∈ H a set S h of representants in H, different vertices h having disjoints representative sets 2. Ensuring that each of these sets is connected (can be contracted)

3. If there is an edge between h 1 and h 2 in H, there must be an edge between the corresponding representative sets

Here is how we will address these constraints :

1. Easy 2. For any h, we can find a spanning tree in S h (an acyclic set of |S h | -1 edges)

3. This one is very costly.

To each directed edge g 1 g 2 (I consider g 1 g 2 and g 2 g 1 as different) and every edge h 1 h 2 is associated a binary variable which can be equal to one only if g 1 represents h 1 and g 2 represents g 2 . We then sum all these variables to be sure there is at least one edge from one set to the other.

LP Formulation :

• Maximize : nothing

• Such that :

-A vertex g ∈ V (G) can represent at most one vertex h ∈ V (H) ∀g ∈ V (G), h∈V (
∈ V (G), e∈E(G) e∼v x e,k,v ≤ 2 - 2 |V (G)|
arc (g1,g2),(h1,h2) can only be equal to 1 if g 1 g 2 is leaving the representative set of h 1 to enter the one of h 2 . (note that this constraints has to be written both for g 1 , g 2 , and then for g 2 , g 1 ) ∀g 1 , g 2 ∈ V (G), g 1 = g 2 , ∀h 1 h 2 ∈ E(H) arc (g1,g2),(h1,h2) ≤ rs h1,g1 and arc (g1,g2),(h1,h2) ≤ rs h2,g2

-We have the necessary edges between the representative sets ∀h 1 h 2 ∈ E(H) ∀g1,g2∈V (G),g1 =g2 arc (g1,g2),(h1,h2) ≥ 1

• Variables

rs h,g binary (does g represent h ? rs = "representative set") t e,h binary (does e belong to the spanning tree of the set representing h ?) x e,v real positive (flow sent from edge e to vertex v) arc (g1,g2),(h1,h2) binary (is edge g 1 g 2 leaving the representative set of h 1 to enter the one of h 2 ?) corresponding patch : http://trac.sagemath.org/sage_trac/ticket/8404

x

  real positive variable (flow sent by the edge) b e binary (is the edge in the solution ?) corresponding patch : http://trac.sagemath.org/sage_trac/ticket/7529

--

  x e,k,u + x e,k,u = 2b e,k * Vertices receive strictly less than 2 ∀i ∈ [1, . . . , k], ∀v ∈ V (G), e∈E(G) e∼v x e,k,v ≤ 2 -2 |V (G)| • Variables b e,k binary (is edge e in set k ?) x e,k,u positive real (flow sent by edge e to vertex u in set k) corresponding patch : http://trac.sagemath.org/sage_trac/ticket/7476 Finding a spanning tree in a Graph G can be done in linear time, whereas computing a Steiner Tree is NPhard. The goal is in this case, given a graph, a weight function w : E(G) → R and a set S of vertices, to find the tree of minimum cost connecting them all together. Equivalently, we will be looking for an acyclic subgraph Hof G containing |V (H)| vertices and |E(H)| = |V (H)| -1 edges, which contains each vertex from S Each vertex from S is in the tree ∀v ∈ S, e∈E(G) e∼v b e ≥ 1 c is equal to 1 when a vertex v is in the tree ∀v ∈ V (G), ∀e ∈ E(G), e ∼ v, b e ≤ c v The tree contains |V (H)| vertices and |E(H)| = |V (H)| -* Each edge sends a flow of 2 if it is taken ∀e = uv ∈ E(G), x e,u + x e,u = 2b e,k * Vertices receive strictly less than 2

  An edge e can only belong to the tree of h if both its endpoints represent h ∀e = g 1 g 2 ∈ E(G), t e,h ≤ rs h,g1 and t e,h ≤ rs h,g2 -In each representative set, the number of vertices is one more than the number of edges in the corresponding tree ∀h,

			rs h,g ≤ 1
			H)
	-g∈V (G)	rs h,g -	e∈E(G)	t e,h = 1
	-No cycles in the union of the spanning trees

* Each edge sends a flow of 2 if it is taken ∀e = uv ∈ E(G), x e,u + x e,v = 2 h∈V (H) t e,h * Vertices receive strictly less than 2 ∀v


Houda Zaidi

Laurent Santandrea

Guillaume Krebs

Yann Le Bihan

Edouard Demaldent

FEM Technique for Modeling Eddy-Current Testing of Ferromagnetic Media With Small Skin Depth

The modeling of eddy-current testing problems, by means of finit elements, requires us to carefully consider the skin effect. When the material to inspect has a small skin depth compared with the other dimensions of the problem, which is a typical situation in the case of ferromagnetic media, the classical FEMs can become unsuitable in terms of mesh quality or computation time. In this paper, the overlapping FEM is implemented to effi iently consider the skin effect. 3-D problems are studied and the validity of the proposed approach is shown. Index Terms-Eddy currents, f nite element analysis, magnetic materials, magnetic multilayers, non-destructive testing.

I. INTRODUCTION

T HE non-destructive testing by means of the ECT is nowadays frequently used to evaluate the properties of conducting materials. The interpretation of the signal variations received by the ECT probe (due for instance to a crack) involves the use of mathematical models. These models permit also to study conf gurations that are diff cult to reproduce experimentally.

The FEM is often chosen due to its high adaptability when using simplex elements [START_REF] Hamia | Eddy-current nondestructive testing using an improved GMR magnetometer and a single wire as inducer: A FEM performance analysis[END_REF] (triangles in 2-D or tetrahedra in 3-D). However, a correct modeling of ECT requires some rules: 1) accurate meshes for the probe and the surrounding air; 2) fi e mesh in the liftoff (distance between the probe and the conducting part); and 3) usually (considering f rstorder elements) at least three layers of f nite elements have to be used in the skin depth of the inspected conductor. This last condition becomes critical when considering ECT of ferromagnetic materials having a low skin depth. Indeed the meshing of these materials will result either in highly distorted simplex fi ite elements or in a high density mesh leading to a huge number of degrees of freedom (DOF). This paper proposes an eff cient technique for modeling ECT of ferromagnetic materials based on the overlapping method. In the previous works, numerical techniques based on this method have been proposed to consider a thin geometrical domain thanks to a layer of overlapping elements [START_REF] Zaidi | Eff cient FEM modelling of printed coils for eddy current testing[END_REF]- [START_REF] Krebs | Overlapping finit elements used to connect non-conforming meshes in 3D with a vector potential formulation[END_REF]. In this paper, we propose to adapt this technique to model eff ciently the skin effect in low skin depth domains using several layers of overlapping elements. The overlapping method is used to obtain a prismatic mesh in the f rst skin depths of the ferromagnetic material and to connect this mesh to the initial tetrahedral mesh of the remaining of the study domain.

The f rst part of this paper describes brief y the used f nite element formulations. In the second part, the overlapping FEM is presented. Then, some explanations are given about the proposed approach to consider the skin effect. Finally, two test cases are considered and results obtained with the implemented method are compared with analytical and experimental data. The f rst test case includes a ferromagnetic piece free from defect, while the second one is dedicated to model a fl wed ferromagnetic piece.

II. FINITE ELEMENT FORMULATIONS

There are two families of formulations for solving Maxwell's equations. One is based on the calculation of the electric f eld e and the other on the calculation of the magnetic fiel h [START_REF] Ren | Computation of 3-D electromagnetic f eld using differential forms based elements and dual formulations[END_REF]. In either electric or magnetic formulation, potentials are generally used to describe the f eld. In the magnetic formulation, the magnetic fiel h is described as a function of the electric vector potentials t and t 0 and with the help of the magnetic scalar potential ϕ h = t + t 0 -grad ϕ.

(1)

The electric current vector potentials t and t 0 come from the relations curl t = j and curl t 0 = j 0 ,w he rej is the density of eddy currents and j 0 is the density of source current. In the electric formulation, with the introduction of a magnetic vector potential a and an electric scalar potential ψ, the electrical fiel e is expressed as

e =-∂ t (a + grad ψ). (2)
The potential a comes from the relation b =curl a,w h e r e b is the magnetic fl x density. The potentials are discretized by Whitney elements. In this paper, the 3-D geometry is discretized using tetrahedral. The modeling is made in harmonic regime. The obtained algebraic system is solved using a Jacobi preconditioned conjugate gradient method.

III. OVERLAPPING FEM

The overlapping FEM has been f rst introduced in 2-D and applied recently in 3-D for ECT applications in [START_REF] Zaidi | Modeling of thin conducting and magnetic layers in eddy current testing by overlapping f nite elements[END_REF] and [START_REF] Tsukerman | Overlapping finit element for problems with movement[END_REF]. Its principle is to connect two domains separated by an initially unmeshed region [see Fig. 1(a)]. The shape functions belonging to the nodes i and j are def ned only in D 1 and D 2 , This double projection creates overlapping elements. These elements are the support of shape functions, which permit to ensure the continuity of the interpolated values. The projected nodes [stars in Fig. 1(b)] are virtual and are not the support of DOF. They are just used to determine integration areas [hatched area in Fig. 1(b)], where the different integral terms of the FEM formulation are evaluated. New shape functions are introduced in D 0 . These functions are represented in Fig. 1(b) by quadrilateral elements (or prisms in 3-D) that overlap in D 0 . In this region, integral terms are computed in the area where shape functions are commonly non-zero [see the hatched area in Fig. 1(b) related to the nodes i and j ] and added to the matrix system. Depending on the used formulation, new vertical edges have to be added in the problem.

For each quadrilateral element (or prismatic element in 3-D), nodal and edge unknowns are considered. The nodal unknowns are the four nodes of the two edges from which derives the quadrilateral element (or the six nodes of the two triangles from which derives the prism). For the hatched area in Fig. 1(c), the nodal unknowns are associated to N 1 , N 2 , N 3 ,a n dN 4 . In a vector potential formulation, unknowns on vertical edges are added [START_REF] Krebs | Overlapping finit elements used to connect non-conforming meshes in 3D with a vector potential formulation[END_REF]. For the hatched area in Fig. 1(c), the edge unknowns are associated to the edges e 1 , e 2 , e 3 ,ande 4. It follows that the edge unknowns are composed of the two edges of the two lines that derives the quadrilateral element (or the six edges of the two triangles) in addition the two vertical edges of the quadrilateral (or the three vertical edges of the prism). This is shown, in the 3-D case, in Fig. 2.

IV. MODELING OF THE SKIN EFFECT

In the ECT of ferromagnetic materials, because of the high value of the magnetic permeability, the skin depth is usually small compared with other geometrical parameters. The eddycurrent density is important in the skin depth and decreases quickly and approximately exponentially. A f ne mesh is then required in the skin depth of the ferromagnetic material. However, since the skin depth is very small, it is diff cult to generate a specif c mesh that contains enough elements in skin depth without a very high total number of DOF in the mesh and, above all, without reducing the quality of the mesh. Indeed, the meshing of the skin depth with several layers of equilateral elements leads to a large number of unknowns. On the contrary, a less dense mesh can lead to highly distorted elements and sometimes to an insuff cient discretization to model induced currents.

The proposed approach considers the f rst skin depths with multiple layers of overlapping elements. With this approach, it is possible to consider the high skin effect occurring in ferromagnetic materials without the previously presented mesh diff culties. Indeed, the principle of this method is to mesh the area with prismatic elements which, unlike the tetrahedral elements, can be flattene to a large extent without critically degrading their quality. It is possible to use overlapping element layers with a non-uniform thickness distribution to consider the f eld evolution in the depth. For example, an exponential distribution according to the depth, while respecting the requirement to have at least three or four layers of elements in the skin depth, should allow to discretize f nely the material surface and more coarsely its depth with the purpose to limit the number of DOF. Fig. 3 shows the overlapping layer distribution adopted in this paper.

With the overlapping method, different triangular discretization can also be adopted from one layer to the other, notably in the case of materials of non-canonical shape. This justif es the choice of the overlapping method compared with other methods such as shell element method.

V. T EST CASE AND COMPARISONS

A. Studied Configu ation

In this paper, two test cases to validate the implementation of the overlapping approach for the ECT of ferromagnetic materials will be treated. The fir t problem concerns the ECT of a ferromagnetic medium free from defect. The aim of the second problem is to obtain the response of an ECT probe to a f aw (slot) in a ferromagnetic material. Fig. 4(a) [respectively Fig. 4(b)] shows the conf guration of the f rst test case (respectively second). The parameters of the two problems are given in Table I [7].

B. Flaw-Free Material

The studied problem is a simple coil above a ferromagnetic plate without defect. It comes from a conf guration proposed in [START_REF] Ditchburn | Eddy-current nondestructive inspection with thin spiral coils: Long cracks in steel[END_REF] with associated experimental and analytical results. The main parameters of the problem are given in Table I (except for the f aw parameters, which are irrelevant here).

The working frequency range is very large from 0.1 kHz to 100 MHz. Consequently, the skin depth in the plate varies considerably from 2.9 mm to 2.9 µm. It was chosen to work with 10 overlapping layers whose thicknesses are equal to (1 ≤ i ≤ 10)

t i = (z i -z i-1 ) with z i = ln(11/(11 -i ))δ (3)
where δ is the skin depth of the plate. By this way, the fi st 2.39 δ of the plate depth are discretized with the overlapping method (Fig. 3). The remainder of the plate (inside) is meshed with tetrahedral elements. For the lowest frequencies (≤10 kHz), the thicknesses of the overlapping layers have been kept constant to avoid elements of too large thickness. For these low frequencies, the classical FEM could be used.

In the considered problem, the thickness of the liftoff (coil plate distance) is not very small compared with the other dimensions of the problem; hence, the liftoff domain is meshed with tetrahedral elements. The proposed approach has been implemented and the impedance variation Z = R + j X of the coil due to the presence of the plate is calculated. The evolution of R and X normalized by the reactance X 0 of the coil in air is shown in Fig. 5 over the frequency range. These results are compared with analytical and experimental data issued from [START_REF] Ditchburn | Eddy-current nondestructive inspection with thin spiral coils: Long cracks in steel[END_REF].

The two formulations presented previously have properties of duality [START_REF] Ren | Computation of 3-D electromagnetic f eld using differential forms based elements and dual formulations[END_REF]. To obtain better results for both magnetic and electric f elds, it is recommended to solve a given problem with both formulations and use averaged values.

It can be noticed that the results are in good agreement. The method allows extending very signif cantly the frequency range of application of the FEM by comparison with classical FEM. Limitations of the classical FEM with a meshing constituted only of tetrahedral elements (and of other numerical techniques) can be found in [START_REF] Zorni | Modelling eddy current testing of ferromagnetic medium[END_REF], where the same test case was considered. This reference shows that the classical FEM results begin to diverge for high frequency (above several kilohertz) due to the increasing skin effect. This demonstrates the interest of the approach developed in this paper.

C. Flawed Plate

In this section, the overlapping method is applied to simulate the ECT signal induced by a f aw in a ferromagnetic plate. The probe is centered on the fl w that is constituted by a long slot whose parameters are presented in Table I. The other parameters of this second problem are the same as the previous one.

The coil impedance variation resulting from the presence of the slot is calculated for the same wide range of frequency than previously (from 0.1 kHz to 100 MHz).

For this example, the top part of the f aw is also considered with overlapping elements. To do this, either S 1 or S 2 (boundary surfaces of the overlapping domain, Fig. 2) has to contain a 2-D mask of the slot [Fig. 6(a)]. The discretization of the two surface S 1 and S 2 can be different. As previously, the f rst 2.39 δ of the plate depth and of the slot is meshed with the overlapping method. The remaining depth of the plate and of the slot is meshed with tetrahedral elements.

The average value (calculated from the two formulations) of impedance variation Z = R + j X of the coil due to the defect is calculated. Fig. 7 shows the evolution of R and X normalized by the reactance of the coil in air, in function of the frequency. These results are compared with analytical and experimental data issued from [START_REF] Ditchburn | Eddy-current nondestructive inspection with thin spiral coils: Long cracks in steel[END_REF].

These results show a good agreement between modeling results and measurements of [START_REF] Ditchburn | Eddy-current nondestructive inspection with thin spiral coils: Long cracks in steel[END_REF] and the results provided by the overlapping method. However, it can be noticed that the overlapping results begin to deviate from the analytical ones above several megahertz. From the overlapping method point of view, this discrepancy could be due to the triangular discretization of the overlapping domain in its horizontal plane, which limits the smallness of the f nite elements around the f aw. Nevertheless, the applicability of the FEM for the modeling of f aw detection in ferromagnetic media is considerably increased using the overlapping method.

VI. CONCLUSION

In this paper, a new approach based on the overlapping method for modeling the ECT of thin skin depth media, such as ferromagnetic media, has been proposed. This approach has been applied into two test cases involving f aw free and fl wed ferromagnetic plate, can concluded from the tests that the overlapping method allows an eff cient modeling of the ECT of low skin depth media. This method can be also applied to consider non-homogeneous media such as stratif ed media. It allows meshing thin regions with adapted elements with a great fl xibility, thanks to the non-conformity of the mesh provided by the method. It can also be noticed that with this method the displacement of the probe can be considered without remeshing.

Fig. 1 .

 1 Fig. 1. (a) Two domains independently meshed. (b) Nodal projections. (c) Nodal and edge unknowns associated with a quadrilateral element.

Fig. 2 .

 2 Fig. 2. (a) Two meshes separated by a free volume D 0 . (b) Two surfaces that bound D 0 .( c )M e s ho fD 0 . (d) Unknowns associated with an integration area (prism).

Fig. 3 .

 3 Fig. 3. Discretization of the fir t skin depths of the ferromagnetic domain in several layers of overlapping elements.

Fig. 4 .

 4 Fig. 4. Configur tions of the ECT test cases. (a) Defect-free ferromagnetic medium. (b) Flawed ferromagnetic medium.

Fig. 5 .

 5 Fig. 5. (a) Normalized reactance and (b) resistance variation resulting from the presence of the plate as a function of the frequency.

Fig. 6 .

 6 Fig. 6. (a) Traces of the tetrahedral meshes on the top and on the bottom of the overlapping domain. (b) Prismatic mesh of plate and the slot.

Fig. 7 .

 7 Fig. 7. (a) Normalized reactance and (b) resistance variation resulting from the presence of the slot in function of the frequency.

TABLE I PARAMETERS

 I OF THE TEST CASE[START_REF] Ditchburn | Eddy-current nondestructive inspection with thin spiral coils: Long cracks in steel[END_REF]

ACKNOWLEDGMENT

This work was supported by DIGITEO (Île-de-France Region) and CEA LIST.


large ohmic drop proportional to the inter electrode distance [START_REF] Rozendal | Towards practical implementation of bioelectrochemical wastewater treatment[END_REF]. (ii) its low soluble COD load (150e400 mg/L [START_REF] Verbanck | Major ion content of urban wastewater : assessment of per capita loading[END_REF]), which limits the kinetics of the bioanode when no artificial carbon source, such as acetate, is added to the electrolyte [START_REF] Peng | Nernst-pingpong model for evaluating the effects of the substrate concentration and anode potential on the kinetic characteristics of bioanode[END_REF][START_REF] Blanchet | Comparison of synthetic medium and wastewater used as dilution medium to design scalable microbial anodes: application to food waste treatment[END_REF]. (iii) its pH, often close to neutrality (6.5e7.5 [START_REF] Dong | A 90-liter stackable baffled microbial fuel cell for brewery wastewater treatment based on energy selfsufficient mode[END_REF][START_REF] Escapa | Scaling-up of membraneless microbial electrolysis cells (MECs) for domestic wastewater treatment: bottlenecks and limitations[END_REF]), which is unfavourable for the electrochemical reduction of water into hydrogen. Commercial electrolysers operate more conven tionally with concentrated alkaline electrolytes (alkaline electrolysis [START_REF] Zeng | Recent progress in alkaline water electrolysis for hydrogen production and applications[END_REF]) or strongly acidic pH (polymer electrolyte membrane fuel cell technologies [START_REF] Nuttall | General Electric's Solid Polymer Electrolyte Water Electrolysis, Direct Energy Convers[END_REF][START_REF] Carmo | A comprehensive review on PEM water electrolysis[END_REF]). (iv) the microbial production of gaseous by products (H 2 S, CO 2 , ammonia, CH 4 ), which are related to the microbial activity of the aerobic, anaerobic and fermenting consortia and degrade the purity of the hydrogen gas produced. Using a membrane to separate anodic and cathodic compartments limits gas mixing but also creates a barrier to the exchange of ions between compartments, increasing the ohmic drop and creating a pH gradient in the cell [START_REF] Zeppilli | Anion vs cation exchange membrane strongly affect mechanisms and yield of CO2 fixation in a microbial electrolysis cell[END_REF].

Many areas of improvement have been proposed to reduce the U cell of MECs fed with domestic wastewater. The anode over potential can be reduced by using 3D electrodes or porous elec trodes that promote biofilm development thanks to their huge specific surface area [START_REF] Xie | Design and fabrication of bioelectrodes for microbial bioelectrochemical systems[END_REF]. The cathode overpotential can be lowered by using metallic catalysts, e.g. Ni based catalysts [START_REF] Mitov | Nickel-based electrodeposits as potential cathode catalysts for hydrogen production by microbial electrolysis[END_REF][START_REF] Selembo | Hydrogen production with nickel powder cathode catalysts in microbial electrolysis cells[END_REF], applied to the cathode surface or biological catalysts naturally forming from wastewater [START_REF] Bergel | Catalysis of oxygen reduction in PEM fuel cell by seawater biofilm[END_REF][START_REF] He | Application of bacterial biocathodes in microbial fuel cells[END_REF]. The ohmic drop in the MEC can be decreased by designing specific reactor geometries that minimize the distance between electrodes [START_REF] Hou | Microbial electrolysis cell with spiral wound electrode for wastewater treatment and methane production[END_REF] or by using membrane electrode assem blies [START_REF] Zhang | Treating refinery wastewaters in microbial fuel cells using separator electrode assembly or spaced electrode configurations[END_REF].

Homogeneous catalysis has, so far, been little exploited for reducing the cathode overpotential of MECs. Weak acids have been demonstrated to catalyse hydrogen evolution reaction on metallic materials [START_REF] Silva | Electrochemical deprotonation of phosphate on stainless steel[END_REF][START_REF] Silva Munoz | Combining phosphate species and stainless steel cathode to enhance hydrogen evolution in microbial electrolysis cell (MEC)[END_REF]. The mechanism of water reduction catalysed by weak acids has been investigated in depth with phosphate species [START_REF] Silva Munoz | Hydrogen production by electrolysis of a phosphate solution on a stainless steel cathode[END_REF] and has been described as a four step process (HA and A À are the non dissociated and dissociated forms, respectively):

HA þ e À 4 H ad þ A À (2)
H ad þ HA þ e À 4 H 2 þ A À (3)
2H ad 4 H 2 (4)

A À þ H 2 O 4 HA þ OH À (5)
The use of carbonate instead of phosphate as a homogeneous catalyst would have many environmental and economic advan tages. Carbonate has less effect on the environment than phos phate, which causes problems of eutrophication and dystrophication in aquatic environments, and carbonate is nearly 50% cheaper than phosphate.

The purpose of this work was to assess the capacity of carbonate to catalyse H 2 evolution on metallic cathodes. The effect of the carbonate concentration on the kinetics of water reduction was studied. The kinetics of water reduction in carbonate solutions and in domestic wastewater were compared in order to estimate the cathode potential (E c ) gain due to the carbonate solution. Experi ments were conducted at pH ranging from 7.0 to 12.0 so as to un derstand which carbonate species (CO 3 2À , HCO 3 À , H 2 CO 3 ) were involved in the mechanism of water reduction catalysis. Finally, various grades of stainless steels and Ni alloys were compared to graphite as cathode material to determine which one would be the most suitable to enhance the catalytic effect of carbonate.

Materials and methods

Chemicals

The chemical products used in the experiments were: potas sium hydrogen carbonate (KHCO 3 , 99.5%; Sigma Aldrich), potas sium carbonate (K 2 CO 3 , 99þ%; Acros Organics), potassium chloride (KCl, 99þ%; Acros Organics), potassium hydroxide (KOH, 85%; Acros Organics), and hydrochloric acid (HCl, 37%; Acros Organics).

All solutions were prepared with deionized water (Elga Purelab Option R; 15 MU cm). The pH of the solutions was measured with a pH meter (Eutech Instruments) and, in some cases, adjusted using 3 M solutions of KOH or HCl. The conductivity of the solutions was measured with an EC meter (RadioMeter) and, in some cases, adjusted to a given value by adding KCl.

Domestic wastewater

Domestic wastewater was collected at the wastewater treat ment plant of Castanet (France) and stored at 4 C before use.

Materials

Working electrodes (cathodes) were rectangular plates (1.5 cm  2 cm). Current densities were expressed with respect to the surface area of the two sides of the electrode, i.e. 6 cm 2 . Seven materials were tested: stainless steels AISI 316L, 304, 310, 430 and 254SMO; nickel alloy (Ni 80%, Fe 20%), and graphite. The compo sitions of the materials were determined by energy dispersive X ray spectroscopy (EDX). Before each experiment, the electrodes were ground with abrasive discs (P800, P1200, P2400; Presi). No chemical treatment was applied to the electrodes. The current collectors were threaded titanium rods insulated along their length with heat shrinkable tube.

Voltammetric study

Electrochemical experiments were carried out in a 150 mL three electrode cell using a VSP2 potentiostat (Bio Logic SA, France) controlled by the EC Lab software. The reference electrode was a saturated calomel electrode (SCE) and the counter electrode was a DSA plate 1.5 Â 2 cm 2 purchased from Magneto (Netherlands).

The current/potential curves were recorded with the linear sweep voltammetry technique of EC Lab. E c was ramped from the open circuit potential to 1.3 or 1.5 V vs. SCE at 10 mV/s. Each curve was recorded 3 times to ensure the reproducibility of the results. The curves were superimposable so only one was chosen to be displayed in the results.

Unless otherwise stated, carbonate solutions were deaerated by bubbling nitrogen gas at a rate of 10 mL/s for 15 min. The pH was measured after air removal so that its variation due to the N 2 bubbling could be taken into account. The pH of all carbonate so lutions was pH 8.7 after deaeration. The nitrogen flux was main tained above the solution during the experiments.

Results

Preliminary calculations to determine the targeted MEC electrolysis voltage and the average current density

Wastewater treatment by MEC technology can be considered as energy neutral when the energy recovered from the hydrogen produced is equal to the energy supplied to the MEC to produce the hydrogen. The energy supplied (E sup ) depends on the cell voltage (U cell ) according to Equation ( 6):

E sup U cell N eÀ F (6)
where U cell is the cell voltage, N e-is the quantity of electrons extracted by the bioanode and F is the Faraday constant ( 96500 C/ mol).

Assuming that the faradaic yield of hydrogen evolution at the cathode is 100%, the energy recovered (E r ) from the combustion of hydrogen is:

E r N eÀ 2 LCV H2 (7)
where LCV H2 is the net calorific value of hydrogen ( 242.7 kJ/mol). By combining Equations ( 6) and ( 7), the limit value of U cell for E sup E r should not exceed 1.2 V and this voltage does not depend on the number of electrons extracted by the bioanode. Several studies conducted with domestic wastewater fed MECs have dis played an applied voltage lower than 1.2 V but the production rates of hydrogen reported were far below those of conventional in dustrial water electrolysis cells (Table 1).

Basically, bioanodes operate at potentials ranging from 0.4 V vs. SCE to þ0.6 V vs. SCE [START_REF] Rimboud | Electroanalysis of microbial anodes for bioelectrochemical systems: basics, progress and perspectives[END_REF]. In this range of anode potential, the best performances reported in the literature in terms of current density obtained with bioanodes fed with real wastewater were about 10 A/m 2 , which corresponds to a hydrogen production of approximately 4 L/h/m 2 cathode at 100% coulombic efficiency [START_REF] Silva Munoz | Combining phosphate species and stainless steel cathode to enhance hydrogen evolution in microbial electrolysis cell (MEC)[END_REF][START_REF] Pannell | Effect of fed-batch vs. continuous mode of operation on microbial fuel cell performance treating biorefinery wastewater[END_REF][START_REF] Yu | Power densities and microbial communities of brewery wastewater-fed microbial fuel cells according to the initial substrates[END_REF]. This highest value of current density obtained from the oxidation of real domestic wastewater will be used further in this study as a reference of current density for optimizing the cathode reaction.

Influence of the carbonate concentration on the water reduction kinetics

Current/potential reduction curves were recorded with a 254SMO stainless steel cathode immersed in potassium carbonate solutions (Fig. 1). In the first set of experiments, carbonate con centrations ranged from 0.01 to 0.2 M and the solutions contained 0.1 M of KCl as supporting electrolyte (Fig. 1. A&B). In the second set (Fig. 1. C&D), a wider range of carbonate concentrations was tested, up to the saturation concentration (2.8 M). To avoid any possible bias caused by a difference in conductivity, the conductivity of all the solutions was adjusted with KCl to the same constant value of 152 mS/cm. This value corresponds to the conductivity of a satu rated carbonate solution (2.8 M at 20 C). A current/potential curve was also recorded in domestic wastewater without any addition of KCl (conductivity 0.9 S/m) as is commonly done in the literature on wastewater treatment by MEC.

In the concentration range of 0 Me2.8 M of carbonate, the water reduction kinetics improved when the carbonate concentration rose (Fig. 1A and1C).

For carbonate concentrations ranging from 0.01 M to 0.2 M, the current density at 1.3 V vs. SCE was directly proportional to the carbonate concentration (correlation coefficient 0.989) (Fig. 1B). These results are consistent with those obtained by De Silva et al. with phosphate solutions having concentrations ranging between 0.01 M and 0.5 M on a 316L SS cathode at pH 8.0 [START_REF] Silva Munoz | Hydrogen production by electrolysis of a phosphate solution on a stainless steel cathode[END_REF].

In the concentration range 0.2e2.8 M, the relation between carbonate concentration and water reduction rate was no longer linear (Fig. 1D). The electro catalytic mechanism proposed by Da Silva et al. [START_REF] Silva | Electrochemical deprotonation of phosphate on stainless steel[END_REF] involves the deprotonation of the weak acid and the adsorption of the H atom on the cathode surface (reaction (2)). The nonlinear behaviour observed in Fig. 1D may have been due to the saturation of H atom adsorption sites on the cathode surface or the fact that there were no longer any carbonate ion transport related limitations when the carbonate concentration was higher than 0.2 M.

At a current density of 10 A/m 2 , E c in 1 M and 2.8 M carbonate solutions was 1.04 V vs. SCE. This represents a 380 mV gain in E c compared with that obtained at the same current density in do mestic wastewater ( 1.42 V vs. SCE). Working with a 1 M carbonate solution as catholyte instead of domestic wastewater would reduce the value of U cell by 32% compared to the 1.2 V MEC limit voltage calculated in 3.1. This reduction of E c is due to both the increase of the electrolyte conductivity and the catalytic effect of the carbonate species.

Influence of the pH of the carbonate solution on water reduction catalysis

Linear sweep voltammetries were performed on a 254SMO stainless steel electrode in aqueous 1 M carbonate solutions where the pH was adjusted in the range of 7e12 (Fig. 2). The highest water reduction rate was obtained at pH 8. The water reduction rates obtained at pH 7, pH 9 and pH 10 were lower than at pH 8 but still satisfactory. At pH 11 and pH 12, the water reduction rates were significantly lower than at more acidic pH.

According to the predominance diagram of carbonic acid (Fig. 3A), the HCO 3 À form is predominant in the pH range from 6.3 to 10.3 (85% at pH 7; 97% at pH 8; 90% at pH 9; 52% at pH 10). Fig. 3B plots the cathode potential at 10 A/m 2 versus pH (Fig. 3B, tri angles) and the percentage of the HCO 3 À species at the same pH values (Fig. 3B, diamonds). The two curves are almost identical in shape, which implies that HCO 3 À is the key species involved in the water reduction catalysis. a Applied voltages calculated from energy consumption of industrial water electrolysis cells (kW=m 3 H2 ) and maximum H 2 production.

consistent with many studies in the literature in which different associations of Cr, Ni, Mo and Mn are used as alloy elements to catalyse the hydrogen evolution reaction in alkaline media. Nickel alloy seems to catalyse hydrogen production even more efficiently than pure nickel. Bachvarov et al. [START_REF] Bachvarov | Electrodeposited NiFeCo and NiFeCoP alloy cathodes for hydrogen evolution reaction in alkaline medium[END_REF] tested a nickel alloy elec trodeposit, containing iron, cobalt and phosphorus, as a catalyst for hydrogen production through alkaline water electrolysis. The presence of iron, cobalt and phosphorus in nickel alloys signifi cantly decreased the overpotential of the hydrogen evolution re action compared to pure nickel. Gonzales Buch et al. [START_REF] Gonz Alez-Buch | Study of the catalytic activity of 3D macroporous Ni and NiMo cathodes for hydrogen production by alkaline water electrolysis[END_REF] also studied alkaline water electrolysis with 3D macroporous Ni and NiMo cathodes. NiMo electrodes had higher catalytic activity than Ni for the hydrogen evolution reaction. Tang et al. studied a nanosheet of Co Mn carbonate hydroxide deposited on nickel foam as a cathode material for the hydrogen evolution reaction and observed a reduced cathode overpotential in comparison with any other metal carbonate hydroxide [START_REF] Tang | Electronic and morphological dual modulation of cobalt carbonate hydroxides by Mn doping toward highly efficient and stable bifunctional electrocatalysts for overall water splitting[END_REF]. All these results were ob tained in highly concentrated alkaline solutions. This suggests that there is a parallel between water deprotonation and weak acid deprotonation since the same elements (Cr, Ni, Mo and Mn) seem to intensify the reaction rate in both cases. These elements could intervene in a key adsorption step, probably limiting, of the cata lytic mechanism.

Discussion

4.1. Reducing the cathode potential by optimizing the catalytic effect of carbonate on stainless steel

Concentration of the carbonate solution

The maximum hydrogen evolution on 254SMO SS was obtained with a highly concentrated carbonate solution (1 M and 2.8 M). At 10 A/m 2 , working with a 1 M carbonate solution as the cath olyte instead of domestic wastewater reduced E c by 380 mV, which represents a gain of 32% over the 1.2 V MEC limit voltage calculated in III.1.

In addition, using a highly concentrated carbonate solution instead of domestic wastewater as catholyte in an MEC also in creases the ionic conductivity in the cathode compartment very significantly (by a factor of 10). High conductivity facilitates ion transport in the solution and reduces the internal resistance of the MEC [START_REF] Oliot | Ion transport in microbial fuel cells: key roles, theory and critical review[END_REF].

Cathode material

Platinum is known to be a great heterogeneous electrocatalyst for hydrogen evolution but its high price discourages its use in industrial scale electrolysis cells. Nickel is also known to be an efficient catalyst for water reduction reaction and nickel based materials are viewed as a cheaper alternative to platinized cath odes [19,20,46e49]. For example, Raney nickel is a Ni Al based powder used as a catalyst in many industrial processes and espe cially in water electrolysis [START_REF] Damien | Hydrog ene par electrolyse de l'eau[END_REF]. Also, De Silva et al. have demon strated that SS AISI 316L is a more stable cathode material than platinum when associated with the homogeneous catalysis of wa ter reduction by phosphate [START_REF] Silva Munoz | Hydrogen production by electrolysis of a phosphate solution on a stainless steel cathode[END_REF].

316L and 254SMO SS proved to be the most appropriate to implement in the carbonate catalysis of water reduction. The choice of the MEC cathode material must be relevant in terms of perfor mance and also in terms of cost, availability and durability. 254SMO is a high quality SS, often used in media containing chlorides, e.g. marine environments [START_REF] Rajala | Influence of chlorination and choice of materials on fouling in cooling water system under brackish seawater conditions[END_REF] and HCl solutions [START_REF] De Micheli | Electrochemical behaviour of 254SMO stainless steel in comparison with 316L stainless steel and Hastelloy C276 in HCl media[END_REF], where the other grades of stainless steel are corroded. It is 30% more expensive, on average, than 316L, which is more commonly used in industry. The use of 316L as cathode material is thus a reasonable choice since it showed some of the best performance levels in terms of current density at 1.3 V vs. SCE, is relatively cheap and is easily available. Additional studies will be needed to validate its long term dura bility in MEC operating conditions. 4.2. Impact of using a concentrated carbonate solution as catholyte on the MEC design

Presence of a membrane

In an MEC in which the anode compartment contains domestic wastewater and the catholyte is a carbonate solution, a membrane is necessary to avoid mixing the two solutions. Since domestic wastewater has a relatively low ion concentration [START_REF] Verbanck | Major ion content of urban wastewater : assessment of per capita loading[END_REF], the use of a concentrated carbonate solution could create overpressure or overflow due to water osmosis from the wastewater to the car bonate solution. This issue may be easily solved by accepting some slight overpressure in the catholyte compartment. The best compromise should be found between improving the performance and minimizing technical constraints due to water osmosis.

pH balance within the cell

The water reduction reaction produces OH À ions and so the pH increases near the cathode during electrolysis. In the case of a single chamber separator less MEC operating in batch mode, the global pH is balanced because H þ ions are produced at the anode. In a dual chamber MEC, the separator limits the movement of ions within the cell. If a solution containing predominantly HCO 3 À ions is used as the catholyte, a cation exchange membrane should be used to prevent loss of these anions to the anolyte via diffusion and migration. In this case, electricity transport through the cell is mainly ensured by the motion of cations through the membrane, which do not contribute to the pH balance and result an enhanced pH increase in the cathode compartment [START_REF] Oliot | Ion transport in microbial fuel cells: key roles, theory and critical review[END_REF]. At 10 A/m 2 , the water reduction potential is significantly more negative at pH 11 and pH 12 than at pH between 7 and 10 (Fig. 2). To ensure stable and fast water reduction, the pH of the cathode compartment should be regulated to stay between 7 and 10.

A sustainable way to buffer the pH in the cathode compartment is to pass CO 2 gas into the catholyte. Bubbling gaseous CO 2 into an aqueous solution:

CO 2 þ H 2 O 4 H 2 CO 3 (CO 2 , H 2 O) (10)
shifts the pH towards lower values, according to the acid/base equilibrium (Equation ( 8)). Bubbling CO 2 into an MEC catholyte will lead to saturation of the carbonate solution. The precipitated car bonate salts can be recovered and reused [START_REF] Chang | Calcium carbonate precipitation for CO2 storage and utilization: a review of the carbonate crystallization and polymorphism[END_REF]. This is a way to store CO 2 in solid form, which may be less hazardous for the environ ment than gaseous storage in deep geological formations, with risk of CO 2 leakage [START_REF] Liu | State of the art review of the environmental assessment and risks of underground geo-energy resources exploitation[END_REF].

Conclusion

The catalytic effect of HCO 3 À for hydrogen evolution on stainless steel has been demonstrated and optimized. Using a 1 M potassium carbonate solution instead of domestic wastewater as the catholyte enabled the cell voltage to be reduced by 32%, at 10 A/m 2 , in rela tion to the calculated objective of 1.2 V. The water reduction reac tion reached its highest rate at pH 8 but the kinetics remained satisfactory in the pH range from 7 to 10. Testing various materials as cathodes highlighted that the electro catalytic effect of HCO 3 À was enhanced on 316L and 254SMO SS. These results are a first step towards the design of an energy neutral MEC for wastewater treatment. The next challenge is to associate a carbonate catalysed cathode with a microbial anode to evaluate and fix the possible issues connected with osmosis and pH drift.

Table 1

 1 Applied voltages and operating conditions reported in the literature for domestic-wastewater-fed MECs and abiotic industrial electrolysers for water electrolysis.

		Applied voltage	Current density	H 2 production	COD removal	Reference
	MECs	1 V	0.25 A/m 2	0.15 L/g-COD	45%	[28]
		0.41 V	0.47 A/m 2	0.045 L/d	87-100%	[29]
		1.1 V	0.27 A/m 2	0.088 L/d	Av 34%	[30]
		1 V	0.2 A/m 2	0.03 L/d	80%	[31]
		0.75 V	0.42 A/m 2		76%	[32]
		0.7 V	0.45 A/m 2		max 92%	[13]
		0.9 V	3.77 A/m 2	0.18 L/d	81-86%	[33]
		1 V	0.2 A/m 2		85%	[34]
		0.9 V	1.1 A/m 2		75%	[35]
	Industrial electrolysers	2 V a		5 Nm 3 /h	N/A	BPMP 5000 from Sagim-gip [36]
		2V a		12 Nm 3 /h	N/A	S18 MP from McPhy [37]
		2 V a		1100 Nm 3 /h	N/A	Electrolyser producing H 2 for the Ariane V rocket [38]

Acknowledgements

This work was funded by the French ANR within the framework of the WE MET project (ERANETMED 2015 European call).

Declaration of conflict of interest

None.