Frédéric Dufour

Gilles Pijaudier-Cabot

Marta Choinska

Antonio Huerta

Extraction of a crack opening from a continuous approach using regularized damage models

Crack opening governs many transfer properties that play a pivotal role in durability analyses. Instead of trying to combine continuum and discrete models in computational analyses, it would be attractive to derive from the continuum approach an estimate of crack opening, without considering the explicit description of a discontinuous displacement field in the computational model. This is the prime objective of this contribution.

The derivation is based on the comparison between two continuous variables: the distribution if the effective non local strain that controls damage and an analytical distribution of the effective non local variable that derives from a strong discontinuity analysis. Close to complete failure, these distributions should be very close to each other. Their comparison provides two quantities: the displacement jump across the crack and the distance between the two profiles. This distance is an error indicator defining how close the damage distribution is from that corresponding to a crack surrounded by a fracture process zone. It may subsequently serve in continuous/discrete models in order to define the threshold below which the continuum approach is close enough to the discrete one in order to switch descriptions. The estimation of the crack opening is illustrated on a onedimensional example and the error between the profiles issued from discontinuous and FE analyses is found to be of a few percents close to complete failure.

I. Introduction

For many concrete structures, crack opening is a key parameter needed in order to estimate durability. Cracks are preferential paths along which fluids or corrosive chemical species may penetrate inside concrete structural elements. For structures such as confinement vessels for instance, tightness to gas or liquids is a major serviceability criterion that is governed by Darcy's relation in which permeability of the material is involved.

The material permeability is strongly related to the amount of cracking in concrete: permeability grows significantly as distributed microcracking develops (see e.g. [START_REF] Choinska | Effects and interaction of temperature and stress level related damage on permeability of concrete[END_REF]) and it jumps several orders of magnitude upon macrocracking (Sugyiama et al. 1996, Hearn and[START_REF] Hearn | Measurement of Permeability under Uniaxial Compression-A Test Method[END_REF]. According to Poiseuille's law, the permeability of a cracked structure (with a single crack) is proportional to the square of the crack opening.

Hence, the prediction of the durability of structural components requires models that describe failure, crack locations and crack openings in the present example too when damage has localised.

• Enhanced continuum and integral damage models are capable of representing diffuse damage, crack initiation and possibly crack propagation [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF], [START_REF] Peerlings | Gradient enhanced damage for quasibrittle materials[END_REF]).

They regard cracking as an ultimate consequence of a gradual loss of material integrity. These models, however, do not predict crack opening as they rely on a continuum approach to fracture.

• Fictitious crack models are based on an explicit description of the discontinuity within the material (e.g. cohesive crack model of [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]). They relate the crack opening to the stress level and they are based on the linear elastic (or plastic) fracture mechanics. Cohesive crack models needs proper algorithms for crack propagation, and more importantly they are not capable of describing crack initiation.

Ideally, the prediction of durability that involves inception of failure, crack location, propagation and crack opening would require to merge the continuum damage approach and the discrete crack approach into a single, consistent, computational model bridging the continuous and discrete approaches. Bridges between damage and fracture have been devised in the literature (see e.g. [START_REF] Mazars | From damage to fracture mechanics and conversely: a combined approach[END_REF], [START_REF] Planas | Cohesive cracks versus nonlocal models: Closing the gap[END_REF]).

They rely on the equivalence between the dissipation of energy due to damage and the energy dissipated in order to propagate a crack. Given the energy dissipated in the damage process, the equivalent crack length is computed, knowing the fracture energy. Generally, the entire energy that is dissipated in the fracture process zone is "converted" into a crack length [START_REF] Mazars | From damage to fracture mechanics and conversely: a combined approach[END_REF]). Some part of this energy may be dissipated in the process zone outside from the crack and it follows that the crack length is probably overestimated, and the crack opening is overestimated too.

The strong discontinuity approach initiated by [START_REF] Simo | An analysis of strong discontinuities induced by strain softening in rate-independent inelastic solids[END_REF] and widely used over the last decade (e.g. [START_REF] Oliver | From continuum mechanics to fracture mechanics: the strong discontinuity approach[END_REF], [START_REF] Oliver | Computational modelling of cracking of concrete in strong discontinuity settings[END_REF], [START_REF] Larsson | Finite element embedded localization band for finite strain plasticity based on a regularized strong discontinuity[END_REF]) offers the possibility of merging in the same formulation a continuous damage model for the bulk response and a cohesive model for the discontinuous part of the kinematics. It is certainly a combination of continuum -discrete modelling that is sound from a theoretical point of view and appealing from the point of view of the physics of fracture. The issue in combining the continuum based model for crack initiation and then a discrete crack model for propagation is, however, the threshold upon which one switches from one analysis to the other. Usually, it is considered that the discontinuity appears when damage, stresses or strain energy reach a certain threshold fixed beforehand, which remains arbitrary [START_REF] Comi | An extended FE strategy for transition from continuum damage to mode I cohesive crack propagation[END_REF], [START_REF] Simone | From continuous to discontinuous failure in a gradient enhanced continuum damage model[END_REF]). As we will see further, one of the outcome of the present paper is to provide an indicator on the basis of which the appearance of a discontinuity during a damage process can be defined, with a given accuracy.

Instead of trying to combine continuum and discrete models in computational analyses, it would be attractive to derive from the continuum approach an estimate of crack opening, without considering the explicit description of a discontinuous displacement field in the computational model. This derivation could be based on some postprocessing of the distribution of strain and damage in the considered structure. The main purpose of this paper is to present such an estimate of crack opening derived from a continuum model description. First, we recall the continuum approaches that will be considered: the (integral) non local damage model and the gradient damage model. Nonlocal models are known to possess shortcomings such as spurious boundary effects on fracture propagation [START_REF] Jirasek | Size effect on fracture energy induced by nonlocality[END_REF]) or incorrect initiation of damage at a crack tip [START_REF] Simone | Incorrect initiation and propagation of failure in non-local and gradient-enhanced media[END_REF]). Still these defects in the model formulation do not alter their ability to capture a fully localised (mode I) crack. The estimate of crack opening is discussed in the second part. We start from an analytical expression of a regularized strong discontinuity and compare the induced strain field with results from the continuum based model. The comparison of the strain profiles provides the crack opening, and at the same time some indication of the quality of the estimate. Results obtained according to the integral and gradient damage models are illustrated on the one dimensional example of a tension bar.

II. Non local damage approach

Damage model

The scalar isotropic damage model due to Mazars [START_REF] Mazars | Continuum damage theory: application to concrete[END_REF]) will be used in the finite element computations for representing the progressive failure. In this model the stress (σ) -strain (ε) relationship is expressed as follows:

( ) ε σ : 1 C D - = ( 1 )
where D is the damage scalar variable and C is elastic stiffness of the material. Damage is a combination of two components:

which is damage due to tension based loads and which is damage due to compression: we shall be interested damage due to tension loads only, and then Eq. ( 2) reduces to:

( ) ( [ ) ] 0 1 1 0 D t Y Y B t t D e A Y A Y D - - - - = ( 3 )
where , and are model parameters and Y is defined by:

t A t B 0 D Y ( ) eq Y Y ε , max = ( 4 )
with initially, and the equivalent strain is defined as:

0 D Y Y = ∑ = + = 3 1 2 i i eq ε ε ( 5 ) + i
ε denotes the positive part of the principal strain i ε . This constitutive relation exhibits strain softening and two regularization techniques shall be considered in order to avoid mesh dependency and ill-posedness of the governing equations of equilibrium.

Integral non local model

In the non local damage model, the effective strain (Eq. ( 5)) is substituted by the average equivalent strain eq ε in the equations governing the growth of damage [START_REF] Pijaudier-Cabot | Nonlocal damage theory[END_REF]). The non local average is defined as

( ) ( ) ( ) ∫ ∫ Ω Ω - - = ds s x ds s s x eq eq φ ε φ ε ( 6 )
Several weight functions exist in the literature, we choose the most often used, i.e. the Gaussian function:

( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - = - 2 2 exp c G l s x s x φ ( 7 )
∑ = + = 3 1 2 i i eq ε ε ( 8 )
where i ε is a regularized principal strain that is solution of the equation

ε ε ε = ∇ -2 c ( 9
)
where c [m] provides a characteristic size to the model. Actually, this gradient enhanced regularization technique is a particular case of the integral approach with a specific, Green weight function that can be written in the 1D case [START_REF] Peerlings | A critical comparison of non local and gradient enhanced softening continua[END_REF]), as:

( ) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ - - = - c c GF l s x l s x exp 2 1 φ ( 1 0 )
One can observe that the slope at the centre (x=s) of the Green function GF φ is finite but non zero, whereas for the Gaussian function G φ this slope is zero. Therefore, the Green weight function is sharper at the centre, it decreases more rapidly providing less weight for neighbours of point x than the Gaussian weight function.

III. Estimation of crack opening

Let us consider the case of a 1D bar loaded in tension. The bar of length L is clamped at and a constant velocity v is applied at the other end ( ). We assume that the loading history is such that failure occurs at point (Fig. 1a) and we look at the ultimate failure state, after the separation of the bar in two pieces, now

0 = x L x = 0 x x =
assumed to be free of load. Note that for simplicity, the crack location is assumed to be known in our approach.

The displacement field will be constant piecewise with a step at the crack position 0 x x = .

The displacement profile (see Fig. 1b) is described according to a strong discontinuity formulation:

( ) [ ] ( 0 0 , x x H U x x u sd - = Γ ) ) ( 1 1 )
where subscript stands for strong discontinuity, is the Heaviside function and [ is the displacement jump across the crack. From Eq. ( 11), we may derive the strain field (see Fig. 1c) with the classical symmetric gradient operator :

sd Γ H ] U s ∇ [ ] ( 0 x x U u sd s sd - = ∇ = δ ε ( 1 2 ) where ( ) x δ
is the Dirac function. According to the above regularized damage models, discontinuous displacements may not necessarily exist, or may represent a limit case that is reached upon infinite strain. Hence, trying to compare the strain field in the bar to a discontinuous strain field may be a difficult task. An approximation may be to compare the profiles of the variable that controls the damage growth instead of comparing the strain itself. According to the non local model, this profile is continuous. 12) into Eq. ( 13), and using the properties of the Dirac function, we get:

