Convergence of a finite volume scheme for a parabolic system with a free boundary modeling concrete carbonation

In this paper we deﬁne and study a ﬁnite volume scheme for a concrete carbonation model proposed by Aiki and Muntean in [1]. The model consists in a system of two weakly coupled parabolic equations in a varying domain whose length is governed by an ordinary diﬀerential equation. The numerical sheme is obtained by a Euler discretisation in time and a Scharfetter-Gummel discretisation in space. We establish the convergence of the scheme. As a by-product, we obtain existence of a solution to the model. Finally, some numerical experiments show the eﬃciency of the scheme.


Introduction
The carbonation phenomenon in reinforced concrete (a material widely used in civil engineering for the construction of buildings, factories, bridges, roads, etc) is a physico-chemical reaction which is the main cause of concrete structure degradation. There exists a wide literature on the modeling of this reaction, see [15,17,18,19] and all references therein. From a physical point of view the carbonation process can be explained as follows: CO 2 from the atmosphere enters in the concrete via the unsaturated porous media matrix where it is quickly transformed in CO 2 in the aqueous phase. This reaction can be described by CO 2 (g) ←→ CO 2 (aq).
This reaction facilitates a drop of the pH inside the material and allows the corrosion process to damage the metallic reinforcement bars. These damages deteriorate the concrete and reduce the durability of the structure. From a mathematical point of view, one of the main source of interest for this problem is the existence of a moving domain (like for instance in the corrosion theory, see [5]). Indeed, the carbonation process produces a moving interface (front position in one dimension) which splits the concrete in two regions: the carbonated one, which grows in time, and the uncarbonated one. In this paper we consider a free-boundary system proposed by Aiki and Muntean in [1] where the varying space domain represents the carbonated zone. The unknowns u and v represent the mass concentration of CO 2 respectively in aqueous (u) and gaseous (v) phase and s represents the penetration depth which measures the size of the carbonated zone. We denote the carbonated domain Q s (T ) defined by Q s (T ) = {(t, y) : 0 < t < T < +∞, 0 < y < s(t)} .
Then, the system considered by Aiki and Muntean writes: s (t) = ψ(u(s(t), t)) for 0 < t < T, u(0, t) = g(t) for 0 < t < T, v(0, t) = r(t) for 0 < t < T, −κ u ∂ y u(s(t), t) − s (t)u(s(t), t) = ψ(u(s(t), t)) for 0 < t < T, u(y, 0) = u 0 (y) for 0 < y < s(0), v(y, 0) = v 0 (y) for 0 < y < s(0), Existence and uniqueness of a global solution to (1) are established in [1]. As in the theoretical analysis, we consider that the following assumptions are satisfied: (A1) ψ : R −→ R represents the kinetics of the carbonation reaction and is defined by ψ(x) = αx with α > 0, (A2) f : R 2 −→ R is given by the Henry's law and is defined by f (u, v) = β(γv − u) with β and γ two real constants, (A3) g and r belong to W 1,2 (0, T ), (A4) u 0 and v 0 belong to L ∞ ([0, s 0 ]), (A5) the diffusive coefficients κ u and κ v are two positive constants, (A6) s 0 > 0, (A7) there exist two positive constants g * and r * with g * = γr * such that 0 ≤ g ≤ g * , 0 ≤ r ≤ r * on [0, +∞), Remark 0.1. In [1], the authors consider, for the kinetics of the reaction, ψ(x) = α(x + ) p where x + = max(x, 0) and with p ≥ 1. As we will see later, the positive part is not necessary in the definition of ψ due to the nonnegativity of u. We note that it is also possible to consider a nonlinear source term f , corresponding to a nonlinear Henry's law as in [4]. However, in this paper we focus on the numerical analysis of a finite volume scheme for (1) in the linear case: p = 1 and with a linear f . The main difficulty is due to the coupling between the reaction-diffusion equations and the ordinary differential equation governing the evolution of the domain. The techniques we develop could be adapted to the nonlinear case.
There exist few convergence results for numerical schemes dedicated to carbonation models. In [16], Muntean develop a semi-discrete finite element scheme in space for a more complex model than (1) written on a varying space domain. He proves that the semi-discrete scheme converges with an order one. In [21], the authors show the convergence of a numerical scheme obtained by a semiimplicit discretisation in time and a mixed finite element method in space for a concrete carbonation model on a fixed space domain. We also refer to [20] where the authors prove the convergence of a Galerkin method for a class of multiscale reaction-diffusion systems, still defined on a fixed space domain. Finally in [23], the author defines a finite volume scheme which approximate the solutions of a nonlinear parabolic system modeling the carbonation process in a porous concrete material occupying a fixed space domain.
In this paper, we propose a numerical scheme obtained by a Euler discretisation in time and a finite volume discretisation in space. Our purpose is to study the convergence of the scheme towards the weak solution of (1). As it is usual in the finite volume framework, the convergence proof is based on some energy estimates, see Propostion 3.1 and Proposition 3.2, which yield the compactness of the sequence of approximate solutions. However, even if the method is classical, the originality of the proof is due to the coupling between the two partial differential equations (1a) and (1b) and the ordinary differential equation (1c) which governs the growth of the domain. The paper is organized as follows. In Section 1, we first rewrite the system (1) on a fixed domain and we introduce the numerical scheme. Then, we state the main results: Theorem 1.1 gives the existence and the uniqueness of the solution to the scheme, while Theorem 1.2 gives the convergence result. Section 2 is devoted to the proof of Theorem 1.1. Then, in Section 3, we establish some energy inequalities satisfied by the approximate solutions. The proof of Theorem 1.2 is given in Section 4. Finally, in Section 5 we discuss some numerical experiments. We show some profiles and illustrate the convergence of the scheme. Moreover, we investigate the long time behavior of the scheme: the penetration depth follows a √ t-law of propagation, which fits experimental observations and which has been established in [2,3].
1 Numerical scheme and main results

