Analysis of a fully discrete approximation to a moving-boundary problem describing rubber exposed to diffusants

We present a fully discrete scheme for the numerical approximation of a moving-boundary problem describing diffusants penetration into rubber. Our scheme utilizes the Galerkin finite element method for the space discretization combined with the backward Euler method for the time discretization. Besides dealing with the existence and uniqueness of solution to the fully discrete problem, we derive a \textit{a priori} error estimate for the mass concentration of the diffusants, and respectively, for the position of the moving boundary. Numerical illustrations verify the obtained theoretical order of convergence in physical parameter regimes.


Introduction
We study the fully discrete approximation of a one-dimensional moving-boundary problem describing the penetration of diffusants into rubber.The model presented here was proposed recently in [24], where the simulation output was compared to experimental data.In this framework, our interest is focused on the numerical analysis of the model.Relying on previous mathematical analysis work done for an adsorption model with moving swelling interfaces (see [17]), which shares the structure of the equations with our current moving-boundary model, we have provided in [25] an analysis of the control of the errors produced by a FEM semi-discretization of our model equations.In this paper, we turn our attention to estimating the errors produced by combining time and space discretizations.Our analysis of the fully discrete approximation to the model equations relies on our previous results [17,25] and should be seen as a natural continuation of the work.Browsing the existing literature, one can find a lot of information regarding the rigorous error analysis of semi-discrete approximation of free-and moving-boundary problems.However, much less seems to be known what concerns the analysis of fully discrete approximation schemes even for one dimensional formulations where the moving interface is in fact only a moving point (with a priori unknown location).The references [2,19] were particularly useful for our investigation.In [19] H. Y. Lee develops a fully discrete scheme for a Stefan problem with non-linear free boundary condition and investigates the order of convergence of the scheme.In ref. [2], the authors construct and analyze fully discrete methods for a free boundary problem arising in the polymer technology.By using the Galerkin finite element formulation in space and a backward Euler scheme in time, the authors were able to prove the a priori error estimate for the concentration of the solvent and for the position of moving boundary.At the technical level, we were very much inspired by the technique that has been used in [2] to get the a priori error bound.It is also worth mentioning that the main difference between the problem considered in these papers and our problem lies in the choice of the boundary conditions.Both of the cited papers impose at least an homogeneous Dirichlet boundary condition at one of the boundaries, while we impose flux boundary conditions at both boundaries that bring in boundary terms that need a careful handling.
The problem setting we are studying here is as follows: For a fixed given observation time T f ∈ (0, ∞), let the interval [0, T f ] be the time span of the physical processes we are considering.Let x ∈ [0, s(t)] and t ∈ [0, T f ] denote the space and respectively, the time variable.Let m(t, x) be the concentration of diffusant placed in position x at time t.The diffusants concentration m(t, x) acts in the region Q s (T f ) defined by Q s (T f ) := {(t, x)|t ∈ (0, T f ) and x ∈ (0, s(t))}.
The problem reads: Find m(t, x) together with the position of the moving boundary (interface) x = s(t) for t ∈ (0, T f ) such that the couple (m(t, x), s(t)) satisfies the following evolution problem: − D ∂m ∂x (t, s(t)) = s (t)m(t, s(t)) for t ∈ (0, T f ), s (t) = a 0 (m(t, s(t)) − σ(s(t)) for t ∈ (0, T f ), ( 4) s(0) = s 0 > 0 with 0 < s 0 < s(t) < L, where a 0 > 0 is a kinetic coefficient, β is a positive constant describing the capacity of the interface at x = 0, D > 0 is the effective diffusion constant, H > 0 is the Henry's constant.Additionally, σ is a real function, b is a given boundary concentration on [0, T ], s 0 > 0 is the initial position of the moving boundary, while m 0 represents the initial concentration of the diffusant.This model reminds of the work by Astarita and collaborators (compare [6] and follow-up papers) on free boundary problems posed in the context of polymeric materials; see [10,12,23,9] for classical older works on the topic.We refer the reader to [7] to a very recent collection of modeling, mathematical analysis, and numerical simulation aspects of moving-boundary problems arising in fluid-structure interaction scenarios.It is worthwhile to note that in our model the effect of the structure part, i.e. either rubber's mechanics (see e.g.[3]) and/or material's capacity to perform capillary transport (see e.g.[22]), is incorporated in the shape of the nonlinearity σ(•).
The paper has the following structure: In Section 2 we specify the used notation, technical assumptions, as well as a couple of of useful basic inequalities.The weak formulation of our model together with the FEM discretization in space are included in Section 3. We recall here also some results obtained earlier by us concerning the semi-discrete FEM approximation.The bulk of the paper is the error analysis of the fully discrete approximation of our concept of solution.This is the purpose of Section 4. Section 5 contains a couple of numerical experiments confirming the theoretical convergence rates.Our conclusions on the obtained estimates and ideas for further work for are listed in Section 6.

Notation. Basic inequalities. Technical assumptions
In this framework, standard notations for Sobolev and Bochner spaces are used.An introduction to Sobolev and Bochner spaces as well as the usual notations, definition of norms and inner products can be found, for instance, in [1,16].For the convenience of writing, we denote by • and (•, •) the norm, and respectively, the inner product in L 2 (Ω).Furthermore, • ∞ refers to the norm of L ∞ (Ω).We also use the notation ( ) to indicate the derivative with respect to time variable.For the benefit of reader, we collect a few elementary inequalities that we frequently use in this context.
Throughout this paper, the involved parameters are assumed to fulfill the following conditions: (A1) a 0 , H, D, s 0 , T f are positive constants.
, where b * and b * are positive constants.
Proof.We refer the reader to Theorem 3.3 and Theorem 3.4 in [17] for a way to ensure the global solvability of the problem and continuous dependence estimates of the solution with respect to the initial data..
We discretize the fixed domain Ω = (0, 1) as follows.Let N ∈ N be given.We set 0 = y 0 < y as a finite dimesnional subspace of H 1 (0, 1).Here P 1 represents the set of polynomials of degree one.We define the interpolation operator where {φ i } N −1 i=0 are a set of basis functions for the space V k .Here the function I k u is called the Lagrange interpolant of u of degree 1; for more details see e.g.[18,5].Lemma 3.1.Take θ ∈ [ 1 2 , 1) and ψ ∈ H 2 (0, 1).Then there exist strictly positive constants γ 1 , γ 2 and γ 3 such that the Lagrange interpolant I k ψ of ψ satisfies the following estimates: Proof.The inequalities (i) and (ii) are standard results, see for instance p.61 in [18] and p. 90 in [5] for details on their proof.The proof of (iii) and (iv) follows from using the interpolation inequality, (i) and (ii).We refer the reader to Lemma 2.1 in [25] for details on their proof..
The continuous in time finite element approximation u k and h k of u and h is now represented by the mappings Our concept of solution to the semi-discrete problem is defined next.Definition 3.2.We call the couple (u k , h k ) a weak solution to the semi-discrete formulation if and only if there is a S T := (0, T ) (for some T > 0) such that and for all τ ∈ S T it holds Theorem 3.2.Assume (A1)-(A5) hold.Then it exists a unique solution on a time S T := (0, T ) for T ∈ (0, T ] in the sense of Definition 3.2.Furthermore, there exists a constant c > 0 Proof.We refer the reader to Theorem 4.2 in our previous work [25] for details of the proof. We anticipate in Proposition 3.1 a regularity result that is going to be useful when proving convergence rates for the proposed fully-discrete approximation.Mind though that the convergence of our scheme holds for much less regularity than stated.Preposition 3.1.(Higher regularity).Assume that (A1)-(A5) hold together with b ∈ W 2,2 (0, T ) and u 0 ∈ H 2 (0, s 0 ).Then where p := h , w := ∂ τ u, and (u, h) is the weak solution in the sense of Definition 3.
Proof.One differentiates with respect to time ( 10)-( 15) and proves for the resulting system of equations a statement as Theorem 3.1 in [17].Note that the embedding

Fully discrete error analysis
In this section, we present firstly a fully discrete numerical scheme of the problem ( 10)-( 15) and then perform the error analysis.Let n, M ∈ N. Concerning the discretization in time, we decompose the interval (0, T ] in to M subintervals.Let ∆τ := T /M be a step size of the time variable.Define At any time level τ n , we denote u(τ n , x) for x ∈ [0, 1] and h(τ n ) by u n and h n respectively.Furthermore, we denote the fully discrete approximation of u n and h n by U n and W n respectively.We use the following notations for the time derivative: Using these notations, the fully discrete problem is formulated as follows: Find the pair (U n+1 , W n+1 ) ∈ V k ×R + such that the following system holds for all n ∈ {0, • • • , M − 1}: where U 0 (0) is an appropriate approximation of the initial condition u 0 in V k .
Lemma 4.1.Assume that (A1)-(A5) hold.Then there exist T ∈ (0, T ] and positive constants K, u, ū and W such that for all n ∈ N the following inequalities hold true: Proof.The property (i) is built in the concept of solution cf.Definition 3.2, while (iii) is a direct consequence of (i) and (ii).The fact that the solution to our problem satisfies a weak maximum principle makes us confident that (ii) holds as well.One way to prove such statement directly would be to ensure that a discrete maximum principle holds in our situation; we refer for instance to the working techniques used in [11].
Preposition 4.1.Assume that (A1)-(A5) hold.Then for sufficiently small ∆τ the solution (U n , W n ) of ( 26)-(29) satisfies the following estimates: The constants K * and K 1 (K * , A 0 , W 0 ) depend on data, parameters, as well as on the quality of the approximation of the term and Cauchy-Schwarz's inequality, we get on both sides and multiplying by 2∆t gives where , the inequality (31) can be written as follows: where Using a 3 > 1 in (32), we get the following inequalities: . . .
Adding (33) − (36), we obtain For sufficiently small ∆τ and n ∈ {0, 1, • • • , M − 1}, we can write It follows from (37) This completes the proof of (i).To prove (ii), we use (27) to get Furthermore, we get Repeating the same procedure to bound |W n | gives Remark 4.1.The statement (i) in Proposition 4.1 can be seen as the discrete version of (25).It provides information regarding the stability of our fully discrete scheme.It turns out that the result holds for a sufficiently small ∆τ for a fixed n.The restriction on time step size indicates the our discretization scheme is possibly conditionally stable.
We now discuss the existence of solution to the fully discrete problem ( 26)-(29).Proof.Starting from W 0 and U 0 , the existence of W 1 comes from the existence of σ(W 0 ).For the existence of the solution for the concentration profile U n+1 , we use the application of the Brouwer fixed point theorem.
We refer the reader to p. 206 in [14], Lemma 4.2 in [4], as well as to Lemma 1.4 in [26] for more functional analytic details.The space V k is a separable as a subspace of H 1 (0, 1).Let {w 1 , w 2 , • • • , } be an orthogonal basis of H 1 (0, 1) and an orthonormal basis for L 2 (0, 1).For each fixed integer m ≥ 1, we define an approximation solution U n+1 of (26) by We define where U n+1 is defined in (39).Thus, (40) will have a solution if there exists a ζ such that for the first term, integration by part formula for the second term on right hand side of (41), we get, We choose U n+1 2 = R 2 large enough such that (F (ζ), ζ) > 0. This completes the proof.Proof.Assume by contradiction that (U n+1 , W n+1 ) and ( Ū n+1 , W n+1 ) are two solutions satisfying ( 26)-( 29) with the same initial data.Then we get the following holds for all ϕ k ∈ V k : Furthermore, we have We will show that these two solutions must coincide.We use the method of induction to prove U n+1 = Ū n+1 and W n+1 = W n+1 for all n ∈ {−1, 0, 1, •, }.Obviously, it holds for n = −1.Assume that the statement holds for an arbitrarily fixed n ∈ N, i.e, U n = Ū n .It remains to show U n+1 = Ū n+1 and W n+1 = W n+1 .Subtracting (45) from (44), and using the induction hypothesis gives By repeating the same process, it yields It now remains to show U n+1 = Ū n+1 .We subtract (43) from (42) and use the induction hypothesis to obtain Using Lemma 4.1 and Young's inequality, we get Choosing ξ ≤ ( W 0 K W ) 2 and ∆τ ≤ 1/c ξ in (47) gives U n+1 − Ū n+1 2 ≤ 0. It implies U n+1 = Ū n+1 a.e. in (0, 1).This completes the proof.
We next analyze error estimates of our fully discrete scheme for the concentration profile and the position of the moving boundary.To estimate the errors e n := U n − u n and e n 1 := W n − h n , we decompose e n into two parts: with ψ n := U n − I k u n and ρ n := I k u n − u n .Here I k u n is a Lagrange interpolation of u n defined in Lemma 3.1.In the rest of section we derive error estimates for the fully discrete scheme ( 26)-(29).To begin with, we perform the error bound for ψ n+1 in the following theorem.Then there exists a constant K > 0 such that the following inequality holds for sufficiently small ∆τ : Proof.Subtracting ( 21) from ( 26), we obtain the following identity: which holds for all ϕ k ∈ V k .Using U n = ψ n + ρ n + u n and arranging conveniently the terms in (50) yields where we introduce the following notations: , Before proceeding further, we collect two useful estimates in following auxillary Lemma 4.2.
We now prove (ii).Using the interpolation inequality (8) yields We now estimate the terms I 1 − I 5 as follows.By adding and subtracting appropriate terms, we get Using Lemma 4.2, it yields After re-arranging the term conveniently, I 2 becomes It holds The bound on |I 3 | follows from the Cauchy-Schwartz inequality, Young's inequality and Lemma 3.1 To deal with I 4 , we start by re-arranging the term in a more convenient way Using Lemma 4.2, we obtain By adding and subtracting appropriate terms in I 5 , we get We now claim the following holds: By the Taylor expansion of u n+1 around τ n with integral reminder yields To prove (53), the fundamental theorem of calculus gives Integrating by parts in the last term of (54) gives Using again the fundamental theorem of calculus for the last but one term in (55) leads to This proves (53).It is worth mentioning that (53) resembles the application of Taylor's approximation with integral reminder for the function u n around τ n+1 under the assumption u ∈ L 1 (τ n , τ n+1 ) (see Theorem 1.3 in [8]).With the help of (53) and Preposition 3.1, we can estimate the second term as follows: The last but one term in I 5 can be estimated as follows: Finally, we bound I 5 by , we also note the following estimate holds: We now consider the equations corresponding to the position of the moving boundary.We write Multiplying ( 59) by e n+1 Using (58) in (60), we get By using Taylor series expansion and Preposition 3.1, we estimate the second term in (61) as follows: Using Lemma 4.2, (62) and Young's inequality in (61) yields Taking ϕ k = ψ n+1 in (51) and adding the result to (63), we get the following estimate: Adding (1 − K 6 (ξ, ξ) − K 5 (ζ, ξ)) ∂ψ n ∂y to both sides of (64) and multiplying the result by 2∆τ , one gets With the notation: )) the inequality (65) can be rewritten as By dividing both sides of (66) by (1 − ∆τ K 4 ), we get the following inequalities: . . . where Adding (67)-(70), we obtain . For sufficiently small ∆τ , we can write This completes the proof.
Theorem 4.4.Assume the assumption of Theorem 4.3 holds.Then there exists a constant Proof.Relying on Theorem 4.3 and using the relation (a + b) 2 ≤ 2(a 2 + b 2 ), it holds for all 0 ≤ n ≤ M , where K 7 is as used in (71).
Remark 4.2.If we assume only u 0 ∈ H 2 (0, s 0 ) and b ∈ W 1,2 (0, T ), then we can only prove at most √ ∆τ convergence rate in time.This fact is a consequence of handling the estimate (56).

Numerical experiments
In this section, we present numerical results to substantiate the theoretically obtained order of convergence in space and time proposed in Section 4. We solve ( 21)-( 24) by using the method of lines; for more details see, for instance, [18].We refer the interested reader to [24] for more technical information on experimental data, implementation of the numerical method, and additional simulation results.The reader may also consult [25] for a priori and a posteriori error estimates of our semi-discrete finite element approximation.Here we take the final time T f to be 10 minutes.We discretize the domain in a uniform mesh size and use piecewise linear functions as basis for the subspace V k .The values of parameters are taken to be s 0 = 0.01 (mm), m 0 = 0.1 (gram/mm 3 ) and b = 1 (gram/mm 3 ).We take the value 3.66 × 10 −4 (mm 2 /min) for the diffusion constant D, 0.564 (mm/min) for absorption rate β and 2.5 for Henry's constant H.We choose σ(s(t)) = s(t)/10 (gram/mm 3 ) and a 0 = 50 (mm 4 /sec/gram).This specific choice of parameters is taken from [24,25].To test numerically the convergence order in space, we fix the uniform time step size ∆t = 0.0001.As we do not know the exact solution, we compute the finite element approximation on a fine mesh size with the total number of node N to be 1280 for a reference solution.To compute errors, we compare our reference solution against finite element approximations corresponding to different mesh sizes with increasing total number of nodes as 20, 40, 80, 160, 320, and 640.We calculate the convergence order in space based on any two consecutive calculations of discrete errors.The obtained errors and convergence orders in space are listed in Table 1.We show in Figure 1 the computed convergence order in space for the approximation of the moving boundary position and of the concentration profile.
To capture numerically the convergence order in time, we fix the total number of space node N to be 320 and choose ∆t = 10 −3 .We compute the finite element approximation on a fine time step size ∆t = ∆t/64 for a reference solution.We then compute approximations for different time step sizes ∆t, ∆t/2, ∆t/4, ∆t/8, ∆t/16, ∆t/32 and then calculate the order of convergence in time.The errors and convergence orders in time are listed in  These numerical results are in agreement with the convergence orders proven in Section 4.

Conclusion
We shown a fully discrete scheme for the numerical approximation of a moving boundary problem describing diffusants penetration into rubber.The proposed scheme utilizes the Galerkin finite-element method in space and the backward Euler method in time.By using Brouwer's fixed-point theorem, we were able to prove the existence of solution to the fully discrete problem.As main result, we obtained a priori error estimates for the mass concentration of diffusants as well as for the position of the moving boundary.The convergence turns to be of first order in both space and time for the approximation of the mass concentration of diffusants as well as for the approximation of the position of the moving boundary.Finally, we illustrated numerically the order of convergence in space and time to confirm the theoretically obtained results.It could be that the order of convergence in time can be improved by selecting other, perhaps better suited, time discretization schemes than the backward Euler one.Because of the presence of the moving boundary, the convergence order is space is lower than what we would expect for the finite element approximation of standard linear parabolic problems.This is in agreement with what is stated in the literature concerning the numerical approximation of one-dimensional moving-boundary problems.

Figure 1 :
Figure 1: Convergence order in space when time step size ∆t = 10 −4 is fixed.Dash lines are lines of slope 1. Left: Log log scale plot of error on the boundary max 0≤n≤M |W n − h n |.Right: Log log scale plot of error on the concentration max 0≤n≤M U n − u n L 2 (0,1) .

Figure 2 :
Figure 2: Convergence order in time when space mesh size is fixed with N = 320.Dash lines are lines of slope 1. Left: Log log scale plot of error on the boundary max 0≤n≤M |W n − h n |.Right: Log log scale plot of error on the concentration max 0≤n≤M U n − u n L 2 (0,1) .

Table 1 :
1) Errors and convergence orders in space with fixed ∆t = 10

Table 2 .
We show in Figure2the computed convergence order in time for the approximation of both the position of the moving boundary and corresponding concentration profile.

Table 2 :
Errors and convergence orders in time with fixed N = 320.