[ ] ( ) [ ] ( ) ( ) ∫ Ω - - = ds s x x x U x x U sd φ φ ε 0 0 , , ( 1 4 )
which is a non local measure of the local displacement jump which has the shape of the weighting function. This error indicator defines how close the effective non local strain profile is from that corresponding to a discontinuous displacement distribution. It provides some information on the failure process and at the same time it is an indicator of accuracy of the crack opening estimate: if it is small, it means that complete failure is almost reached, and that the estimate of the crack opening is close to the real displacement discontinuity. This error indicator may subsequently serve in continuous/discrete models in order to define the threshold below which the continuum approach is close enough to the discrete one in order to switch descriptions.

] U Δ [ ] ( ) [ ] ( ) ( 0 0 0 , , , , , x x U x x U x x eq sd ) ε ε - = Δ ( 1
To summarize, we obtain from the comparison of the continuum based non local effective strain and the strong discontinuity based non local effective strain two quantities: (1) the displacement discontinuity (the crack opening), and (2) the error associated to the computation of this jump. The error is defined locally at each point of the bar, and globally after some averaging and normalisation. Note that in the absence of localised damage, the error in Eq. ( 17) is expected to be very high. This is just the consequence that distributed damage is the opposite of localised damage (and cracking) and in the first case, a crack opening does not really exist.

We are going now to apply this scheme to computations performed with the regularized damage models. The integral and gradient damage models will also be compared by means of their capacity to represent the discontinuity jump at failure accurately. In the comparisons, we have set

G φ φ =
in the derivations of the regularized strain fields according to the integral and gradient models. For the gradient model, the weight function is different. It is explicitly known in the 1D case, but it is not the case in a general 3D calculation.

Instead of computing the effective strain distribution that results from the Fredholm equation (Eq. 9) in which the local effective strain is a dirac δ function, we compute a new average, integral, effective strain from the local strain distribution obtained according to the gradient model with the weight function

IV. Numerical Examples

We solve numerically the model problem defined in Fig. 1a. The length of the bar is L=1m. The location of maximum strain and damage (crack position) is forced in the finite element located at the centre of the bar by setting a defect (smaller Young's modulus ) in this element. Damage will be first triggered in this element by a larger strain than elsewhere. The material parameters of the damage model are provided in Table 1 and Fig. 2 shows the material response in tension. for the large FPZ (LFPZ). The c parameter in the gradient enhanced approach was calibrated such that the integral and gradient models give the same FPZ width, that is, the same width of both the damage profiles. We obtained for the small FPZ and for the large one. Two finite element meshes are considered too: a coarse mesh of 31 elements and a fine mesh of 61 elements.

Integral and gradient damage responses

The global response is shown for the gradient model in Fig. 3a and for the non local model in Fig. 3b. These computations have been performed with the smallest value of the internal length in both non local models and for the fine and coarse meshes. We can observe that finite element results are not sensitive on the mesh size, which means that FE convergence with respect to the finite element size has been reached (for both large and small internal lengths). It is worth noting that due to the shape of the weight function (less sharp), the FE convergence is reached for a smaller number of finite element in the case of the integral approach. An arc-length technique has been used to reproduce such a global response which exhibits snap-back. Fig. 3c illustrates the differences between the gradient and integral approaches. Keeping the same size for the FPZ according to both models provides a difference between the load deflection curves. Conversely, [START_REF] Jason | Benchmarks for the validation of a non local damage model[END_REF] calibrated the gradient model in order to obtain the same load v.s. displacement response as for the integral model. In this case, c l c is approximately 2.55.

We plot damage and strain profiles in the post-peak regime for several ratios of the load to the peak load in the post-peak regime. Fig. 4 shows the corresponding evolution of damage profiles. As soon as damage starts to increase the width of the damaged zone is roughly twice the internal length (of the integral model) and it remains constant throughout the fracture process. Although the FPZ width has been calibrated to be the same for both regularization techniques, one can observe that the gradient model presents sharper peaks on the damage profiles compared to the integral approach. Fig. 5 shows the strain profiles along the bar that correspond to the damage profiles in Fig. 4. There is a consistent difference between the integral and the gradient approaches on these curves. The width of the strain profile is smaller according to the gradient approach than that of the integral approach. This is more clearly illustrated with the largest value of the internal length -comparing Fig. 5c and5d. The strain profile is sharper according to the gradient approach, and this provides the difference observed on the load -displacement responses obtained with the two models. These differences can be regarded also as the consequences of a sharper equivalent weight function according to the gradient approach compared to the Gaussian function used in the integral model. approaches and for the large FPZ with the gradient (c) and the integral (d) approaches.

Estimation of crack opening

We are now going to compare the computed regularized effective strain profiles with the profiles derived from the strong discontinuity analysis, derive the estimate of the crack opening and compute the distance between the two profiles which is a failure indicator, or equivalently an error indicator on the accuracy of the crack opening estimate. We consider first a state of damage which is close to complete failure. The ratio of the load to the peak load is 0.01 and the corresponding damage profiles are those drawn with a thick curve in Fig. 4. Figure 6 shows this comparison for the smallest internal length. By definition, the profiles coincide in the centre of the bar, the location with the highest value of damage and effective regularized strain, and the displacement jump [ ] U is computed from this condition. The profiles deduced from FE computations are wider, but overall, the strong discontinuity profiles and the computed profiles are rather close to each other. Figure 7 shows the estimates of the crack opening (Fig. 7a) and the evolution of the error indicator computed according to Eq. ( 17) in the course of the loading history. The load level considered in Fig. 6 corresponds to the black dot in Fig. 7. Once the crack is fully opened the increase of its opening is equal to the displacement applied at the active end of the bar as no further strain are stored in the material (see Fig. 7a). Therefore the plot of the displacement jump versus the imposed displacement tends to the bisecting line. The error indicator (see Fig. 7b) is quite large just after the peak load since the macro-crack is not yet clearly formed. If damage would have been uniformly distributed along the bar, the error would be even higher. Afterwards, the error indicator decreases rapidly to reach a limit value of 4% for the integral approach and 2% for the gradient regularization.

These values represent the level of ability of such regularized damage models to simulate discontinuous failure in terms of kinematic variables. They illustrate also the influence of the weight function. As pointed out in section 2, the major difference between the gradient and the integral approach is the weight function. With a sharper weight function, the estimate of the crack opening seems to be improved. Still, a limit value of the error of a few percent may be acceptable, provided some validation with experimental data is available. Fig. 8 shows that the evolution of the crack opening displacement and the quality of the crack opening estimates are not really dependent on the size of the FPZ. With a small FPZ (small internal length), snap-back on the load displacement response is more severe and limit values of the error are reached for a smaller applied displacement. Nevertheless, the results are quite similar. On this figure, the results obtained with the integral approach have been plotted. Similar observations could be made according to the gradient approach.

V. Conclusions

We have presented in this paper a post processing technique for the evaluation of a crack opening displacement, or its 1D equivalent, a displacement jump, from failure analyses based on regularized damage models. Instead of inserting into the kinematics of the problem the strong discontinuity and the inherent constitutive response, the estimate of the crack opening is based on the comparison of two strain distributions: the regularized effective strain that derives from a non local (integral or gradient) model and the regularized effective strain that is derived from a strong discontinuity that is placed at the expected crack location. By definition, these are two regularized continuous fields that can be approached with standard discretization techniques. If these distributions are close to each other, it means that the variable that control damage has reached a distribution that is close to that of a real crack and it is expected that the corresponding distribution of continuous damage is close to that resulting from a discontinuity (crack) in a regularised, non local, analysis. From the comparison between these two regularized fields, the amplitude of the displacement jump and an error indicator on these estimate, that defines also how far from complete failure the distribution of effective strain is, are computed. Close to complete failure, this post processing technique provides an estimate of the crack opening displacement from a finite element computation with an error of a few percent. Although the method is illustrated on a 1D example, it can be extended to 2D analysis too [START_REF] Legrain | Extraction of crack opening from a non local damage field[END_REF]).

The technique applies readily to the integral non local damage model where the regularization is explicit. In the gradient model, the calculation of the regularised strain distribution from the strong discontinuity displacement field is required. An alternative averaging technique -same as for the integral model -has been devised in order to compare effective strain distributions. The comparison between results obtained according to the non local integral and gradient damage models show that both approaches provide results that are similar. The quality of the crack opening estimate depends on the weight function that enters in the non local expression. The gradient approach is equivalent to a non local average with a sharper distribution, lending less weight to neighbouring points compared to the Gaussian distribution in the classical integral approach. As a consequence, the gradient approach provides better limit values of the quality for a formed crack than with the integral model.

Note that when damage is fully developed, the profiles of damage start to enlarge slightly according to the two models. This spurious feature has been already observed by Geers and co-workers (1998) who devised a possible remedy and restricted the non local averaging process to zones where damage grows. This modification of the non local damage model is expected to provide a better description of complete failure and should result into an error indicator on the quality of the crack opening estimate that is improved. This point remains to be investigated, along the quality of the crack opening estimate performed with existing models in which the internal length is changing in the course of damage, and of course validation using experimental data.

  α depend on the strain tensor. In the case of a uniaxial tensile

  where is the internal length of the model, a new model parameter compared to the classical -localof regularized model used in this paper follows the same principle, except that the nonlocal variable derives from a Helmholtz equation. The non local equivalent strain eq ε replaces its local counterpart. It is defined as:

  If it is close to represent a strong discontinuity, then the distribution of the non local effective strain to the distribution of effective strain that results from the displacement jump [ ] U .In order to perform this comparison, an effective strain distribution . 1d) are computed according to the same procedure as in the non local damage model.

Figure 1 :

 1 Figure 1: Ultimate state of failure of a 1D bar. (a) Loading system of with a default at

  amplitude of the discontinuity [ ] U . Imposing equality between both regularized effective strains distributions at the location of the crack (which is supposed to be known beforehand), one obtains: of the discontinuity -the crack opening -is derived from this condition (Eq. 15). One could have chosen another condition for computing [ ] U , for instance by minimization of the difference between both profiles. It will not change fundamentally the results if the effective non local distribution is close it should affect the numerical values of the displacement jump when the two nonlocal profiles are not close to each other. In other words, it will influence the values of the jump [ ] U prior to failure is reached. and from this comparison, an error indicator can be computed that measures how far from complete failure a given distribution of damage is, and as a consequence the confidence that one may have into the calculation of the displacement jump [ . First, an absolute error field is defined:

  6 )Then, an average relative error is evaluated by integrating the absolute error over the bar and normalizing it with the integral of the continuum non local strain:

  same weighting function. Therefore, the comparison of effective strains distributions that provide the estimate of crack opening is indirect compared to that in the integral formulation, but again, if the local strain distribution is discontinuous, the two averaged fields should be very close to each other.

  Material parameters for finite element computations.

Figure 2 :

 2 Figure 2: Stress-strain relation for the chosen set of material parameters

Figure 3 :

 3 Figure 3: Global response for the integral approach (a), for the gradient approach (b) for two meshes and for the smallest internal length. The global response (c) for the fine mesh, the two regularization techniques and the two FPZ sizes.

Figure 4 :

 4 Figure 4: Evolution of damage for the small FPZ with the gradient (a) and the integral (b) approaches, and for the large FPZ with the gradient (c) and the integral (d) approaches. Thick curves corresponds to damage profiles almost at complete failure, when the load is 1% of the peak load.