The drift-diffusion model on a fixed domain
For numerical reasons, it is convenient to rewrite (1) on a fixed space domain. To this end, we follow the ideas of [5,16]: we consider the following change of variables , t .

Let us consider
The diffusion-reaction equation Multiplying (2) by s 2 (t), we obtain the following convection-diffusion-reaction equation We use the same technique for the equation (1b). If we drop the bars, the system (1) rewrites as The general convection-diffusion fluxes write where w refers to either u or v. We also use this notation w = u or v without further mention in the sequel. Let us now define the notion of weak solution of (3). We consider the functional space H = {z ∈ H 1 (0, 1) : z(0) = 0}, endowed with the H 1 (0, 1) norm. Assuming (A1) − (A7), we say that (s, u, v) is a weak solution of (3) if the following conditions (S1)-(S5) are satisfied: As already mentioned the existence and uniqueness of a weak solution to (3) has been proved in [1]. But, in this paper, Aiki and Muntean considered a different functional space for (S1). Indeed, they add the condition that ∂ t u and ∂ t v ∈ L 2 (0, T ; H * ) with H * the dual space to H. Nevertheless, with the assumptions (A1)−(A7) and the conditions (S1)-(S5) we can prove that ∂ t u and ∂ t v are actually in L 2 (0, T ; H * ). As a consequence, we deduce the equivalence of the two definitions of a weak solution to (3).