Figure 5 :

 5 Figure 5: Strain profiles evolution during loading for the small FPZ with the gradient (a) and the integral (b)

Figure 6 :

 6 Figure 6: Regularized effective strain profiles from the strong discontinuity approach (Dirac) and from FE computations using the integral technique (a) or the gradient technique (b) for the smallest FPZ width.

Figure 7 :

 7 Figure 7: Crack opening as a function of imposed displacement (a) or maximum strain (c) ; relative error as a function of imposed displacement (b) or maximum strain (d). Plots are given for the two regularization techniques with the smallest FPZ width. The black dot corresponds to a load level of 1% of the peak load in the softening regime.

Figure 8 :

 8 Figure 8: Crack opening as a function of imposed displacement (a) or maximum strain (c); relative error as a function of imposed displacement (b) or maximum strain (d). Plots are given for the integral damage model and for the two sizes of FPZ.

Acknowledgement: financial support from the Agence Nationale de la Recherche through the non thematic programme (contract NT05-4_43689) is gratefully acknowledged.


Morgane Hoareau

Luc Etcheverry

Benjamin Erable

Alain Bergel
email: alain.bergel@toulouse-inp.fr

Oxygen supply management to intensify wastewater treatment by a microbial electrochemical snorkel

Keywords:

The microbial electrochemical snorkel (MES) is a low-cost, low-maintenance technology that could considerably reduce wastewater treatment costs by reducing the need for aeration. An MES is a single electrode. The lower part, in the anoxic zone of the bioreactor, oxidizes organic matter by transferring electrons to the electrode. The upper part, in the oxic zone, releases the electrons by reducing dissolved oxygen. This study gives new insights into the correlation between the MES potential, the concentration of dissolved oxygen and COD abatement, thus enabling practical rules to be extracted for running the MES in optimal conditions, particularly for adjusting the aeration zone and frequency.

Here, aeration promoted the cathodic reaction and increased the MES potential to a maximum (0.0 to 0.1 V/SCE). After aeration stopped, the potential dropped to a minimum ( < -0.4 V/SCE). With 26 cm high MESs, sequential aeration at the top was not sufficient to support an effective cathodic zone, while aeration at the bottom was detrimental to the formation of the microbial anode. With MESs 48 cm high and aeration set up in the highest quarter, 30 minutes of aeration every 4.5 hours gave potential variations that were stable for weeks. The MESs continued to oxidize organic matter 30 minutes after the aeration had stopped. In this period, the MESs removed 3 to 6.3 times more COD than the control reactors. Increasing the electrode capacitance is therefore suggested as an effective way to further decrease the aeration cost by decreasing the aeration frequency.

Introduction

Wastewater treatment plants consume large amounts of energy, which represent a significant portion of municipalities', and thus countries', electricity consumption [START_REF] Fang | Linkage analysis for the water-energy nexus of city[END_REF] , reaching, for example, 1% of the total energy consumed by Germany or Italy [START_REF] Molinos-Senante | Energy intensity modeling for wastewater treatment technologies[END_REF] . In this framework, supplying oxygen is the most cost-intensive operation of the aerobic biological treatment plants. The microbial electrochemical snorkel (MES) is a low-cost, low-maintenance technology that could considerably reduce the cost of the oxygen supply by reducing or even avoiding the need for aeration of the bioreactor.

An MES can be defined simply as a single electrode in an environment where the composition of oxidizing and reducing species is spatially heterogeneous. In the context of wastewater treatment, the lower part of the MES is positioned in the lower, anoxic part of the bioreactor and ensures the oxidation of the organic matter [ 3 , 4 ]. Anaerobic bacteria spontaneously form an electroactive biofilm on the surface of that part of the MES and oxidize organic matter by transferring electrons to the electrode. This part of the MES forms the anode, which acts as the electron acceptor for the exoelectrogenic bacteria. The upper part of the MES, located in the aerobic upper zone of the bioreactor, releases the electrons from the electrode by reducing oxygen. It constitutes the cathode. Different types of organic matter can be oxidized at the anode part, such as volatile fatty acids and hydrocarbons [ 5 , 6 ], and other electron acceptors than oxygen, e.g. nitrates [START_REF] Yang | Denitrification of overlying water by microbial electrochemical snorkel[END_REF] , can be used at the cathode part. An MES can be considered as the simplest possible microbial fuel cell with short-circuited anode and cathode. As the anode and cathode are short-circuited, the MES does not produce any power but it supports the highest possible electrochemical reaction rates that an MFC could reach.

The first publication dedicated to an MES applied the concept to wastewater treatment [START_REF] Erable | From microbial fuel cell (MFC) to microbial electrochemical snorkel (MES): maximizing chemical oxygen demand (COD) removal from wastewater[END_REF] . The MES was compared with a microbial fuel cell (MFC) equipped with a resistance of 10 0 0 , a value commonly reported for MFCs. The chemical oxygen demand (COD) was abated by 40-50% with the MFC, while the abatement reached 60-75 % with the MES. MESs in the form of beds of conductive granules have given impressive results for wastewater treatment in constructed wetlands [START_REF] Aguirre-Sierra | A new concept in constructed wetlands: assessment of aerobic electroconductive biofilters[END_REF] . The concept is then named METland (Microbial Electrochemical Technology in constructed wetland). Compared to a conventional treatment lagoon equipped with https://doi.org/10.1016/j.electacta.2021.139103 0013-4686/ non-conductive gravel, a METland designed with conductive granules resulted in COD and BOD 5 reductions 3 and 4.5 times higher, respectively. The MES allowed the legal discharge thresholds of 125 mg.L -1 of COD and 25 mg.L -1 of BOD 5 to be reached with a considerably shorter hydraulic retention time (0.5 day) than the gravel-containing control (3.4 days) [START_REF] Aguirre-Sierra | Microbial electrochemical systems outperform fixed-bed biofilters in cleaning up urban wastewater[END_REF]. In the treatment of pig slurry, the COD reduction reached 90% with METland compared to 81% with a lagoon containing a sand bed [START_REF] Arias Ramírez-Vargas | Microbial community function in electroactive biofilm-based constructed wetlands[END_REF]. This process has also proven to be effective for nitrogen and phosphate treatment [START_REF] Ramírez-Vargas | Electroactive biofilm-based constructed wetland (EABB-CW): A mesocosmscale test of an innovative setup for wastewater treatment[END_REF].

Despite clearly promising successes, MESs, like MFCs, show performance levels that are severely limited by the low oxygen reduction rate at the cathode part [START_REF] Rimboud | Different methods used to form oxygen reducing biocathodes lead to different biomass quantities, bacterial communities, and electrochemical kinetics[END_REF][START_REF] Erable | Microbial Catalysis of the Oxygen Reduction Reaction for Microbial Fuel Cells: A Review[END_REF]. On the one hand, no catalyst is really efficient at pH values around neutrality [START_REF] Rimboud | Multiple electron transfer systems in oxygen reducing biocathodes revealed by different conditions of aeration/agitation[END_REF]. On the other hand, the concentration of dissolved oxygen is very low, around 9.09 mg/L at 20° C in water (0.28 mM), and its transport to the MES surface is consequently slow in the absence of stirring or forced aeration. It is essential to deal with the issue of oxygen supply management if the objective is to scale up MESs for wastewater treatment. It is necessary to find the delicate balance between increasing the mass transfer of oxygen to the MES and keeping the cost of aeration as low as possible.

The purpose of the present work was to optimize the oxygen supply in an MES. The challenge was to provide the MES with the amount of oxygen required for the extraction of electrons without hindering the oxidation reaction, which requires anoxic conditions. Two strategies for doing this were investigated.

In Strategy 1, the objective is to establish a stable profile of dissolved oxygen inside a bioreactor with a lower, anoxic zone and an upper, oxic one. The cathode develops in the upper zone, which is sequentially aerated, while, hopefully, the anoxic bottom is not reached by oxygen. In this way, it should be possible to achieve significant cost savings in wastewater treatment plants by aerating only the upper part of a bioreactor equipped with an MES.

In Strategy 2, the whole bioreactor is aerated sequentially, as is currently done in usual aerobic WW bioreactors. In this case, the whole electrode sequentially behaves as an anode, when it is not aerated, and as a cathode, when it is aerated. When placed in anoxic conditions, the electrode oxidizes organic matter by accepting electrons (anode). Conversely, when placed in oxic conditions, the electrode releases the electrons by reducing oxygen (cathode). So-called bidirectional microbial electrodes have been shown to be efficient for both oxygen reduction and organic matter oxidation [START_REF] Hoareau | Oxygen-reducing bidirectional microbial electrodes: A mini-review[END_REF][START_REF] Hoareau | Oxygen-reducing bidirectional microbial electrodes designed in real domestic wastewater[END_REF]. In Strategy 2, the hope is that the presence of the MES, acting as a bidirectional electrode, will reduce the frequency of aeration compared to current WW treatment bioreactors with no MES.

Besides checking these two strategies, this work demonstrated that the potential of the MES and the soluble oxygen concentration are closely related. Furthermore, the efficiency of the COD abatement by MES also depended on the potential. The MES potential is thus shown to be the key parameter that should be monitored in order to optimize the MES operating protocol. This conclusion has already been speculated theoretically [START_REF] Hoareau | Microbial electrochemical snorkels (MESs): A budding technology for multiple applications. A mini review[END_REF] but is demonstrated experimentally for the first time here. Furthermore, the monitoring of the potential allowed us to define the phase during which the MES is the most efficient in terms of COD abatement and thus gives a practical tool to optimize the operating parameters.

Materials and methods

Media and electrochemical set-up

Domestic wastewater (WW) was collected from the aerobic tank of a WW treatment plant (Castanet-Tolosan, near Toulouse, France) and used without any addition. When indicated, it was inoculated with 5 % (v/v) of activated sludge (AS) collected from the same aerobic tank. In some cases, when indicated, hydrolysed WW was used to stabilize the soluble COD (Chemical Oxygen Demand) content. Hydrolysis consisted of maintaining fresh WW under anaerobic conditions during one week at room temperature. This process increased the soluble COD from 150-250 mg O 2 .L -1 in fresh WW to 250-350 mg O 2 .L -1 in hydrolysed WW.

Two different MES configurations were used, called MES_26 and MES_48, because they were equipped with electrodes 26 cm and 48 cm high, respectively. Both set-ups were made of vertical glass cylinders of small diameter and large height in order to stabilize an aerobic upper zone and an anaerobic lower zone more easily (objective of Strategy 1). In addition, the reactors were filled with a large volume of electrode over their entire height ( Fig. 1 ) in order to approximate the situation of a liquid column in a bioreactor.

Each MES_26 electrode was made of two strips of carbon felt sheet (RVG 40 0 0, Mersen, France) 4 cm wide, 0.5 cm thick and 26 cm high, assembled in a cross shape with a titanium rod (1 mm in diameter) in the centre over the whole height ( Fig. 1 ). The MES_26 reactors were glass cylinders 4.2 cm in diameter, which contained WW to a height of 27 cm, i.e. the solution surface was 1 cm above the electrode.

The MES_48 electrodes were similar to the MES_26 ones but with a height of 48 cm and a width of 3 cm. The MES_48 reactors contained WW to a height of 50 cm, i.e. with 2 cm of solution above the electrode. The small reduction of the diameter and the significant increase in the height compared to the MES_26 were intended to increase the length of the upper cathode part, without disturbing the lower part, following the results obtained with the MES_26 (as explained at the beginning of part 3.2).

Each MES reactor was included in a closed hydraulic loop including a peristaltic pump (Masterflex®, Cole-Parmer, USA) and a 2 L storage tank. The recirculation was set up at 1.5 mL.min -1 . Forced aeration was achieved with air pumps, which bubbled air vigorously into the solution.

The titanium rod in the centre of the electrode was used to connect the MES to the electrical circuit. The potential of each MES was monitored vs. a saturated calomel reference electrode (SCE, Sentek, Australia, 0.24 V/SHE) using a multi-channel potentiostat (MPG, EClab software, Bio-Logic SA). All potentials were expressed vs. SCE.

Biofilm formation and operation

Preliminary oxygen concentration measurements were carried out in batch in an MES_26 reactor filled with WW, without the MES electrode inside. Forced aeration was set up on the top of the reactor, 2 cm from the free surface, for 3 hours.

The first experimental run implemented four MES_26 in parallel, with air bubbling set on the top in two of them (as in the preliminary experiment), and at the bottom in the others. After a stabilization period of 0.7 day without aeration, forced aeration sequencing was applied manually. After 1.6 day, automatic sequential aeration was set up: 2 hours with aeration and 4 hours without. The loop initially contained 2 L of WW inoculated with 5 % (v/v) AS. The storage tanks were first changed after 4 days and then replaced every 3 to 4 days with tanks containing 2L of only WW without AS.

The second experimental run implemented four MES_48 in parallel, two with an MES electrode and two controls with a plastic grid to replace the electrode. The control reactors were similar to the MES reactors and were equipped with the same cross structure with the same dimensions as the MES electrodes. Each of the two parts of the cross was composed of three plastic grids superimposed, so as to have a structure with the same thickness as the MES electrodes. For all 4 reactors, the aeration was set up in the upper quarter of the reactor, 12.5 cm from the free surface. At the beginning of the experiment, the aeration sequence was 2 hours with aeration and 4 hours without. At day 6, it was changed to 30 minutes with aeration and 5.5 hours without. The loop initially contained 2 L of WW inoculated with 5 % (v/v) AS. The storage tanks were changed on days 9 and 16 and replaced with tanks containing 2 L of WW only.

Three batch sessions were used to evaluate the removal of organic matter. For each batch, at the end of the 30 minute aeration period, the hydraulic recirculation pumps were stopped, the reactors were emptied and then refilled with around 300 mL of hydrolysed WW. The emptying and refilling operations were performed under a constant flux of nitrogen in order to avoid parasitic aerobic consumption of the COD. The organic matter removal was then monitored by taking successive samples from the middle of the MES reactor. COD removal rates ( COD) were calculated as:

COD ( % ) = C O D 0 -C O D t CO D 0 100 ( 1
)
where COD t is the COD value at time t and COD 0 the initial value of the first sampling just after medium change.

Analytical procedures

The concentration of dissolved oxygen was measured using selfadhesive oxygen-sensitive spots (SP-PSt3-SA23-D5-YOP-US, Pre-Sens, Germany), which were glued on the glass reactor wall at eights of 5, 10, 15, 20 and 25 cm. The reading was taken through a polymer optical fibre (POF-L2.5-2SMA, PreSens, Germany).

The soluble COD was measured after filtering the samples with a syringe filter (0.22 μm) using the Hach Lange LCK 314 system for concentrations ranging from 15 to 150 mg O 2 .L -1 and LCK 514 microkits for concentrations ranging from 100 to 2,0 0 0 mg O 2 .L -1 .

Results and discussion

MESs can consist of two different electrodes connected by an electric wire [START_REF] Hoareau | Microbial electrochemical snorkels (MESs): A budding technology for multiple applications. A mini review[END_REF] . This configuration allows a clear distinction to be made between anode and cathode, so the effect of the MES on COD removal can be accurately evaluated by comparing the two connected electrodes with the two unconnected electrodes, which then serve as the control.

MESs can also consist of a single piece of conductive material, such as a monolithic electrode or a bed of conductive granules [START_REF] Hoareau | Microbial electrochemical snorkels (MESs): A budding technology for multiple applications. A mini review[END_REF] . In this case, it is difficult to distinguish the anode and cathode zones. These zones can vary depending on oxygen transport in the solution. Moreover, it is difficult to define a relevant control setup because it is necessary to use a different material that is not conductive. For example, MESs in the form of fixed beds of carbon granules have been compared to gravel beds [START_REF] Aguirre-Sierra | Microbial electrochemical systems outperform fixed-bed biofilters in cleaning up urban wastewater[END_REF] or sand filters [START_REF] Ramírez-Vargas | Electroactive biofilm-based constructed wetland (EABB-CW): a mesocosmscale test of an innovative setup for wastewater treatment[END_REF] , which are then used as the control.

Nevertheless, the single piece configuration was implemented here because it is closer to the configuration that could be developed at large size in treatment plants. Here, the control MESs were designed to be as close as possible to the MES conductive set-up by using plastic grids of dimensions identical to those of the carbon electrodes.

Oxic/anoxic profile (Strategy 1) vs. bidirectional electrode (Strategy 2)

Preliminary measurements of the dissolved oxygen concentration profile were carried out in an MES_26 reactor, filled with 27 cm of WW without the electrode inside. Measurements were performed in batch with continuous aeration on the top of the reactor, 2 cm from the free surface. Oxygen concentration reached 6.4 mg.L -1 in 20 minutes at the top of the reactor (height 25 cm). Then the value remained fairly stable at the top with 6.2 ± 0.6 mg.L -1 for more than 2 hours ( Fig. 2 ). At 20 cm height, the oxygen concentration rose to a maximum of 3.1 mg.L -1 at 30 minutes but then dropped again to 0 until the end. All the other measurements at 15, 10 and 5 cm height showed a nil concentration during the 3 hours.

The aeration system used here made it possible to create a stable profile of dissolved oxygen with an upper oxic zone and a lower anoxic part in the MES_26 reactor. It was consequently decided to run two MES_26 with this aeration system, in order to assess Strategy 1 (see Introduction), and two other MES_26 with aeration at the bottom of the reactor to compare with Strategy 2. In both cases, the objective was to optimize the sequence of aeration.

After 0.7 day under constant hydraulic recirculation and without any aeration, the potentials of the four MES_26, which were run in parallel, were stable at -0.42 ±0.01 V/SCE ( Fig. 3 A). The oxygen concentration was monitored in one reactor that was aerated at the bottom. After 0.7 day, the concentration of dissolved oxygen was nil. As the measurement was made at the top of the reactor (25 cm high), 2 cm under the free surface of the solution, this implies that the oxygen concentration was nil in the whole reactor. When the aeration was started at 0.7 day, the potential increased rapidly, reaching -0.13 ±0.01 V/SCE within 10 minutes, and then stayed stable at -0.05 ± 0.07 V/SCE during the aeration period ( Fig. 3 B). The oxygen concentration rose from 0 to 5.8 mg.L -1 within 10 minutes and then stabilized at 7.31 ± 0.1 mg.L -1 during the aeration period. Aeration was then stopped and the oxygen As explained elsewhere [START_REF] Hoareau | Microbial electrochemical snorkels (MESs): A budding technology for multiple applications. A mini review[END_REF] , the potential of an MES is a mixed potential, which results from the large potential gradient along the electrode, from the anode part at the bottom to the cathode part at the top. Here, the variation from the minimal value below -0.4 V/SCE in anoxic conditions, to the maximal value close to 0.0 V/SCE during the forced aeration period, means that the MES was driven by its anode or cathode function, respectively.

Studies on bidirectional microbial electrodes have observed similar potential variations according to aeration conditions with a high potential in the aerated cathodic phase and a low potential in the anoxic anodic phase. A similar amplitude, with a potential ranging from -405 mV/SCE during the anodic phase to 195 mV/SCE during the cathodic phase, has been observed in a solar energy powered MFC [START_REF] Strik | Solar Energy Powered Microbial Fuel Cell with a Reversible Bioelectrode[END_REF] . Greater amplitude, ranging from -450 to 300 mV/SCE, has been reported for dye treatment with an MFC [START_REF] Sun | Sequential decolorization of azo dye and mineralization of decolorization liquid coupled with bioelectricity generation using a pH self-neutralized photobioelectrochemical system operated with polarity reversion[END_REF] . Larger amplitudes of potential can be measured in MFCs because the anode and the cathode are separated, so the potential is measured on an electrode that functions strictly as an anode or strictly as a cathode. In contrast, the potential of an MES is most often a mixed potential.

The objective for WW treatment is that the MES should operate at high potential in order to promote the oxidation of organic matter. The aeration periods must drive the potential to high values. Then the non-aerated periods must not leave the potential at low values for too long because, at low potential values, the MES is no longer active for the oxidation of organic matter.

With this rule in mind, the first two aeration periods, which were managed manually, showed that an aeration period of 2 hours, coupled with a period of 4 hours without aeration seemed to be a suitable sequence. Accordingly, the automatic aeration was started after 39 hours (1.6 days) with this sequence and maintained until the end of the experiment. The maximal and minimal values of potential and dissolved oxygen concentration showed fair stability during the 14 cycles ( Fig. 3 A). The MES worked in an appropriate way, reaching high values of potential during forced aeration and not staying at minimal potential values for too long.

As described for the first aeration period ( Fig. 3 B), the potential response was always delayed with respect to oxygen concentration when the aeration stopped. For instance, during the third aeration period ( Fig. 3 C) the oxygen concentration dropped to 0 within 1 hour, while the potential dropped from -0.01 ±0.01 V/SCE to -0.43 ±0.03 V/SCE within 2.5 hours. The potential of the MES continued to decrease, when there was no longer any oxygen in solution. This showed that the MES remained able to accept electrons from organic matter even when there was no longer oxygen in solution to release them. This capacity of the MES is important because it makes it possible to lengthen the period without aeration, and thus to minimize the aeration cost. Increasing the capacitance of the MES should be a key objective for future research in order to decrease the frequency of the forced aeration phases.