Definition of the numerical scheme
In order to write the finite volume scheme we introduce some notation related to the discretization of [0, 1]. A mesh T , consists in a finite sequence of cells denoted by ( , for 1 ≤ i ≤ l, the length of the i-th cell. The mesh size is defined as h = max{h i , 1 ≤ i ≤ l}. Moreover, for 1 ≤ i ≤ l, we define x i as the center of the cell (x i− 1 2 , x i+ 1 2 ), x 0 = x 1 2 and x l+1 = x l+ 1 2 . We set We define a time step ∆t and we assume that there exists an integer N T such that N T ∆t = T (we can define N T as the integer part of T /N T ). We consider the sequence (t n ) 0≤n≤N T with t n = n∆t. Then, for 1 ≤ i ≤ l and 0 ≤ n ≤ N T − 1, the scheme writes Let us highlight that the time discretisation ensures a decoupling between (4a), (4b) and (4c): knowing (s n , v n , u n ), s n+1 , v n+1 and u n+1 are computed sucessively in an explicit manner. The term G n+1 w,i+ 1 2 is the numerical approximation of J w (x i+ 1 2 , t n+1 ). We choose to discretize simultaneously the diffusive part and the convective part of the fluxes J w . More precisely, we consider the Scharfetter-Gummel fluxes intoduced by Il'in in [13] and Scharfetter and Gummel in [22]. We set where B is the Bernoulli function defined by This choice of numerical fluxes is motivated by the fact that the Scharfetter-Gummel scheme exhibits a second order convergence rate in space, see [10,14]. We supplement the numerical scheme with the discretization of the boundary and of the initial conditions We denote by (S) the scheme (4)-(7).

Main results and outline of the paper
We first state the existence and uniqueness result of a solution to (S). Moreover, as in the continuous case, we prove that the solutions of (S) fulfill lower and upper bounds.
Then, there exists a unique solution to the numerical scheme (S). Furthermore, under the assumptions (A1)-(A7), we have for every n ≥ 0 and i ∈ {0, · · · , l + 1} and For 0 ≤ n ≤ N T , we deduce Remark 1.1. We notice that the inequality (10) is far away from the √ tbehavior mentioned in the introduction. Nevertheless, the estimate (10) suffices fir the convergence analysis of the scheme.
As it is usual for finite volume method (see [11]), we define some piecewise constant functions in space and time. For a given mesh T and a given time step ∆t, we define For these functions we define a discrete derivative operator in space as We can also reconstruct a piecewise affine function s ∆t , for a given time step ∆t, we set The choice of this piecewise affine reconstruction for s m will allow us to apply Ascoli's compactness theorem in the proof of the convergence of the scheme (see Theorem 1.2). Then, we consider a sequence of meshes and time steps (T m , ∆t m ) m such that h m ↓ 0 and ∆t m ↓ 0 as m ↑ ∞. As a consequence of Theorem 1.1, we can define a sequence of solutions to (S), denoted (s m , u m , v m ) m with s m = s ∆tm and w m = w hm,∆tm . We establish the convergence of the sequence (s m , u m , v m ) m , as stated in the following theorem.
Theorem 1.2. Let us assume that the hypotheses (A1)-(A7) are satisfied. Then, the sequence (s m , u m , v m ) m converges to some (s, u, v) with s ∈ W 1,∞ (0, T ) and and (s, u, v) is the weak solution to (3).
We prove Theorem 1.1 in Section 2. The proof is based on some monotonicity properties of the numerical scheme. In Section 3, we establish L 2 (0, T ; H 1 (0, 1)) and H 1 (0, T ; H * ) discrete estimates for the approximate solutions. These estimates allow us to apply a discrete version of Aubin-Simon lemma, see [12], and we deduce some compactness properties of the sequence (s m , u m , v m ) m . In Section 4, we pass to the limit in the scheme and we prove Theorem 1.2. The proof is based on some arguments developed in [6,7]. Eventually, in Section 5, we present some numerical experiments. We draw profiles of u, v and s. Moreover, we investigate the question of the L 2 -convergence rate in space of the numerical scheme.
2 Proof of Theorem 1.1 In this Section, we prove Theorem 1.1. To this end, we proceed by induction on n and we follow some ideas developed in [5]. First, let us note that (s 0 , u 0 , v 0 ) is well defined and satisfies (8)- (10). We now assume that (8)-(10) is verified for (s n , u n , v n ), for some n ≥ 0. The inequalities (9) and (10) are straightforward consequences of the definition of s n+1 . Furthermore, as a by-product of the lower bound of (9), we observe that C n+1 ≥ 0. We first notice, thanks to (6b) and the definition of B, that we can express v n+1 l+1 as As already mentioned, the equations (4) are decoupled. We rewrite (4b) and (4c) as two independent linear systems where we consider that u n+1 = (u n+1 The matrices M n u and M n v are tridiagonal and defined by The vectors b n u and b n v are defined by The matrices M n u and M n v have positive diagonal terms and non-positive offdiagonal terms and are strictly diagonally dominant with respect to their columns. Therefore, M n u and M n v are M-matrices and thus invertible and monotone. We first deduce the existence and uniqueness of a solution to (S). Moreover, as b n u ≥ 0 and b n v ≥ 0, we deduce, thanks to the induction hypothesis, that u n+1 ≥ 0, v n+1 ≥ 0 and by definition of v n+1 l+1 , we conclude that v n+1 l+1 ≥ 0. Let us prove now the upper bounds (8). We notice that the inequalities are true for i = 0 by construction of the scheme (S). Now let R * ∈ M l,1 (R) be the vector with all its components equal to r * . Using the equality B(x) − B(−x) = −x, we obtain that, for 2 ≤ i ≤ l − 1, Thanks to the induction hypothesis and to the identity g * = γ r * , it yields for Then, for i = l, we have Thus, using the induction hypothesis, we obtain and v n+1 l+1 ≤ r * follows from the definition of v n+1 l+1 . Similarly, let G * ∈ M l+1,1 (R) a constant vector of component g * , we show that for i ∈ {1, · · · , l} and (M n u u n+1 ) l+1 = 0. Eventually, (M n u (u n+1 − G * )) l+1 ≤ 0 and as M n u is an M-matrix we deduce the upper bounds of (8). This concludes the proof of Theorem 1.1.

Functional spaces and discrete norms
For the proof of Theorem 1.2, we introduce some discrete functional spaces and norms.
Definition 3.1. Let T a mesh of size h and ∆t a time step. We define the set X T of piecewise constant functions in space by We also define the set X T ,∆t of piecewise constant functions in space and time by For z h ∈ X T , we notice that if || · || 0 denote the L 2 (0, 1) norm, then On X T and X T ,∆t we define some norms which are the discrete version of the H 1 (0, 1), H * , L 2 (0, T ; L 2 (0, 1)), L 2 (0, T ; H 1 (0, 1)) and L 2 (0, T ; H * ) norms. , ∀z h,∆t ∈ X T ,∆t .
In [8] the authors have shown that, for all z h ∈ X T , As a by-product, we obtain the discrete Poincaré inequality, for all z h ∈ X T 3.2 Discrete L 2 (0, T ; H 1 (0, 1)) estimates In this section we prove L 2 (0, T ; H 1 (0, 1)) discrete estimates for the solutions of (S). Let ∆t be a time step. Then, we introduce a piecewise constant reconstruction of the boundary condition where, for 0 ≤ n ≤ N T − 1, the elements g n+1 and r n+1 are defined by (6a). For these functions we define a discrete seminorm, denoted || · || 1,2,∆t , defined by for q ∆t = g ∆t or r ∆t . As g and r ∈ W 1,2 (0, T ), there exist G, R < +∞, not depending on ∆t, such that ||g ∆t || 1,2,∆t < G and ||r ∆t || 1,2,∆t < R.
Proof. Let us show the result for v h,∆t . We multiply (4b) by ∆t s n+1 s n (v n+1 i − r n+1 ) and we sum over n and i, we obtain E + F + G = 0, with We notice that we can rewrite E as E = E 1 + E 2 + E 3 , with For E 1 , we apply the formula (a − b) a ≥ (a 2 − b 2 )/2, to get a telescopic sum. Then, we obtain For E 2 , using the inequality 2ab ≥ −a 2 − b 2 , we deduce Finally, applying the Cauchy-Schwarz inequality, we have for E 3 Thus, from Theorem 1.1 and the Young inequality, we deduce that Hence, we obtain For the term F , we use a discrete integration by parts Following the ideas of [6], we recall the standard decomposition of the numerical fluxes: with Then, applying (16) to F , we have Using B c (x) ≥ 1 for all x ∈ R and reordering the terms, we obtain
Let us prove now the result for u h,∆t . As previously, we multiply (4c) by ∆t s n+1 s n (u n+1 i − g n+1 ) and we sum over n and i, so that we obtain similarly E + F + G = 0. The only difference with the previous case concerns the term F . Indeed, due to the boundary conditions, the discrete integration by parts yields to Thus, we define Using the boundary condition (6c) and Theorem 1.1, we deduce that Then, applying the same techniques as before we conclude that there exists a positive constant C only depending on s 0 , α, ||g ∆t || 1,2,∆t , g * , r * , κ u , β, γ and T such that ||u h,∆t || 0;1,2,T ≤ C.

Discrete H 1 (0, T ; H * ) estimates
In order to show some compactness properties for the sequence of approximate solutions (s m , u m , v m ), we want to apply a discrete version of Aubin-Simon lemma obtained in [12]. In this purpose, we need to establish discrete L 2 (0, T ; H * )-estimates for ∂ t,∆t u m and ∂ t,∆t v m , where ∂ t,∆t denotes the discrete derivative operator in time defined, for all z h,∆t ∈ X T ,∆t , by Proposition 3.2. For a given ∆t, a given mesh T and under the assumptions (A1)-(A7), there exists a positive constant C depending only on s 0 , α, g * , r * , κ w , β, γ and T and independent of h and ∆t such that Proof. As in the previous proof, we first show the result for v h,∆t . Let η h ∈ X T such that η h (0) = 0 and ||η h || 1,2,T ≤ 1. We first notice that For a given n we multiply (4b) by 1 s n+1 s n η i and we sum over i, we obtain Using the Cauchy-Schwarz inequality, Theorem 1.1 and (14), we obtain For G, we use a discrete integration by parts and apply (16), so that G = G 1 + G 2 , with Thanks to Theorem 1.1, the term C n+1 is bounded. Then, we observe that B c (h i+ 1 2 C n+1 x i+ 1 2 /κ v ) is also bounded for 0 ≤ i ≤ l. We conclude, from Theorem 1.1 and the Cauchy-Schwarz inequality that Applying Theorem 1.1, we obtain that Thus, the Cauchy-Schwarz inequality yields to Hence, we conclude that Using the Cauchy-Schwarz inequality, we have Eventually, the L ∞ -estimates and (14) provide Then, from (19), (20) and (21), we deduce: There exists a constant C depending on s 0 , α, g * , r * , κ v , β, γ and T such that Finally, we multiply by ∆t, we sum over n and we get Proposition 3.1 provides the existence of a positive constant C independent of h and ∆t such that Let us prove now the result for u h,∆t . Let η h ∈ X T such that η h (0) = 0 and ||η h || 1,2,T ≤ 1. For a given n we multiply (4c) by 1 s n+1 s n η i and we sum over i, we obtain E = F +G+H. As in the proof of Proposition 3.1, the only difference concerns the term G. Then, after a discrete integration by parts to G, we have Thanks to (6c) and Theorem 1.1, we deduce that Then, applying the same techniques as before we conclude that there exists a positive constant C depending only on s 0 , α, g * , r * , κ u , β, γ and T and independent of h and ∆t such that 4 Proof of Theorem 1.2

Compactness results
We now consider a sequence of meshes and time steps (T m , ∆t m ) m , such that h m ↓ 0 and ∆t m ↓ 0 as m ↑ ∞, and the associated sequence (X Tm ) m . Let us notice that (X Tm ) m is a sequence of finite-dimensional subspaces of L 2 (0, 1). Furthermore, we can endow each X Tm with the || · || 1,2,Tm norm or with the ||·|| −1,2,Tm norm. These norms achieve the hypotheses of Gallouët-Latché lemma (lemma 3.1 in [12]). For a rigorous proof of this fact one can read Lemmas 2.1 and 2.2 in [9]. Eventually, for each m, we define (s m , u m , v m ) a sequence of solution to (S) with s m = s ∆tm and w m = w hm,∆tm .
For the sequence (s m ) m , Theorem 1.1 and Ascoli theorem give the existence of a function s ∈ C([0, T ]) such that, up to a subsequence, Finally, using (9) of Theorem 1.1, we obtain the existence of q ∈ L ∞ (0, T ) with, up to a subsequence, Moreover, in the sense of distribution, s = q.

Convergence of the traces
In the next section we pass to the limit h ↓ 0 and ∆t ↓ 0 in the scheme. Let us first consider the boundary terms. Let us define the trace δz h,∆t , of z h,∆t ∈ X T ,∆t , by

Passage to the limit
In order to conclude the proof of Theorem 1.2, it remains to show that (s, u, v), obtained in Proposition 4.1 is a weak solution to (3). In order to prove that (s, u, v) satisfies (S3), (S4) and (S5), we follow some ideas developed in [6,9]. For sake of simplicity, for a given mesh T and a given time step ∆t, we define: where, for 0 ≤ k ≤ N T , the term w k h is defined by (11). For a sequence of meshes and time steps (T m , ∆t m ) such that h m ↓ 0 and ∆t m ↓ 0 as m ↑ ∞, we notice thatw wherew m =w hm,∆tm ,š m =š ∆tm ands m =s ∆tm . Let φ ∈ D([0, T ) × (0, 1]), we define We set Then, thanks to Proposition 4.1 and Proposition 4.2, we deduce In order to prove (S4), it remains to prove that m → 0 as m ↑ ∞. For this, we multiply (4c) by ∆t m φ n i , where φ n i = φ(x i , t n ), and we sum over i and n. We get Applying the standard method used in [6,9], we must show that for j ∈ {1, 2, 3} we have In our case, we notice that the only term which differ from [6,9] is the term |A u,1 − A u,1 |. For this reason, we only prove that |A u,1 (m) − A u,1 (m)| → 0 as m ↑ ∞ in the sequel. For the other terms we refer to [6,9]. By definition of A u,1 (m), we have For A u,1 (m), using a discrete integration by parts we get Then, using φ N i = 0 and inserting the term (−s n φ n i + s n φ n i ) we obtain Therefore, thanks to the regularity of φ, the inequalities (8), (9) and (10) and the Cauchy-Schwarz inequality, we conclude that there exists a constant C only depending on s 0 , g * , α and T such that Then, we obtain that m → 0 as m ↑ ∞ and (S4) is satisfied. Using the same method we may show that (s, u, v) verifies (S5).
We deduce from this inequality that u−g ∈ L 2 (0, T ; H) and the same arguments show that v − r ∈ L 2 (0, T ; H).
As already mentioned, we can show that u and v belong also to H 1 (0, T ; H * ). Then, thanks to [1], we deduce the uniqueness of the weak solution (s, u, v). We also deduce that the whole sequence (s m , u m , v m ) m converges to the weak solution of (3).
0.1 0.5 1 6.5 7.5 Table 1: Definition of parameters used in the test case behavior of the penetration depth and also in the order of convergence of the scheme. We may compare the results given by the Scharfetter-Gummel scheme presented in the paper and the upwind scheme corresponding to the B function defined by B(x) = 1 + x − with x − = max(−x, 0) (see [10]). Figure 1 shows the different profiles of v and u as functions of x where x ∈ [0, s(t)] for t ∈ {20, 40, 60, 80, 100}. The results are obtained with 1000 cells and ∆t = 10 −2 . We observe profiles similar to those given in [3,15]. In Figure 2, we illustrate the behaviour of s in linear scale on the left and in logarithmic scale on the right, up to time T = 1000. We observe that the penetration depth s follows a √ t-law of propagation (see [3,15]). The results are obtained with 1000 cells and ∆t = 10 −2 . Let us mention that the upwind scheme gives the same results than the Scharfetter-Gummel scheme.
Since the exact solutions u and v of (3) are not explicitly known, we compute two reference solutions on a uniform mesh composed of 2560 cells and with ∆t = (1/2560) 2 in the Scharfetter-Gummel case and ∆t = 1/2560 in the upwind case, in order to investigate the question of the convergence rate in space for (S). We impose this condition on ∆t because the Euler discretisation in time exhibits a first order convergence rate while we expect a second order convergence rate in space for the Scharfetter-Gummel scheme and a first order for the upwind scheme. Then, we compute approximate solutions on a uniform mesh made of respectively 10, 20, 40, 80, 160, 320, 640 and 1280 cells. Finally, we compute the L 2 -norm of the difference between the approximate solution and the average of the reference solution at different final times: T = 0.1, 1 and 5. The results are shown on Figure 3 for the Scharfetter-Gummel scheme and on Figure 4 for the upwind scheme. As expected (see [14]), the Scharfetter-Gummel scheme converges with an order around 2, while the upwind scheme converges with an order around 1.

Conclusion
In this paper, we have proven the convergence of an efficient finite volume scheme for the carbonation model introduced by Aiki and Muntean in [1]. As a by-product, we obtain a new proof of existence for this free-boundary system. Moreover, this scheme gives profiles of u and v similar to those given in the literature [3,15] and the simulations support the idea that the approximate penetration depth s ∆t follows a √ t-law of propagation as showed in [2,3]. Nevertheless, a rigorous justification of this result needs further investigations.