The potential of the MESs monitored for 15 days showed some changes of the value of the maximal potential just after the storage tank replacement ( Fig. 4 ). For the MESs with aeration at the top, after the tank change, the maximal potential dropped by 0.10 ±0.03 V/SCE compared to the values before the change. For the MESs with aeration at the bottom, it only dropped by 0.04 ±0.02 V/SCE ( Table 1 ).

After each tank change, the MESs with aeration at the top had more and more difficulty in recovering high maximal values. For instance, after the second tank replacement, about 6 aeration cycles were necessary to recover a maximal potential around 0.0 V/SCE, whereas 9 were necessary after the third tank replacement ( Fig. 4 ). This means that the cathode had more and more difficulty in raising the potential to high values when the MES was supplied with organic matter. This was due to the progressive improvement of the anode part, which became increasingly efficient and pulled the potential down when the MES was provided with organic matter. Such evolution was not observed with the MESs aerated at the bottom, probably because the oxic conditions throughout the reactor hindered the formation of an efficient microbial anode.

It should be noted that the literature reports several bidirectional microbial electrodes [START_REF] Hoareau | Oxygen-reducing bidirectional microbial electrodes: A mini-review[END_REF] that have shown efficient anode properties and accepted full exposure to aerated conditions. However, all these bidirectional microbial electrodes were formed with an initial anode phase of several days, sometimes several weeks, performed under anaerobic conditions. Here, in the MESs aerated at the bottom, the phases of full aeration started from the beginning of the experiment. Such early aeration phases were certainly not appropriate to develop efficient microbial anodes. Aeration from the bottom could remain a possible alternative after a first phase of microbial anode formation, which should be performed with aeration at the top. Nevertheless, this would make the set-up and protocol much more complex, particularly if the development of full sized MESs in WW treatment bioreactors is envisioned. This led us to abandon Strategy 2 based on full aeration of the MES.

Domestic wastewater treatment

Strategy 1 was pursued with longer MESs (MES_48). The aeration zone was enlarged in order to obtain stronger cathodes, able to balance even fully developed anode parts. Actually, the results obtained with the MES_26 showed that the cathodes formed by aeration 2 cm below the free surface were no longer efficient enough to drive the potential to high values when the anode became more efficient. The MES_48 were consequently aerated at the upper quarter of the reactor, 12.5 cm from the free surface.

Until day 6, the aeration sequence included 2 hours with aeration and 4 hours without, based on the results of MES_26. After day 6, the aeration period was decreased to 30 minutes as this was sufficient to reach the high potential. The cycles then included a period of 30 minutes with aeration and 4.5 hours without. The potential showed the same general trend as with the MES_26, with fair stability ( Fig. 5 ). A slight drop in potential after the storage tank changes showed that the anode parts were efficient. Nevertheless, the cathode parts managed to maintain the potential above 0.0 V/SCE almost all the time, which showed that enlarging the aeration zone had a positive effect.

At the end of the experiment, oxygen concentration was measured as a function of time at different heights of the MES. At the end of the period without aeration, the oxygen concentration was zero in the whole reactor. After 30 minutes of aeration, the concentration remained zero at the bottom of the MES, but had increased in the middle and reached 6.7 mg.L -1 from 25 cm high to the top ( Fig. 6 ). It can be concluded that the 30 minutes of aeration managed to establish the desired profile of dissolved oxygen concentration, with an upper, fully oxic zone and a lower zone remaining anoxic.

MES_48s showed high stability during the 27 days of operation ( Fig. 5 ). The aeration periods allowed the MESs to reach maximal potential between 0 V/SCE, immediately after the tank changes, and + 0.1 V/SCE when the substrate was low. The minimal potential of -0.45 ±0.07 V/SCE was very stable during the anoxic phases. Even the pauses in aeration and recirculation, which were inserted to monitor the COD abatement in batch, did not cause any variation in the maximal and minimal potential values, which returned to the same values at the end of the batch phases.

Batch phases were set just after the changes of the storage tanks in order to monitor COD abatement in batch conditions. In the first batch phase, the initial COD measurement ( DCO 0 in equation 1 ) was made just after the tank change (t = 1 hour in Fig. 7 A) and the second 5 hours later (t = 6 hours in Fig. 7 A). COD removal was 2 ±4 % in the controls equipped with plastic structures and 32 ±0 % in the MES_48. The snorkel effect was thus obvious. Nevertheless, COD abatement did not increase afterwards. On the contrary, it tended to decrease, with a value of 25 ±0.3 % in the MES_48 after 24 hours (less than 4 ±6 % in the control reactors). This apparent decrease of COD removal may have been due to the solubilization of solid COD or release of organic matter from the biofilm. Looking at the COD abatement in parallel with the potential evolution ( Fig. 7 A) reveals that there was no further abatement from the time the minimal potential was reached. The second batch therefore focused on this 4 hour period of time and the third batch focused on the first hour. Batch 2 showed that there was no longer any COD abatement 2 hours after aeration was stopped. At this time, the potential just recovered its minimal value ( Fig. 7 B). Batches 2 and 3 gave similar COD removal rates 1 hour after the aeration was stopped (2 ±1 % and 12 ±2 % in the controls and the MES_48, respectively). Batch 3 showed that the period of the main COD abatement was about 30 minutes after the aeration was stopped. During these 30 minutes, the potential had completed about 90% of its drop from its maximal to its minimal value ( Fig. 7 C). The potential drop was directly correlated to the MES activity.

Physically, the electrode is emptied of electrons during the aeration phase and thus reaches its maximal potential value. When the aeration phase is stopped, the electrode collects electrons by oxidizing organic matter. The MES is thus active and abates the The three batches showed the efficiency of the MES_48. COD removal rates were always significantly higher than in the control reactors equipped with a plastic structure. For Batch 1, the MESs removed 6.3 times more COD than the controls did (25 ±0.3 % vs 4 ±6 %). The result was slightly lower for Batch 2, with 5.7 times greater removal rates (17 ±2 % vs 3 ±2 % 4 hours after aeration was stopped). Batch 3 led to 3 times greater removal rates with the MES (12 ±1 % vs 3 ±0.1 %) 1 hour after aeration was stopped.

These COD abatement rates seem modest compared to some high values reported in the literature, for instance, the 60-75% obtained with an MES for wastewater treatment [START_REF] Erable | From microbial fuel cell (MFC) to microbial electrochemical snorkel (MES): maximizing chemical oxygen demand (COD) removal from wastewater[END_REF] . However, the data presented here are very encouraging considering that the MES_48 worked for only around 2 hours after the cessation of aeration and achieved the majority of the COD abatement in only 30 minutes. Here, the objective was to advance in deciphering the MES behaviour. For this reason, the COD abatement was monitored in batch, during the period of non-aeration, in order to really identify the MES effect and to distinguish it from the aeration effect. To be fair, the comparison had to be made between the same reactors, operated in the same way and equipped or not with an MES. This comparison showed that the MES multiplied the removal rates by a factor of 3 to 6.3 depending on the batch and the time it lasted.

The approach taken here affords useful information for progress to be made in optimizing the operating conditions of MESs. First of all, it should be emphasized that measuring the MES potential must now be considered as the essential parameter to control the MES operation. A clear correlation was pointed out here between the MES potential, oxygen concentration and MES activity. This correlation now gives clear practical rules for running MES reactors in conditions that are as close as possible to optimal, including optimization of the aeration frequency and of the aeration location. These practical rules can be summarized as follows:

-Aeration must pull the MES potential to its maximal value of about 0.1 V/SCE and it can be stopped when this potential is reached. -During the non-aeration periods, the MES is active for organic matter oxidation as long as the minimal value of the potential is not reached. When the minimal value is reached (-0.4 to -0.5 V/SCE), aeration must be restarted. -Aeration must be performed in the upper zone of the MES and keep an anoxic zone at the bottom (Strategy 1). Actually, aeration of the whole MES length (Strategy 2) is detrimental to the initial formation of the microbial anode. -A significant decrease of the value of maximal potential when organic matter is supplied to the MES indicates that the cathode is not strong enough to drive the potential towards high values and that the upper aeration zone must be enlarged.

In addition, it was pointed out that the activity of the MES for organic matter oxidation during the non-aerated periods depends on its capacity to accumulate electrons. Consequently, the more capacitive the electrode is, the more capable it will be of storing electrons, enabling longer non-aerated periods. Increasing the capacitance of the MES material should consequently be a way to space out the aeration periods, and thus to decrease the aeration cost. The interest of capacitive materials has been widely investigated in the domain of microbial fuel cells [START_REF] Deeke | Influence of the thickness of the capacitive layer on the performance of bioan-odes in Microbial Fuel Cells[END_REF][START_REF] Borsje | Performance of single carbon granules as perspective for larger scale capacitive bioanodes[END_REF][START_REF] Li | Surface modification of microbial fuel cells anodes: approaches to practical design[END_REF][START_REF] Santoro | Threedimensional graphene nanosheets as cathode catalysts in standard and supercapacitive microbial fuel cell[END_REF][START_REF] Caizán-Juanarena | Combination of bioelectrochemical systems and electrochemical capacitors: Principles, analysis and opportunities[END_REF] . The considerable experience gained in this context should be very useful in improving MESs by reducing the frequency of aeration.

Conclusion

Potential has been demonstrated to be the key parameter that should be monitored to manage an MES here, particularly when setting the aeration zone and controlling the aeration frequency. To the best of our knowledge, this is the first time that a correlation has been shown between the MES potential and the dissolved oxygen concentration. Practical rules have thus been stated for using the potential measurement to control the aeration frequency and set the aeration location in the MES reactors. The period of activity of the MES after the end of the aeration was about 30 minutes here. During these 30 minutes, the MESs abated 3 to 6.3 times more COD than the controls. This active period without aeration corresponds to the time required for the electrode to reach its maximal level of accumulation of the electrons coming from the oxidation of organic matter. Increasing the capacitance of the MES material should therefore be a research priority to improve the system further by decreasing the aeration frequency.

Fig. 1 .

 1 Fig. 1. Diagram of MES reactor and MES electrode

Fig. 2 .

 2 Fig. 2. Concentration of dissolved oxygen as a function of time at different heights (from the bottom) of an MES_26. Batch mode, aeration 2 cm below the free surface.

Fig. 3 .

 3 Fig. 3. Potential and oxygen concentration versus time in an MES_26 with aeration at the bottom. (A) first four days with switch to automatic aeration at day 1.8, (B) focus on the first aeration period with manual control, and (C) focus on the first aeration period with automatic control; the black arrow indicates replacement of the storage tanks.

Fig. 4 .

 4 Fig. 4. Potential of the four MES_26 versus time. Two MES_26 had aeration at the top (blue) and two at the bottom (orange); black arrows indicate the tank changes

Fig. 5 .

 5 Fig. 5. Potential of the two MES_48 versus time. On day 6 (dashed line), the aeration sequence changed from 2 hours with air bubbling / 4 hours without to 30 minutes with / 5.5 hours without. COD monitoring in batch took place at days 7, 15 and 27 (dashed arrows).

Fig. 6 .

 6 Fig. 6. Oxygen concentration versus time at different heights (from the bottom) of an MES_48 during the aeration period

Fig. 7 .

 7 Fig. 7. Potential and average COD removal rate in the MES_48 and controls versus time in batches 1, 2 and 3 (A, B and C respectively). Aeration started at 0.5 hours (arrow up) and stopped 30 minutes later (arrow down). The medium was then changed, recirculation was stopped and COD measurements began.

Table 1

 1 Maximal potential before and after tank changes

	Moment	Measurement	Tank change	Aeration at top	Aeration at the bottom
	Max.	Before	1	0.02 ± 0.01	0.04 ± 0.01
	Po-	tank	2	0.04 ± 0.03	0.04 ± 0.01
	ten-	change	3	0.03 ± 0.02	0.04 ± 0.01
	tial	After	1	-0.12 ± 0.04	-0.02 ± 0.01
	(V/SCE)	tank	2	-0.04 ± 0.01	0.01 ± 0.01
		change	3	-0.06 ± 0.01	0.02 ± 0.01

Acknowledgements

This work was supported by the French Agence Nationale de la Recherche (ANR), within the framework of the Biotuba project ( ANR-17-CE06-0015 ).

Author's role

Morgane HOAREAU: Conceptualization, Methodology, Investigation, Writing -Original Draft.

Luc ETCHEVERRY: Conceptualization, Methodology.

Benjamin ERABLE: Conceptualization, Methodology, Validation, Supervision.

Alain BERGEL: Conceptualization, Validation, Writing -Original Draft, Writing -Review & Editing, Supervision, Project.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Theocharis Baxevanis
email: theocharis@tem.uoc.gr

Gilles Pijaudier-Cabot
email: gilles.pijaudier-cabot@ec-nantes.fr

Frédéric Dufour
email: frederic.dufour@ec-nantes.fr

Bifurcation and creep effects in a viscoelastic non-local damageable continuum

Keywords: Bifurcation, non-local constitutive laws, size effects, creep, damage

The conditions for localization in a material described by a non-local damage-based constitutive relation coupled with a Kelvin type creep relation are derived in a closed form. The inception of a localized mode is considered as a bifurcation into

Introduction

The growth of micro-cracks in progressively fracturing rate-independent materials is a fairly distributed process which is associated with stable material response. For certain stress trajectories, however, a different deformation mode may prevail, consisting of formation of discrete failure planes (macro-cracks).

In the latter case, the mechanical response as observed at the macroscale becomes, in general, unstable. The inception of a localized mode may be considered as a bifurcation problem related to the loss of positive definiteness of the tangent material stiffness operator governing the homogeneous deformation [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF][START_REF] Rice | The localization of plastic deformation, theoritical and applied mechanics[END_REF]. In mathematical terms, the elliptic character of the set of partial differential equations governing equilibrium is lost (ill-posedness of the related boundary value problem). This result was derived for the linearized rate equation problem, considering what is commonly denoted as a linear comparison solid [START_REF] Hill | Some basic principles in the mechanics of solids without natural time[END_REF]. Loss of ellipticity may result in a discontinuous rate of deformation and constant total energy consumption.

Path stability considerations indicate that this particular solution is expected and consequently failure occurs without energy dissipation [START_REF] Bažant | Stable states and paths of structures with plasticity or damage[END_REF]. In dynamics, the situation is similar. The hyperbolic differential equations of motion may become parabolic or elliptic and consequently the initial-boundary value problem ill-posed [START_REF] Hadamart | Leçons sur la Propagation des Ondes[END_REF].

Non-conventional constitutive relations, called localization limiters, based on sophisticated techniques that enrich the standard continuum, have been adopted by many authors to eliminate the deficiencies caused by ill-posedness of boundary value problems. A wide class of localization limiters is based on the assumption that the stress of a material point is not only determined by the history of strain at this point but also by interactions with other material points (non-local integral continuum). The distance over which interactions are important is related to an internal length that prevents energy dissipation to become zero as failure occurs [START_REF] Simo | On a stress resultant geometrically exact shell model, part i, formulation and optimal parametrization[END_REF][START_REF] Bažant | Nonlocal continuum damage, localization instability and convergence[END_REF]. Another class of localization limiters relies on enrichment by terms that contain gradients of the state variables. Some of these models deal with strain gradients, i.e. higher-order gradients of the displacement field. In the above theories the rate-independent nature of the material model is preserved. Incorporation of rate-dependent viscous terms can limit localization due to softening [START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF][START_REF] Loret | Dynamic strain localization in fluid saturated porous media[END_REF][START_REF] Dubé | Rate dependent damage model for concrete in dynamics[END_REF]. The size of the localization zone, however, is controlled by this rate effect and it may not be possible to fit, with the same expression of the evolution of damage or viscoplastic strains, the rate effect on the stress-strain response and the width of the localization zone that is consistent with experimental observations. In this contribution, attention is focused on concrete structures and we are interested in predicting the failure due to localization of damage of those structures. In common practice, it is usually assumed for concrete that lin-ear viscoelasticity takes place for low stress levels where the instantaneous mechanical behavior is elastic. In the contrary, for high stress levels the linearity hypothesis fails. Bazant and co-workers [START_REF] Bažant | Current status and advances in the theory of creep and interaction with fracture[END_REF] have shown that creep at high stress levels cannot be captured merely by a non-linear generalization of viscoelastic stress-strain relation. Experimental evidence associate creep strains at high stress levels to microcracking nucleation and growth, i.e. to damage (Rüsch was the first one to conduct such experiments [START_REF] Rüsch | Research toward a general flexural theory for structural concrete[END_REF]). As recalled by Bazant and Planas [START_REF] Bažant | Fracture and size Effect in Concrete and Other Quasi-Brittle Materials[END_REF], timedependent fracture of concrete is caused by viscoelasticity of the material, and (or) breakage of bonds in the fracture process zone. The difference between these two effects can be observed on structural responses and size effect. In the range of quasi-static loading, that is, in the absence of inertia forces and wave propagation effects, viscoelasticity in the bulk (linear creep) causes an increase of the bearing capacity for increasing loading rates, with a decreasing peak displacement. For geometrically similar structures of different sizes and mode I crack propagation, the set of data points corresponding to the nominal stresses at failure shifts toward increasing brittleness on the size effect plot [START_REF] Bažant | Fracture and size Effect in Concrete and Other Quasi-Brittle Materials[END_REF] The size of the fracture process zone decreases (see e.g. [START_REF] Haidar | Correlation between internal length, fracture process zone and size effect in mortar and model materials[END_REF]). This shift has been observed experimentally by Bazant and Gettu [START_REF] Bažant | Rate effects and load relaxation: static fracture of concrete[END_REF][START_REF] Bažant | Cohesive crack with rate-dependent opening and viscoelasticity: I. mathematical model and scaling[END_REF] Many writers have used coupled models to capture mechanical, viscoelastic and non-linear instantaneous behaviors. In these models, linear viscoelasticity is coupled to a rate independent elasto-plastic model [START_REF] De Borst | Computational issues in time-dependent deformation and fracture of concrete[END_REF], to a smeared crack model [START_REF] Rots | Computational Modelling of Concrete Fracture[END_REF] or to a damage model [START_REF] Mazzotti | Nonlinear creep damage model for concrete under uniaxial compression[END_REF]. The aim of this paper is to investigate the condition of localization for such a coupled model and to see the influence of creep on failure.

Since the influence of viscoelasticity on the strength of concrete structures may be of importance (decrease of bearing capacity, increase of brittleness) it is important for engineering practice to check that available constitutive models are capable of predicting such effects.

In the present study we will consider as an example a viscoelastic model coupled to rate independent damage inspired from Omar et al. [START_REF] Omar | Etude comparative du couplage endommagement -fluage[END_REF], in which the relationship between the effective stress, defined in a standard way according to continuum damage [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF], and the strain follows a Kelvin chain. The analysis follows the method used by Benallal [START_REF]Localization phenomena in thermo-elastoplasticity[END_REF] for thermo-mechanical problems and Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Strain localization and bifurcation in a nonlocal continuum[END_REF] for the case of a rate independent non-local model. This paper is organized as follows: In Section 2, the creep and damage models are briefly recalled and the rate formulation for creep-damage model is derived based on the small strain assumption. In Section 3, the strain localization analysis is performed. The rate constitutive relations for the linear comparison solid, as well as the bifurcation conditions under the assumption that the solid considered is either infinite or at least large enough so that boundary layers effects can be neglected, are derived.

Moreover, the dependence of the admissible wavelength on the parameters of constitutive law and on the state variables is discussed. An illustration is presented in the case of a uniaxial response in Section 4.

Constitutive relation

Creep model

Consider first a quasi-brittle material, such as concrete, exhibiting non-ageing linear viscoelasticity and characterized by the uniaxial compliance function J(t, t ), representing the uniaxial strain at age t caused by a unit stress enforced at any age t (the linearity hypothesis agrees very well with test results in which no strain reversals have taken place, in the case of basic creep, and for stress levels less than 40% of the strength limit). The assumption of a non-ageing material will be used throughout this paper for the sake of simplicity, and causes J(t, t ) to be a function of just the time lag (t -t ).

By approximating the compliance function by a Dirichlet series

J(t, t ) 1 E + N i=1 1 E i 1 -exp - t -t τ i , (1)
where τ i , i = 1, 2, ....N , are fixed parameters called retardation times, E is

Young's modulus and E i , i = 1, 2, ...N , are age-independent moduli which can be determined by least-square fitting to the "exact" compliance function, it can be proved that

tot (t) = el,d (t) + N i=1 cr i (t), i = 1, 2, ...N, (2)
where the strains el,d and cr i , i = 1, 2, ...N, are governed by the following equations

E el,d (t) = σ(t), (3)
and

E i cr i (t) + τ i E i ˙ cr i (t) = σ(t), i = 1, 2, ...N. (4)
The total stress σ in ( 4) is expressed as the sum of two terms. The first term,

E i cr i
, corresponds to the stress of an elastic spring of stiffness E i . The second term, τ i E i ˙ cr i , is the stress generated by the strain rate ˙ cr i in a linear dashpot of viscosity η i = τ i E i . This means that in the non-ageing case the approximation of the compliance function by the Dirichlet series corresponds to a Kelvin chain with constant properties of individual elements. The initial conditions for cr i are cr i (0) = 0.

Using (2), ( 3) and ( 4) , we finally obtain the rate form

σ = E ˙ tot -E N i=1 σ τ i E i - cr i τ i . (5)
A generalization of this uniaxial stress-strain constitutive relation is performed under the simplifying assumption that the Poisson ratio for creep is approximately time-independent and about the same as the elastic Poisson ratio ν.

Thus, in view of the assumed isotropy and based on the principle of superposition, the following relation holds:

ε tot (t) = C -1 ν J(t, t )σ(t ) + C -1 ν t t J(t, τ )dσ(τ ), (6)
where σ(t) is the column matrix of stress at time t, ε tot (t) is the column matrix of engineering strain components, and

C -1 ν =                    1 -ν -ν 0 0 0 1 -ν 0 0 0 1 0 0 0 2 (1 + ν) 0 0 sym. 2 (1 + ν) 0 2 (1 + ν)                   
, is the unit elastic compliance matrix, corresponding to Young's modulus E = 1 (for a review see [START_REF] Bažant | Current status and advances in the theory of creep and interaction with fracture[END_REF]) and leads to the following incremental form

σ = EC ν εtot + EC ν N i=1 ε cr i τ i -E N i=1 σ τ i E i . ( 7
)

Non-local damage model

We will use in the following the scalar continuum damage model proposed by Pijaudier-Cabot and Bazant (Pijaudier-Cabot and Bažant, 1987). The description of the model is that of [START_REF] Pijaudier-Cabot | Strain localization and bifurcation in a nonlocal continuum[END_REF]. Despite its simplicity (namely the fact that the model response is the same in tension and compression), this model bears the essential characteristics pertaining to a non-local integral model. The constitutive relation reads:

σ = (1 -D)EC ν ε el,d , (8)
where D is the damage variable and ε el,d is the (elastic) strain.

The growth of damage is defined by a loading function f

f (ȳ, D) = ȳ 0 F (z)dz -D, (9)
where F is a function deduced from experimental data and ȳ(x) is the average energy release rate due to damage at point x of the solid

ȳ(x) = V ψ(x -s)y(s)ds. (10)
V is the volume of the structure, ψ(xs) is a normalized weighting function and y(s) is the energy release rate due to damage at point s defined by

y(s) = 1 2 ε el,d (s)EC ν ε el,d (s). ( 11
)
The evolution law is prescribed according to the standard format of non-associated irreversible processes:

Ḋ = φ ∂g ∂ ȳ , (12)
with the classical Kuhn-Tucker conditions φ ≥ 0, f ≤ 0 and φf = 0. The dot over a variable indicates differentiation with respect to time, g is the evolution potential controlling the growth of damage and φ is the damage multiplier. In this paper for the sake of simplicity, we assume g = ȳ.

The rate form of ( 8) reads

σ = (1 -D)EC ν εel,d - Ḋ 1 -D ε el,d . ( 13
)

Rate-type creep-damage law

The proposed rate-type formulation of a viscoelastic damageable non-ageing continuum assumes a deteriorating Young's modulus E. In fact, the effective stress σ t is defined first according to:

σ t = σ 1 -D (14)
Then, it is assumed that the relationship between the effective stress and the total strain follows a linear viscoelastic model. This is a rather classical approach, inspired from Lemaitre and Chaboche [START_REF] Lemaitre | Mécanique des matériaux solides[END_REF], which turns out to induce a decrease of the Young's modulus in the constitutive relations in our case (without affecting the characteristic times in the Kelvin chain). Consequently, we obtain the following stress-strain relation

σ = (1 -D)EC ν εtot -ḊEC ν ε el,d + (1 -D)EC ν N i=1 ε cr i τ i -(1 -D)E N i=1 σ τ i E i . (15)
The above model in tensorial notation reads

σ = (1 -D)C : εtot -ḊC : ε el,d + (1 -D)C : N i=1 ε cr i τ i -(1 -D)E N i=1 σ τ i E i , (16)
where the colon denotes the contracted tensorial product, C is the elastic material stiffness tensor, and by X we denote the tensorial transcription of the engineering matrix X.

The present model is also based on the assumption that only a fraction of the energy release rate due to creep contributes to the damage evolution with time. The motivation for introducing this parameter into the constitutive law is that for low stress levels although creep strain can be large (even larger than that corresponding to peak stress for short term loading) there is no significant variation of the elastic modulus; that is no significant damage [START_REF] Mazzotti | Experimental study of nonlinear creep of concrete at high stress levels[END_REF]. Hence, we assume that the function y(s) defined in ( 11) is now equal to

y(s) = 1 2 σ t (s) : ε el,d (s) + ασ t (s) : N i=1 ε cr i (s) , (17)
or equivalently:

y(s) = 1 2 ε el,d (s) : C : ε el,d (s) + αε el,d (s) : C : N i=1 ε cr i (s) , (18)
where the creep strain is

N i=1 ε cr i = ε tot -ε el,d , (19)
and the coefficient α is such that 0 ≤ α ≤ 1.

Consider now an initial state of equilibrium at time t 0 denoted by the state variables ε el,d 0 and D 0 . The rate constitutive relation describing the behavior of the material from this state is

σ = (1-D 0 )C : εtot -ḊC : ε el,d 0 +(1-D 0 )C : N i=1 ε cr i 0 τ i -(1-D 0 )E N i=1 σ 0 τ i E i . ( 20
)
3 Strain localization analysis

Equations of motion

The equations of motion are a set of non-linear integro-differential equations since the constitutive relations themselves are non-linear. Upon linearization of the equation of motion about the initial state ε tot 0 , ε el,d 0 , D 0 , the momentum

equations become div σ(x) = ρ ∂ 2 v ∂t 2 , ( 21
)
where v is the time derivative of the perturbation applied to the initial state.

Equations ( 21) are still non-linear due to the constitutive relations. Linearization is now performed under the assumption ḟ = 0. This assumption is classical in the analyses of localization. The solid that follows such a constitutive relation is called linear comparison solid [START_REF] Hill | Some basic principles in the mechanics of solids without natural time[END_REF]. Under this assumption and using the damage law (12) and equation ( 18), the equations of motion become div

  (1 -D 0 ) C : εtot -C : ε el,d 0 F (ȳ 0 ) V 1 2 ψ(s) 2ε el,d 0 (x + s) : C : εel,d (x + s) + α ε el,d 0 (x + s) : C : M i=0 εcr i (x + s) + M i=0 ε cr i 0 (x + s) : C : εel,d (x + s) d(s) + (1 -D 0 )C : N i=1 ε cr i 0 τ i -(1 -D 0 )E N i=1 σ 0 τ i E i   = ρ ∂ 2 v ∂t 2 . ( 22
)
The initial state of deformation and damage is assumed homogeneous through out the solid of volume V . The volume of the solid is assumed to be large enough so that boundary layer effects introduced by spatial averaging can be neglected.

Let us now consider the propagation of a harmonic wave in the direction

defined by n v = Ae -iξ(n.x-ct) , ( 23
)
where ξ is the wave number, c is the phase velocity, A is the amplitude of the perturbation and i is the imaginary constant such that i 2 = -1. The resulting rate of deformation is

εtot = - 1 2 iξ (A ⊗ n + n ⊗ A) e -iξ(n.x-ct) , (24)
where ⊗ denotes the tensorial product.

This harmonic perturbation is admissible if it satisfies the rate equation of equilibrium. Under the assumption that the state variables ε el,d 0 and D 0 are constant throughout the solid, taking into account (19) and noticing that εcr i defined by

εcr i = σ 0 τ i E i - ε cr i 0 τ i (25)
depends on the initial state only and is homogeneous over the solid, the following relation is obtained:

n.(1 -D 0 )C.n -ψ(ξn)F (ȳ 0 ) n. 1 2 2-α C : ε el,d 0 ⊗ε el,d 0 : C.n+n. 1 2 αC : ε el,d 0 ⊗ε tot 0 : C.n .A = ρc 2 A, (26)
where ψ(ξn) is defined by

ψ(ξn) = V ψ(s)e -iξn.s d(s). (27)
In order to obtain (26), one must substitute ( 23), ( 24) and ( 25) into (22).

Since the solid is assumed large, ψ(ξn) reduces to the Fourier transform of the weighting function. By considering the weighting function as isotropic, ψ(ξn) becomes independent of the direction n, hence ψ(ξn) ≡ ψ(ξ). Equation ( 26) can be written equivalently as an eigenvalue problem

[n.H * .n -ρc 2 I].A = 0, ( 28
)
where I is the second order identity tensor and

n.H * .n = n.(1 -D 0 )C.n -ψ(ξn)F (ȳ 0 ) n. 1 2 2 -α C : ε el,d 0 ⊗ ε el,d 0 : C.n + n. 1 2 αC : ε el,d 0 ⊗ ε tot 0 : C.n . ( 29
)
Equation ( 26 , 1993). Now, we restrict our analysis for the sake of simplicity to the uniaxial case or to loading paths such that

ε el,d 0 = κ 0 ε tot 0 , (32)
where κ 0 is a real such that 0 ≤ κ 0 ≤ 1 (Appendix A). Using (32), we obtain after some algebra

1 2 (2 -α)κ 2 0 + ακ 0 n.C : ε tot 0 . (n.C.n) -1 .ε tot 0 : C.n = (1 -D 0 ) ψ(ξn)F (ȳ 0 . ( 33
)
Let now,

ψ(x) = ψ 0 exp - x 2 2l 2 c , (34)
where l c is the internal length of the non-local continuum and ψ 0 is a normalizing factor. The Fourier transform of the above is

ψ(ξn) = exp - ξ 2 l 2 c 2 . ( 35
)
Solving (31) consists in finding the normal n and the wavelength 2π/ξ satisfying this equation. A geometrical approach proposed by Benallal [START_REF]Localization phenomena in thermo-elastoplasticity[END_REF] for the analysis of localization phenomena in thermo-elasto-plasticity is used here. Under the assumption of constant Poisson ratio

C ≡ λδ kl δ mn + µ (δ km δ ln + δ kn δ lm ) , (36)
where

λ = νE (1 + ν)(1 -2ν) , µ = E 2(1 + ν) , (37)
are the Lamé constants. One then gets

(n.C.n) -1 = I µ - λ + µ µ (λ + 2µ) n ⊗ n, (38)
C : ε tot 0 .n = λ tr(ε tot 0 )n + 2µε tot 0 .n, (39a)
n.C :

ε tot 0 .n = λ tr(ε tot 0 )n + 2µn.ε tot 0 .n. ( 39b
)
Substitution of ( 38), ( 39) into (33) yields after some algebra

4µ(E tot ) 2 + 4µ 2 λ + 2µ e tot + λ 2µ tr(ε tot 0 ) 2 = 1 -D 0 1 2 (2 -α)κ 2 0 + ακ 0 ψ(ξ)F (ȳ 0 ) , (40)
where we have set

e tot = n.ε tot 0 .n, (41a)
(E tot ) 2 = ε tot 0 .n . ε tot 0 .n -n.ε tot 0 .n 2 . ( 41b
)
E and e are respectively the tangent and normal components of the total strain vector in the direction n. Therefore, the initial state of deformation ε n is unique and it is obtained after ( 40)

l(n) = 2π ξ = πl c √ 2 / log 4µ(E tot ) 2 + 4µ 2 λ + 2µ e tot + λ 2µ tr(ε tot 0 ) 2 1 2 (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) (1 -D 0 ) . ( 42
)
This wavelength is related to the width of the localization zone.

[Fig. 2 about here.]

The above relation can be written also in terms of the elasto-damage strain

l(n) = 2π ξ = πl c √ 2 / log 4µ(E el,d ) 2 + 4µ 2 λ + 2µ e el,d + λ 2µ tr(ε el,d 0 ) 2 1 2 2 -α + α κ 0 F (ȳ 0 ) (1 -D 0 ) . (43)
[Fig. 3

about here.]

Creep effect on the width of the localization zone

In order to study the effect of creep on the width of the localization zone, we consider two different load histories. The state variables will be denoted with the index sc for the first one, and with index mc for the second one, in which the ratio κ 0 of elasto-damage to total strain is smaller κ sc 0 > κ mc 0 .

We assume that at the bifurcation points of the two different histories, the elasto-damage strain is the same ε

(el,d),mc 0 = ε (el,d),sc 0
, an assumption that yields

D mc 0 > D sc 0 .
Moreover, since the initial elasto-damage strain is the same for both loading paths one may prove that the normal n will be the same as well (Appendix B). Considering α constant and from (43), the wavelengths of the localized mode are such that l mc (n) < l sc (n) as long as

F (ȳ mc 0 ) 1 -D mc 0 > F (ȳ sc 0 ) 1 -D sc 0 , (44)
or equivalently, in account of ( 9), as long as

d log F dȳ ≤ -1. (45)
Hence, provided the above conditions are fulfilled, the larger the creep strain at constant total strain, the smaller the wavelength of the localized mode. This decrease of wavelength provides also a decrease of the fracture process zone, an increase of brittleness and a shift on the size effect plot for geometrically similar specimens that is consistent with the experimental results in Fig. (B.1).

Role of α on the width of the localization zone

The role of α, defined in Eqn. ( 18) is quite straightforward. If we replace α in ( 43) by α such that α ≤ α, then we conclude that the width of the localization zone increases. In other words, if the contribution of the creep strain in the growth of damage decreases, the wavelength of the localized mode increases. This produces an increase of the size of the fracture process zone and an increase of ductility. For geometrically similar specimens, this increase of ductility is materialized by a shift of the experimental data point towards a strength criterion, to the left.

When damage depends solely upon the elastic energy (α = 0), the bifurcation criterion becomes a function of only the elastic strain and creep effects can no longer be reproduced.

Numerical example

Consider the uniaxial initial state of strain

ε tot 0 =        tot 01 0 0 0 -ν tot 01 0 0 0 -ν tot 01       
.

The critical wavelengths can be obtained when the ellipse corresponding to the above initial state of deformation is tangent to the largest Mohr circle leading to the localization criterion

E(1 + ν) 2 tot 01 sin(2θ) 2 + E(1 -2ν)(1 + ν) 4(1 -ν) 1 + cos(2θ) 2 ( tot 01 ) 2 - (1 -D 0 ) 1 2 (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) ψ(ξ) = 0, ( 46
)
where the angle θ defines the normal n = (cos θ, 0, sin θ) to the localization band and satisfies

tan 2 (2θ) = 2 E(1-ν)(1+ν) 1-2ν (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) ψ(ξ)( tot 01 ) 2 -(1 -D 0 ) (1 -D 0 ) -E(1 + ν) (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) ψ(ξ)( tot 01 ) 2 , (47)
as it can be proved with Mohr analysis.

[Fig. 4 about here.]

We particularize function F , for illustration purposes, as

F (ȳ) = b 1 + 2b 2 (ȳ -y 0 ) 1 + b 1 (ȳ -y 0 ) + b 2 (ȳ -y 0 ) 2 2 , (48)
where the numerical values of the model parameters are E = 32.000 M P a ,

ν = 0.2, b 1 = 605 M P a -1 , b 2 = 5.
42 10 4 M P a -2 and y 0 = 6 10 -5 M P a.

Damage is then given by

D(ȳ) = 1 - 1 1 + b 1 (ȳ -y 0 ) + b 2 (ȳ -y 0 ) 2 . ( 49
)
As we have an infinite body ȳ = y = 1 2 1 -α + α/κ 0 E tot Solving (47) for ψ(ξ) and substituting in ( 46) yields

E(1 -2ν)(1 + ν) 4(1 -ν) 1 + cos(2θ) 2 -4 E(1 -ν)(1 + ν) 1 -2ν cos 2 (2θ) + 5E(1 + ν) 2 sin 2 (2θ) = 0. (50)
Thus, the direction of the normal n to the localization band does not depend, in the uniaxial case, upon the state of strain or on the damage or on the yield function of damage but only upon the elastic parameters of the material. Solving (50), one gets θ = π/4. This angle corresponds to a direction of localization in the (x 1 , x 2 ) or (x 1 , x 3 ) planes.

[Fig. 6 The proposed coupled model preserves localization limiting properties, same as in the rate independent case (Pijaudier-Cabot and [START_REF] Pijaudier-Cabot | Strain localization and bifurcation in a nonlocal continuum[END_REF] and it is consistent with the effect of creep observed experimentally. Namely, the minimum wavelength of the localization modes cannot be zero in accordance with energy considerations and the admissible wavelengths decrease when the creep effect (strains) increases, at least under a certain condition on the growth of the yield function of damage. This last property should induce a variation on size effect data for geometrically similar specimens: a shift in the size effect plot towards increasing brittleness for an increasing creep strain. Moreover, when the fraction of creep strain which contributes to damage is decreasing, a shift on the size effect plot occurs towards the left hand side, i.e. towards the strength criterion and an increase of ductility is expected. The model fails to predict the influence of creep on the wavelength of localized modes when damage is considered to depend solely upon the elastic energy.

Such analytical results are expected to be recovered in computational failure analyses. They are expected to have important practical consequences:

as creep develops, the residual capacity of a given structure may decrease;

the structure is becoming more brittle, safety margins may be reduced. This of course should hold for severely loaded concrete structures, where damage interacts with creep. It may typically occur in some regions of prestressed concrete structures subjected to stress concentrations.

A

Consider a time-independent process until the column matrix stress takes the value σ ind and then the loading state remains fixed for the time interval [t ind , t 0 ], where t 0 is the time at which bifurcation occurs. For this loading path relations (4) would give

ε cr 0 = N i=1 ε cr i 0 = N i=1 1 -e -t/τ i E i C -1 ν σ ind . (A.1) Since ε el,d 0 = 1 (1 -D 0 )E C -1 ν σ ind , (A.2)
one may conclude that (32) holds true for

κ 0 = 1/ (1 -D 0 )E 1/ (1 -D 0 )E + N i=1 1 -e -t/τ i /E i . (A.3)
In the same spirit, relation (32) can be proved for more general loading paths for which the loading state is either fixed or changes through a time-independent process.

In the uniaxial loading now, at the onset of bifurcation, relation (6) reads

   tot 01 -ν tot 01 -ν tot 01 0 0 0    T = aC -1 ν    σ 01 0 0 0 0 0    T , (A.4) that is ε tot 0 = C -1 ν aσ 0 , (A.5)
where for simplicity we assumed t = 0 and a is a number such that

aσ 0 = t 0 0 J(t, τ )dσ(τ ).
Moreover,

ε el,d 0 = 1 (1 -D 0 )E C -1 ν σ 0 , (A.6)
holds true. Using now (A.5) and (A.6), we deduce (32) with κ 0 = 1/ (1 -D 0 )Ea .

B

The normal n = (cos θ, 0, sin θ) to the localization band can be obtained when the ellipse corresponding to the initial state of strain ε tot 0 is tangent to the largest Mohr circle leading to system It is then obvious from (B.4) that, as long as (32) holds true, the normal n remains fixed for fixed initial elasto-damage strain ε el,d 0 .

tan 2 (2θ) = 2(λ + 2µ) (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) ψ(ξ)( tot 01 -tot 03 ) 2 -(1 -D 0 ) (1 -D 0 ) -2µ (2 -α)κ 2 0 + ακ 0 F (ȳ 0 ) ψ(ξ)( tot 01 -tot 03 ) 2 , ( B

  , upon a decrease of the loading rate. It means that the set of data points are getting closer to the LEFM criterion (Figure (B.1)).

  , as shown in Figure (B.1). On the other hand, if the rate dependence is caused only by bond breakages, the peak loads correspond to increasing displacements and there is no shift of this kind on the size effect plot. Note that ductility and brittleness are defined here, and in the remaining, as the position of the experimental data with respect to the two asymptotic criteria: strength of materials and LEFM (see Figure (B.1)). Ductility increases when the set of data points shift to the left and brittleness increases when they shift to the right. [Fig. 1 about here.]

  ) admits non-trivial solutions if and only if det[n.H * .n -ρc 2 I] the localization condition in a non-local continuum derived by Pijaudier-Cabot and Benallal (Pijaudier-Cabot and Benallal

  tot 0 maps the (e, E) plane into the Mohr circles of deformation for all the directions n and the bifurcation criterion maps into a set of ellipses in this plane. The size of the ellipses increases as the wave number ξ increases. Before localization the smallest possible ellipse contains the largest Mohr circle. Bifurcation will occur for the first time when the smallest possible ellipse is tangent to the largest Mohr circle of deformation (the smallest possible ellipse corresponds to ψ(ξ) = 1 or equivalently to ξ = 0). The normal vector n can be calculated geometrically from standard Mohr analysis (see Figure (B.2). After the first occurrence of bifurcation the ellipse will intersect the largest Mohr circle. In this case there is a finite set of normal vectors n such that bifurcation is possible. However, each vector n corresponds to the same ellipse whose radius is defined by ξ (Figure (B.3)). The wavelength corresponding to each vector

F

  satisfies criterion (45) for y ≤ y cr = 2.995. y cr yields the damage value D cr 0.999998 (Figure (B.4)). [Fig. 5 about here.]

  Figure (B.5) shows the admissible wavelengths normalized to the character-

Fig. B. 1 .Fig

 1 Fig. B.1. Nominal strengths of 4 groups of 3 specimens of different sizes tested at 4 different times to peak (after[START_REF] Bažant | Rate effects and load relaxation: static fracture of concrete[END_REF]

This work was partially supported from the EU project 'Degradation and Instabilities in Geomaterials with Application to Hazard Mitigation' (DIGA-HPRN-CT-2002-00220) in the framework of the Human Potential Program, Research Training Networks.