Mathematical analysis of some diffusion problems associated to the modeling of surfactant compounds at the air-water interface

MATHEMATICAL ANALYSIS OF SOME DIFFUSION PROBLEMS ASSOCIATED TO THE MODELING OF SURFACTANT COMPOUNDS AT THE AIR-WATER INTERFACE Cristina N´u˜nez Garc´ıa Departamento de Matem´atica Aplicada Universidade de Santiago de Compostela Ph.D. Dissertation September 2013 ii Don Jos´e Ram´on Fern´andez Garc´ıa, profesor Titular de Universidad del Departamento de Matem´atica Aplicada I de la Universidad de Vigo, y Do˜na Mar´ıa del Carmen Mu˜niz Casti˜neira, profesora Titular de Universidad del Departamento de Matem´atica Aplicada de la Universidad de Santiago de Compostela, informan de que la memoria titulada: MATHEMATICAL ANALYSIS OF SOME DIFFUSION PROBLEMS ASSOCIATED TO THE MODELING OF SURFACTANT COMPOUNDS AT THE AIR-WATER INTERFACE fue realizada bajo su direcci´on por Do˜na Cristina N´u˜nez Garc´ıa, estimando que la interesada se encuentra en condiciones de optar al grado de Doctor en Ciencias Matem´aticas, por lo que solicitan que sea admitida a tr´amite para su lectura y defensa p´ublica. En Santiago de Compostela, a 30 de Septiembre de 2013. Los directores: Prof. Dr. Jos´e R. Fern´andez Garc´ıa Prof. Dra. M. del Carmen Mu˜niz Casti˜neira La doctoranda: Cristina N´u˜nez Garc´ıa A mis padres, Antonio y Elena. Agradecimientos Quisiera comenzar agradeciendo profundamente a mis directores de tesis, Jose Ram´on Fern´andez y Mar´ıa del Carmen Mu˜niz, por abrirme las puertas al mundo de la investigaci´on y por haberme dado la oportunidad de trabajar a su lado. Valoro de manera inestimable sus ´animos, su dedicaci´on, sus ense˜nanzas y, sobre todo, su infinita paciencia. Tambi´en me gustar´ıa tener unas palabras especiales de agradecimiento para los profesores Stanislaw Mig´orski y Piotr Kalita, por todo lo que aprend´ı de ellos durante mi estancia en Polonia. Adem´as, quisiera destacar que, gracias a su hospitalidad y simpat´ıa, me sent´ı como en casa en todo momento. Quisiera mostrar mi gratitud al profesor Luis Garc´ıa R´ıo, del Departamento de Qu´ımica F´ısica de la Universidad de Santiago de Compostela, por su amabilidad y su generosidad en las explicaciones e interpretaciones qu´ımicas de este trabajo. Tambi´en me gustar´ıa agradecer a su alumno de tesis, Serxio Ith Arias, por la ayuda que me ha prestado en numerosas ocasiones. Me gustar´ıa hacer extensible este agradecimiento a todos los miembros del Departamento de Matem´atica Aplicada de la Universidad de Santiago de Compostela. Ha sido un privilegio poder compartir con vosotros esta experiencia. No quisiera olvidarme en este apartado de todos mis compa˜neros de facultad: Javi, Miguel, ´ Angel, Jorge, Saray, Miriam, Nizom,... en especial de los del despacho por todo lo que vivimos juntos. Ser´an inolvidables todas las aventuras que compart´ı con vosotros, desde los momentos de trabajo hasta los caf´es de la ma˜nana, las cenas, las tardes de piscina en casa de ´ Angel, e incluso los acalorados debates. Tambi´en quisiera darles las vii viii gracias a todos los compa˜neros con los que empec´e mi “andaina polas matem´aticas”, a los que acabaron la carrera y a los que se quedaron por el camino: Lara, Bea, Elena, Alba, Javi, Jacobo, Natalia, Andrea, D´eborah,... y al resto de amigos que siempre est´an ah´ı: Marta, Iv´an, Luis, Isra, Rebe, Noa,... Finalmente, me gustar´ıa dirigir estas ´ultimas palabras de agradecimiento a los verdaderos art´ıfices de este trabajo, mi familia. Por eso quiero agradecerles a mis padres, Antonio y Elena, el apoyo que siempre me han brindado; a mi madrina, a la par que abuela, el haberme cuidado tanto; a mis segundos padres, mis primos Antonio y Mari, los buenos consejos que siempre me han dado. Tambi´en quiero darles las gracias a mis primos Ana y Toni, especialmente a Ana porque siempre tiene palabras de ´animo y tranquilizadoras en los momentos que m´as lo necesito. Adem´as, quiero mencionar especialmente a Ant´ıa, mi ahijada, a la que tuve que esperar tres largos meses para conocer ya que su nacimiento “me pill´o” por Polonia. No me quisiera olvidar de dedicar unas palabras de gratitud a Samu, por todo su apoyo incondicional. Sin todos ellos, realizar este trabajo hubiera sido imposible. Unha aperta para todos! Contents Introduction 1 1 Mixed kinetic-diffusion model for the Henry isotherm 11 1.1 The mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . 12 1.2 Weak formulation of the problem . . . . . . . . . . . . . . . . . . . . . 14 1.3 Existence and uniqueness result . . . . . . . . . . . . . . . . . . . . . . 16 1.4 Fully discrete approximations: numerical analysis . . . . . . . . . . . . 21 1.5 Numericalresults.............................. 35 1.5.1 First example: numerical convergence . . . . . . . . . . . . . . . 36 1.5.2 Second example: simulation of hexanol . . . . . . . . . . . . . . 38 1.5.3 Third example: simulation of heptanol . . . . . . . . . . . . . . 42 2 Diffusion-controlled model with Langmuir isotherm 45 2.1 The mathematical model . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.2 Weak formulation of the problem . . . . . . . . . . . . . . . . . . . . . 49 2.3 Existence and uniqueness results . . . . . . . . . . . . . . . . . . . . . . 49 2.4 Analysis of a semi-discrete problem . . . . . . . . . . . . . . . . . . . . 75 2.5 Fully discrete approximations: a priori error estimates . . . . . . . . . . 81 2.6 Numericalresults.............................. 91 2.6.1 First example: numerical convergence . . . . . . . . . . . . . . . 93 2.6.2 Second example: simulation of propanol . . . . . . . . . . . . . 93 2.6.3 Third example: simulation of sodium dodecylsulfate . . . . . . . 95 ix 6Introduction we have the following formulation of the problem (see [7, 35]): ∂c ∂t(t, x)−D∂2c ∂x2(t, x)=0, t ∈(0, T), x ∈(0, l),(1) with boundary conditions: D∂c ∂x(t, 0) = dΓ dt (t), t > 0,(2) c(t, l) = cb, t > 0,(3) and initial conditions: c(0, x) = c0(x), x ∈(0, l),(4) Γ(0) = Γ0.(5) In this system of equations, the positive constants cband Ddenote the bulk concentration and the diffusion coefficient, respectively. We note that, in this thesis, we deal with surfactant solutions below their cmc, so a constant diffusion coefficient can be assumed (see [46]). Besides, c0(x) is a function defined in [0, l], which equals cbat x=l and Γ0is a nonnegative constant. Then, the transport of molecules from the bulk of the solution to the surface is modeled by equations (1)-(5). Equation (1), that describes the diffusion in the bulk of the solution considering a finite diffusion length, is obtained from the general transport equation by neglecting the convective term since, as we said previously, we are in the framework of quiescent surfactant solutions. Boundary condition (2) describes the surfactant flux from the subsurface to the surface (adsorption) and vice versa, from the surface to the subsurface (desorption). Moreover, we assume that the boundary x=lis kept at a constant concentration, cb, during the process, so a Dirichlet boundary condition, given by expression (3), is imposed there. In terms of the initial conditions, we assume that, at the beginning, the surface concentration is equal to Γ0and the concentration in the bulk of the solution is given by the function c0(x), x ∈[0, l]. Introduction 7 Therefore, given l, T, D, cb, c0(x) and Γ0, the problem consists in finding both the surface and subsurface concentrations. Since the surface concentration, Γ(t), is also an unknown of the system, an additional condition must be given in order to close the problem. In this sense, the additional condition, that is coupled to the system of equations (1)-(5) by means of the boundary condition at the subsurface, is established by the adsorption mechanism; so either the diffusion-controlled model or the mixed kinetic-diffusion one has to be used. When considering a diffusion-controlled model, a thermodynamic adsorption isotherm states the dependence between the surface and subsurface concentrations. Three of the most commonly used isotherms in the literature are (see [7, 16, 28, 29, 34]): •The Henry isotherm: it is the simplest isotherm and it is only valid for low surface concentrations since it does not take into account interactions between adsorbed molecules. Besides, in this model there is no limit on the surface concentration (see [7, 16]). This isotherm establishes a linear dependence between the surface and subsurface concentrations, Γ(t) = KHc(t, 0), t ≥0,(6) where KHis the Henry equilibrium adsorption constant, which is a measure of the surface activity of the surfactant. •The Langmuir isotherm: it defines a nonlinear dependence between the surface and subsurface concentrations, Γ(t) = Γm KLc(t, 0) 1 + KLc(t, 0), t ≥0,(7) being Γmand KLthe maximum surface concentration and the Langmuir equilibrium adsorption constants, respectively. •The Frumkin isotherm: as in the case of the Langmuir isotherm, this expression states a nonlinear dependence between both surface and subsurface concentra- 8Introduction tions, Γ(t)=Γm KFc(t, 0) e−AΓ(t) Γm+KFc(t, 0) , t ≥0,(8) where KFis the Frumkin equilibrium adsorption constant and Ais a parameter that indicates if the adsorption is anticooperative or not. In the case of negative A, then the adsorption is anticooperative; that is to say, it becomes more difficult as the coverage of the surface increases. The case of Apositive describes the existence of cohesive intermolecular forces, which increases the surface coverage and which makes that the desorption rate decreases (see [29]). Finally, if A= 0, then Langmuir isotherm is actually recovered. On the other hand, in mixed kinetic-diffusion models, a kinetic expression identifies the rate of change of the surface concentration with the balance between the adsorption and desorption rates. The most studied equations are (see [7]): •The linear kinetic model: in which the rate of adsorption is proportional to the subsurface concentration, while the rate of desorption is proportional to the surface concentration, dΓ dt (t) = ka Hc(t, 0) −kd HΓ(t), t > 0,(9) where ka Hand kd Hare the adsorption and desorption constants, respectively. •The Langmuir-Hinshelwood kinetic model: in which the rate of adsorption depends on the subsurface concentration, but also on the fraction of empty space at the surface, dΓ dt (t) = ka Lc(t, 0)1−Γ(t) Γm−kd LΓ(t), t > 0,(10) being ka Land kd Lthe adsorption and desorption constants for the LangmuirHinshelwood kinetic model, respectively. •The modified Langmuir-Hinshelwood kinetic model: the modification of the Lang- Introduction 9 muir-Hinshelwood equation was proposed by Chang and Franses in 1992 (see [6]), because the previous kinetic equation did not fit the experimental data of some surfactants well. With this modification, better results for those surfactants were obtained (see [6, 7]), dΓ dt (t) = ka Lc(t, 0)1−Γ(t) Γme−BΓ(t) Γm−kd LΓ(t)e−BΓ(t) Γm, t > 0,(11) where the real constant Bis an empirical parameter. In the case that B= 0, the classical Langmuir-Hinshelwood expression (10) is recovered. Once the problem is solved, then both surface and subsurface concentrations are known. Now, the surface tension is calculated through the so-called equation of state, that is deduced from the Gibbs adsorption isotherm. This equation relates the surface concentration, Γ(t), to both surface tension, eγ(t), and subsurface concentration, c(t, 0), at a constant temperature, and its expression is given by: Γ(t) = −1 n R θ ∂eγ(t) ∂ln c(t, 0)θ , where Ris the gas constant, nis a constant which is equal to 1 for non-ionic surfactants and θis the temperature. As we mentioned previously, in the chemical literature, we can find several publications devoted to solve numerically some of these problems (see, for instance, [5, 6, 7, 30]). However, to our knowledge, none of those works deal with the study of the existence and uniqueness of weak solutions for both the variational formulations of the problems and their discrete approximations. Furthermore, error estimates for the differences between the continuous solutions and the discrete ones, as well as convergence order results, were not introduced yet. Therefore, the main contribution of this thesis is that we perform variational and numerical analyses of a diffusion problem coupled with a dynamical boundary condition. The outline of this Ph.D. Thesis is as follows. In Chapter 1, we concern on the variational and numerical formulation of the diffusion problem considering the linear 10 Introduction kinetic model. The existence of solution to the weak problem is proved by formulating an auxiliary problem followed by the application of the Banach fixed-point theorem. Uniqueness of solution of the weak problem is also shown. Moreover, fully discrete approximations of the problem are obtained by using the finite element method and a combination of both backward and forward Euler schemes. Under additional regularity conditions, an error estimate result is obtained from which the linear convergence is deduced. Some numerical simulations, in order to show the accuracy of the algorithm and its behavior for a commercial surfactant, are provided. We point out that this chapter has given rise to [20] and [22]. In Chapter 2, we focus on the diffusion problem regarding the Langmuir isotherm. For this problem, the existence of a unique weak solution is proved by using the Rothe’s method and fixed-point techniques. Besides, a semi-discrete problem in time is analyzed for which we get the linear convergence under additional regularity conditions. Also, following the same ideas as in Chapter 1, fully discrete approximations of the problem are presented. An error estimate result is proved from which the linear convergence is followed under suitable regularity conditions. We indicate that the work introduced in this chapter has been collected in [11] and [21]. Finally, in Chapter 3, the diffusion problem together with either the LangmuirHinshelwood equation or the modified Langmuir-Hinshelwood one is taken into account. The main results of this chapter deal with the existence and uniqueness of weak solution for the truncated versions associated to both problems. Their proofs are obtained by dividing the truncated problems into two auxiliary problems and the application of the Schauder fixed-point theorem. A numerical analysis is performed using some of the ideas already applied in the previous chapters, and numerical simulations that exhibit the accuracy and the behavior of the algorithm for some surfactants are also shown. The work presented in this chapter has led to the manuscripts [23] and [24]. Chapter 1 Mixed kinetic-diffusion model for the Henry isotherm In this chapter, we describe the adsorption-desorption dynamics of a surfactant solution at the air-water interface by considering a linear mixed kinetic-diffusion model. It is given by the simplest kinetic expression modeling this behavior, which establishes that the rate of change of the surface concentration is related to the balance between the amount of surfactant molecules that migrate from the subsurface to the surface and the amount of surfactant molecules that move from the surface to the subsurface. From a mathematical point of view, the whole dynamic process is modeled by a coupled nonlinear system of a parabolic equation, for the description of the diffusion dynamics, and an ordinary differential equation, for the adsorption-desorption mechanism. Here, we prove the existence and uniqueness of solution to the weak problem by using classical results for linear parabolic equations and fixed-point techniques. Then, fully discrete approximations are obtained by using the finite element method and a hybrid combination of both backward and forward Euler schemes. An a priori error estimates result is presented from which, under adequate additional regularity conditions, the linear convergence of the algorithm is derived. Finally, numerical simulations are introduced to demonstrate the accuracy of the algorithm and the behavior for two commercially available surfactants. 11 12 Mixed kinetic-diffusion model for the Henry isotherm 1.1 The mathematical model Let us denote by xthe distance from the interface and c(t, x) the concentration of surfactant at time t∈[0, T] and point x∈[0, l]. The boundary x= 0 of the spatial interval corresponds to the location of the subsurface, an imaginary layer between the region in which only diffusion takes place and the domain where only adsorption-desorption occurs. Denoting by Γ(t) the time-dependent surface concentration and taking into account the Fick’s law, we consider the diffusion partial differential equation: ∂c ∂t(t, x)−D∂2c ∂x2(t, x) = 0, t > 0, x ∈(0, l),(1.1) together with the boundary conditions (see [7, 35]): D∂c ∂x(t, 0) = dΓ dt (t), t > 0,(1.2) c(t, l) = cb, t > 0,(1.3) and the initial conditions: c(0, x) = c0(x), x ∈(0, l),(1.4) Γ(0) = Γ0.(1.5) In equations (1.1)-(1.3), the positive constants Dand cbrepresent the diffusion coefficient and the bulk concentration, respectively. Besides, c0(x) is a function defined in [0, l] which is equal to cbon x=l. We remark that the time-dependent surface concentration, Γ(t), actually becomes an unknown of the system and then an additional condition must be given in order to close the problem. Hereinafter, in this chapter, we consider a linear mixed kinetic-diffusion model, given by the simplest kinetic expression modeling the mass transfer between the surface and subsurface at low concentrations which leads to the following ordinary differential equation (see [7, 44]): 1.1. The mathematical model 13 dΓ dt (t) = ka Hc(t, 0) −kd HΓ(t), t > 0,(1.6) where ka Hand kd Hare the adsorption and desorption constants, respectively. This expression identifies the rate of change of the surface concentration with the balance between the adsorption and desorption rates. Moreover, it leads to the Henry isotherm at equilibrium (see [19]). Remark 1.1 We note that Henry’s isotherm is the simplest equation for describing the adsorption dynamics. It establishes a linear dependence between the subsurface and surface concentrations, assuming that the surface concentration is proportional to the subsurface concentration. Its expression is given by: Γ(t) = KHc(t, 0), t ≥0, KHbeing the Henry equilibrium isotherm. At equilibrium or steady-state, dΓ/dt = 0 and, from equation (1.6), the classical Henry’s isotherm is recovered with KH=ka H/kd H. The study of the surfactant behavior at the air-water interface accounting the Henry isotherm was studied in [19]. In that work, the existence and uniqueness of solution to the weak problem was proved. Moreover, fully discrete approximations of the problem were presented for which error estimates were obtained, and under adequate additional regularity conditions, the linear converge of the algorithm was derived. Finally, some numerical simulations were shown to demonstrate the accuracy of the algorithm and the behavior of the solution. Assuming that the solution is regular enough, the previous ordinary differential equation, (1.6), together with the initial condition (1.5) can be straightforwardly integrated to obtain: Γ(t)=Γ0e−kd Ht+ka He−kd HtZt 0 ekd Hτc(τ, 0)d τ. 14 Mixed kinetic-diffusion model for the Henry isotherm Therefore, boundary condition (1.2) reads D∂c ∂x(t, 0) = ka Hc(t, 0) −φ(t, c(·,0)), t > 0,(1.7) where φ(t, ζ) = kd HΓ0e−kd Ht+kd Hka He−kd HtZt 0 ekd Hτζ(τ)d τ. (1.8) We emphasize that boundary condition (1.7) determines a non-local boundary condition in time since, for the construction of the flux at time t, the values of the subsurface concentration at previous times are also required. We are now concerned in analyzing problem (1.1), (1.3) and (1.4), together with the new boundary condition (1.7). Moreover, for the sake of clarity in the presentation of this chapter, and in order to simplify the calculations of the following sections, we assume, without loss of generality, that cbequals zero and so a homogeneous boundary condition is imposed on the right end of the spatial interval. 1.2 Weak formulation of the problem Before establishing the weak formulation of the problem, we introduce the notation we use hereinafter in this manuscript. Let Vbe the Hilbert space V={v∈H1(0, l); v(l) = 0}, endowed with the inner product and the associated norm given, respectively, by ((v, w)) = Zl 0 ∂v ∂x ∂w ∂x dx, kvkV= ((v, v))1/2. As usual, we denote by V0the dual space to Vand by h·,·i the scalar product for the duality V0,V. Furthermore, we recall the inner product in H=L2(0, l) given by (v, w)H=Zl 0 v(x)w(x)dx, 1.2. Weak formulation of the problem 15 with associated norm kvkH= (v, v)1/2 H. Moreover, we consider the Hilbert space V= L2(0, T;V) with dual space V0=L2(0, T;V0) together with W2(0, T) = {v∈ V;∂v ∂t ∈ V0}, where the time derivative is understood in distributional sense (see [40]). Furthermore, for a Banach space Xand a nonnegative integer r, here Cr([0, T]; X) denotes the space of rtimes continuously differentiable functions from [0, T] to X. We denote by γ0:H1(0, l)→Rbe the trace operator on x= 0 given by γ0(v) = v(0). From the continuity of the trace operator (see, for instance, Theorem 3.9.34 in [12]), it follows that |γ0(v)| ≤ Ctr kvkVfor all v∈Vwith Ctr =kγ0kL(V,R).(1.9) Now, let cbe a smooth function which solves the problem given by equations (1.1), (1.3), (1.4) and (1.7). Multiplying expression (1.1) by a smooth function zdefined in [0, l], such that z(l) = 0, integrating in (0, l) and using the integration by parts formula, we obtain, for a.e. t∈(0, T), Zl 0 ∂c ∂t(t, x)z(x)dx +Zl 0 D∂c ∂x(t, x)∂z ∂x(x)dx +D∂c ∂x(t, 0)z(0) = 0. Using equation (1.7), we find that, for a.e. t∈(0, T), Zl 0 ∂c ∂t(t, x)z(x)dx +Zl 0 D∂c ∂x(t, x)∂z ∂x(x)dx +ka Hc(t, 0) z(0) = φ(t, c(·,0)) z(0). Therefore, we have the following weak formulation of problem (1.1), (1.3), (1.4) and (1.7): Problem PH W. For a given c0∈H, find a function c∈W2(0, T) such that h∂c ∂t(t), viV0×V+D((c(t), v)) + ka Hγ0(c(t)) γ0(v) = φ(t, γ0(c)) γ0(v), for a.e. t∈(0, T),∀v∈V, (1.10) c(0) = c0.(1.11) 22 Mixed kinetic-diffusion model for the Henry isotherm Problem Phk H. Find chk ={chk n}N n=0 ⊂Vhsuch that chk 0=ch 0,(1.20) and, for n= 1, . . . , N and for all vh∈Vh, (δchk n, vh)H+D((chk n, vh)) + ka Hγ0(chk n)γ0(vh) = φhk n−1γ0(vh),(1.21) where ch 0∈Vhis an appropriate approximation of the initial condition c0and φhk n−1=kd HΓ0e−kd Htn+kd Hka Hk n−1 X j=0 ekd H(tj−tn)γ0(chk j).(1.22) For Problem Phk H, we have the following result. Theorem 1.2 Assume that the hypotheses of Theorem 1.1 hold. Then, Problem Phk H has a unique solution. Proof. Let us consider the bilinear form aH:V×V→Rgiven by aH(u, v) = Zl 0 u v dx +D k Zl 0 ∂u ∂x ∂v ∂x dx +k ka Hγ0(u)γ0(v), and the linear form LH:V→Rdefined by LH(v) = Zl 0 chk n−1v dx +k φhk n−1γ0(v). The bilinear form aHis continuous on V×Vand coercive on V. Indeed, taking into account both H¨older and trace inequalities (see (1.9)) and that the norms k·kH1(0,l) and k·kVare equivalent on the space V, we have |aH(u, v)| ≤ kukHkvkH+D kkukVkvkV+k ka HC2 trkukVkvkV ≤max{1, D k, k ka HC2 tr}kukH1(0,l)kvkH1(0,l)≤M∗kukVkvkV,∀u, v ∈V, 1.4. Fully discrete approximations: numerical analysis 23 where M∗is a positive constant large enough, and, aH(v, v) = Zl 0 v2dx +D k Zl 0∂v ∂x2 dx +k ka H(γ0(v))2 ≥min{1, D k}kvk2 H1(0,l)≥αkvk2 V,∀v∈V, αbeing a positive constant small enough. Moreover, considering both H¨older and trace inequalities again, the equivalence between the norms k·kH1(0,l)and k·kVon Vand that the function φhk n−1is bounded for n∈ {1,2, . . . , N}, we deduce that |LH(v)|≤kchk n−1kHkvkH+k|φhk n−1||γ0(v)| ≤ M∗∗kvkV,∀v∈V, where M∗∗ is a positive constant large enough. Then, LHis continuous on V. Consequently, we can apply Lax-Milgram theorem and the result follows.  In the sequel, we will derive an error estimate for the difference cn−chk nassuming the following additional regularity: c∈ C([0, T]; V)∩C1([0, T]; H).(1.23) Taking v=cn−vh∈Vin equation (1.10) at time t=tn, we find that, for n= 1,2, . . . , N, ∂c ∂t(tn), cn−vhH +D((cn, cn−vh)) + ka Hγ0(cn)γ0(cn−vh) = φnγ0(cn−vh), where φn=φ(tn, γ0(c)),and therefore, since the previous expression holds also for vh=chk n, it follows ∂c ∂t(tn), cn−chk nH +D((cn, cn−chk n)) +ka Hγ0(cn)γ0(cn−chk n)−φnγ0(cn−chk n) =∂c ∂t(tn), cn−vhH +D((cn, cn−vh)) +ka Hγ0(cn)γ0(cn−vh)−φnγ0(cn−vh). (1.24) 24 Mixed kinetic-diffusion model for the Henry isotherm On the other hand, taking vh−chk n∈Vhas a test function in (1.21) and writing vh−chk n=vh−cn+cn−chk nwe have, for all vh∈Vh, (δchk n, cn−chk n)H+D((chk n, cn−chk n)) +ka Hγ0(chk n)γ0(cn−chk n)−φhk n−1γ0(cn−chk n) = (δchk n, cn−vh)H+D((chk n, cn−vh)) +ka Hγ0(chk n)γ0(cn−vh)−φhk n−1γ0(cn−vh).(1.25) Subtracting now equations (1.24) and (1.25) and taking into account the linearity of the trace operator, we obtain, for all vh∈Vh, ∂c ∂t(tn)−δchk n, cn−chk nH +Dkcn−chk nk2 V+ka H|γ0(cn−chk n)|2 −(φn−φhk n−1)γ0(cn−chk n) =∂c ∂t(tn)−δchk n, cn−vhH +D((cn−chk n, cn−vh)) +ka Hγ0(cn−chk n)γ0(cn−vh)−(φn−φhk n−1)γ0(cn−vh). Taking into account that ∂c ∂t(tn)−δchk n=∂c ∂t(tn)−δcn+δcn−δchk nand after easy algebraic manipulations we find that, for all vh∈Vh, (δcn−δchk n, cn−chk n)H+Dkcn−chk nk2 V+ka H|γ0(cn−chk n)|2 =∂c ∂t(tn)−δchk n, cn−vhH +D((cn−chk n, cn−vh)) +(φn−φhk n−1)γ0(cn−chk n)−(φn−φhk n−1)γ0(cn−vh) +ka Hγ0(cn−chk n)γ0(cn−vh) + δcn−∂c ∂t(tn), cn−chk nH ,(1.26) where we recall that δcn= (cn−cn−1)/k. 1.4. Fully discrete approximations: numerical analysis 25 Moreover, using the following property of the divided differences: (δan−δbn, an−bn)H= (an−an−1 k−bn−bn−1 k, an−bn)H =1 kkan−bnk2 H−1 k(an−1−bn−1, an−bn)H,(1.27) equation (1.26) reads 1 kkcn−chk nk2 H+Dkcn−chk nk2 V+ka H|γ0(cn−chk n)|2 =∂c ∂t(tn)−δchk n, cn−vhH +D((cn−chk n, cn−vh)) +(φn−φhk n−1)γ0(cn−chk n)−(φn−φhk n−1)γ0(cn−vh) +ka Hγ0(cn−chk n)γ0(cn−vh) + δcn−∂c ∂t(tn), cn−chk nH +1 k(cn−1−chk n−1, cn−chk n)H,∀vh∈Vh.(1.28) Using both Cauchy-Schwarz and Cauchy inequalities, it follows that 1 2kkcn−chk nk2 H+D 2kcn−chk nk2 V+ka H|γ0(cn−chk n)|2 ≤∂c ∂t(tn)−δchk n, cn−vhH +D 2kcn−vhk2 V +(φn−φhk n−1)γ0(cn−chk n)−(φn−φhk n−1)γ0(cn−vh) +ka Hγ0(cn−chk n)γ0(cn−vh) + δcn−∂c ∂t(tn), cn−chk nH +1 2kkcn−1−chk n−1k2 H,∀vh∈Vh.(1.29) We will use now the following technical lemmas. Lemma 1.1 The following estimate holds: |φn−φhk n−1|2≤2I2 n+βk n−1 X j=0 |γ0(cj−chk j)|2, 26 Mixed kinetic-diffusion model for the Henry isotherm βbeing a positive constant independent of k,hand nand Inthe integration error given by In:= ka Hkd He−kd HtnZtn 0 ekd Hτγ0(c(τ)) dτ − n−1 X j=0 k ekd Htjγ0(cj).(1.30) Proof. First, we find that |φn−φhk n−1|=kd Hka He−kd Htn(Ztn 0 ekd Hτγ0(c(τ)) dτ − n−1 X j=0 k ekd Htjγ0(chk j)) =kd Hka He−kd HtnZtn 0 ekd Hτγ0(c(τ)) dτ − n−1 X j=0 k ekd Htjγ0(cj) + n−1 X j=0 k ekd Htjγ0(cj−chk j)), and therefore, reminding the definition of the integration error, In, given in (1.30), we obtain |φn−φhk n−1| ≤ In+kd Hka Hk n−1 X j=0 |γ0(cj−chk j)|. Using the following property (see [38]) (a+b)℘≤2(℘−1)+(a℘+b℘),for a, b ≥0,and ℘ > 0,(1.31) with ℘= 2, we deduce |φn−φhk n−1|2≤2I2 n+ 2 (kd Hka Hk)2 n−1 X j=0 |γ0(cj−chk j)|!2 .(1.32) Since n−1 X j=0 |γ0(cj−chk j)| ≤ n1/2 n−1 X j=0 |γ0(cj−chk j)|2!1/2 , 1.4. Fully discrete approximations: numerical analysis 27 from estimate (1.32) we have |φn−φhk n−1|2≤2I2 n+ 2 (kd Hka Hk)2n n−1 X j=0 |γ0(cj−chk j)|2, and, keeping in mind that k N =T, the result holds.  Lemma 1.2 There exist two positive constants, αand β,α < β, independent of h,k and nsuch that, using the notation an:= kcn−chk nk2 H+k n X j=0 [Dkcj−chk jk2 V+α|γ0(cj−chk j)|2], bn(vh) := 2  ∂c ∂t(tn)−δcn   2 H+kcn−vhk2 H+Dkcn−vhk2 V +β|γ0(cn−vh)|2+βI2 n, dn(vh) := (cn−chk n−(cn−1−chk n−1), cn−vh)H, it follows that an≤an−1+k(bn(vh) + β an)+2dn(vh),∀vh∈Vh, n ≥1.(1.33) Proof. Since ∂c ∂t(tn)−δchk n, cn−vhH =∂c ∂t(tn)−δcn, cn−vhH +1 k(cn−chk n−(cn−1−chk n−1), cn−vh)H ≤1 2  ∂c ∂t(tn)−δcn   2 H+1 2kcn−vhk2 H +1 k(cn−chk n−(cn−1−chk n−1), cn−vh)H,∀vh∈Vh, then estimate (1.29) implies that 28 Mixed kinetic-diffusion model for the Henry isotherm 1 2kkcn−chk nk2 H+D 2kcn−chk nk2 V+ka H|γ0(cn−chk n)|2 ≤1 2  ∂c ∂t(tn)−δcn   2 H+1 2kcn−vhk2 H+D 2kcn−vhk2 V+1 2kkcn−1−chk n−1k2 H +(φn−φhk n−1)γ0(cn−chk n)−(φn−φhk n−1)γ0(cn−vh) +ka Hγ0(cn−chk n)γ0(cn−vh) + δcn−∂c ∂t(tn), cn−chk nH +1 k(cn−chk n−(cn−1−chk n−1), cn−vh)H,∀vh∈Vh.(1.34) Moreover, using both H¨older and Cauchy inequalities, estimate (1.34) leads to the following estimate, for all vh∈Vh, 1 2kkcn−chk nk2 H+D 2kcn−chk nk2 V+ka H|γ0(cn−chk n)|2 ≤  ∂c ∂t(tn)−δcn   2 H+1 2kcn−vhk2 H+D 2kcn−vhk2 V+1 2kkcn−1−chk n−1k2 H +(φn−φhk n−1)γ0(cn−chk n)−(φn−φhk n−1)γ0(cn−vh) +ka Hγ0(cn−chk n)γ0(cn−vh) + 1 2kcn−chk nk2 H +1 k(cn−chk n−(cn−1−chk n−1), cn−vh)H.(1.35) Finally, since we have γ0(cn−chk n)γ0(cn−vh)≤1 2|γ0(cn−chk n)|2+1 2|γ0(cn−vh)|2, (φn−φhk n−1)γ0(cn−chk n)≤1 4ε|φn−φhk n−1|2+ε|γ0(cn−chk n)|2, −(φn−φhk n−1)γ0(cn−vh)≤1 2|φn−φhk n−1|2+1 2|γ0(cn−vh)|2, for a parameter ε > 0 assumed small enough, estimate (1.35) implies that, for all vh∈Vh, 1.4. Fully discrete approximations: numerical analysis 29 1 2kkcn−chk nk2 H+D 2kcn−chk nk2 V+α 2|γ0(cn−chk n)|2 ≤  ∂c ∂t(tn)−δcn   2 H+1 2kcn−vhk2 H+D 2kcn−vhk2 V+1 2kcn−chk nk2 H +β 2|φn−φhk n−1|2+β 2|γ0(cn−vh)|2 +1 2kkcn−1−chk n−1k2 H+1 k(cn−chk n−(cn−1−chk n−1), cn−vh)H, where αand βare generic positive constants, α < β, assumed to be small and large enough, respectively, independent of h,kand nand whose value may vary from line to line. Therefore, multiplying by 2kwe get, for all vh∈Vh, kcn−chk nk2 H+D k kcn−chk nk2 V+αk |γ0(cn−chk n)|2 ≤k2  ∂c ∂t(tn)−δcn   2 H+kcn−vhk2 H+Dkcn−vhk2 V+kcn−chk nk2 H +β|φn−φhk n−1|2+β|γ0(cn−vh)|2+kcn−1−chk n−1k2 H +2(cn−chk n−(cn−1−chk n−1), cn−vh)H.(1.36) Lemma 1.2 is now a consequence of estimate (1.36) and Lemma 1.1. Indeed, adding k n−1 X j=0 (Dkcj−chk jk2 V+α|γ0(cj−chk j)|2), in both sides of inequality (1.36) and using Lemma 1.1, (1.33) holds.  Consequently, from (1.33) we obtain an≤a0+ n X j=1 (k(bj(vh j) + β aj)+2dj(vh j)),∀{vh j}n j=1 ⊂Vh.(1.37) 30 Mixed kinetic-diffusion model for the Henry isotherm Since n X j=1 k bj(vh j)≤k N X j=1 bj(vh j)≤TM, where M= max 1≤j≤Nbj(vh j), estimate (1.37) reads an≤a0+T M + n X j=1 (k β aj+ 2 dj(vh j)),∀{vh j}n j=1 ⊂Vh.(1.38) Taking into account (1.20), we notice now that, for all {vh j}n j=1 ⊂Vh, n X j=1 dj(vh j)=(cn−chk n, cn−vh n)H−(c0−ch 0, c1−vh 1)H + n−1 X j=1 (cj−chk j, cj−vh j−(cj+1 −vh j+1))H ≤εkcn−chk nk2 H+1 4εkcn−vh nk2 H+1 2kc0−ch 0k2 H+1 2kc1−vh 1k2 H + n−1 X j=1 kkcj−chk jk2 H+ n−1 X j=1 1 4kkcj−vh j−(cj+1 −vh j+1)k2 H,(1.39) where ε > 0 is a positive parameter assumed to be small enough. Then, using the fact that a0=kc0−ch 0k2 H, we get, for all {vh j}n j=1 ⊂Vh, n X j=1 dj(vh j)≤ε an+1 4ε+1 2M+1 2a0 + n−1 X j=1 kkcj−chk jk2 H+ n−1 X j=1 1 4kkcj−vh j−(cj+1 −vh j+1)k2 H, and thus, estimate (1.38) can be written as follows, (1 −2ε)an≤2a0+T M +k n X j=1 β aj+1 2εM+M +2 n−1 X j=1 kkcj−chk jk2 H+1 2 n−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H,∀{vh j}n j=1 ⊂Vh, 1.4. Fully discrete approximations: numerical analysis 31 and finally an≤` gn+` k n X j=1 aj, n = 1,2, . . . , N, where `is a positive constant and gn:= a0+M+ n−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H. Applying a discrete version of Gronwall’s inequality with ` k ≤1/2 (see [18]), we find that max 0≤n≤Nan≤(`(1 + ` T e2`T )) max 0≤n≤Ngn. Therefore, we have proved the following result. Theorem 1.3 Under the assumptions of Theorem 1.1 and assuming that regularity condition (1.23) holds, there exists a positive constant β > 0, independent of the discretization parameters hand k, such that the following error estimates are satisfied for all {vh n}N n=1 ⊂Vh, max 0≤n≤Nkcn−chk nk2 H+k N X j=0 kcj−chk jk2 V+|γ0(cj−chk j)|2 ≤βhkc0−ch 0k2 H+ max 1≤n≤N  ∂c ∂t(tn)−δcn   2 H+kcn−vh nk2 V+I2 n + N−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 Hi.(1.40) Estimates (1.40) are the basis for the analysis of the convergence order. From now on and in order to approximate the space V, we consider the finite element space Vh defined in the following form: Vh={vh∈ C([0, l]) ; vh |[xi−1,xi]∈P1([xi−1, xi]),for i= 1,..., ˆ M, vh(l)=0}, (1.41) where the spatial discretization of the interval [0, l] is given by 0 = x0< x1< . . . < 38 Mixed kinetic-diffusion model for the Henry isotherm 1.5.2 Second example: simulation of hexanol As a second problem, we consider a dilute solution of the commercial alcohol hexanol, using the data from references [7] and [44], namely: cb= 3.44 mol/m3, D = 7.16 ×10−10m2/s, l = 10−4m, T= 0.5 s,Γ0= 0 mol/m2. Moreover, the initial condition c0is here defined as c0(x) = cbfor all x∈[0,10−4]. Using the discretization parameters h= 10−8m and k= 10−4s and the adsorption and desorption constants, ka H= 1.73 ×10−4m/s and kd H= 157 s−1, the concentration at final time and the evolution in time of the subsurface concentration are shown in Figure 1.2. Figure 1.2: Concentration at final time (left) and evolution in time of subsurface concentration (MATLAB results). Now, these results are compared to those obtained by using the commercial code COMSOL Multiphysics. Indeed, in Figure 1.3 the concentration at final time and the evolution in time of the subsurface concentration are plotted again. As it can be seen, these results are in good agreement with those obtained with our algorithm. Moreover, in Figures 1.4 and 1.5 we compare the evolution in time of both subsurface and surface concentrations, respectively, obtained with the linear mixed kineticdiffusion model described in this chapter, with that results obtained with the diffusion- 1.5. Numerical results 39 Figure 1.3: Concentration at final time (left) and evolution in time of subsurface concentration (COMSOL results). controlled model for the classical Henry’s isotherm, where the Henry equilibrium adsorption constant KHequals ka H/kd H. As it can be seen in Figures 1.4 and 1.5, the diffusion-controlled model predicts a faster equilibration of both subsurface and surface concentrations than the mixed kinetic-diffusion one. In the case of the latter model, the adsorption-desorption dynamics limits the mass transfer from the bulk solution to the surface due to the existence of an adsorption barrier that surfactant molecules have to undergo in order to move from the subsurface to the surface and viceversa. However, in the diffusion-controlled model, diffusion mechanics limits the entire process since the equilibration between the subsurface and surface layers is assumed to be immediate. Next, we analyze the dependence on the adsorption and desorption rate constants, ka Hand kd H, so we choose different values of these constants as reported in [44], leading to the following six cases: •Case i:ka H= 2.583 ×10−3m/s and kd H= 2348 s−1. •Case ii:ka H= 6.456 ×10−4m/s and kd H= 587 s−1. •Case iii:ka H= 1.73 ×10−4m/s and kd H= 157 s−1. •Case iv:ka H= 1.96 ×10−5m/s and kd H= 18 s−1. •Case v:ka H= 0 m/s and kd H= 0 s−1. 40 Mixed kinetic-diffusion model for the Henry isotherm 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 Time, s Subsurface concentration, mol/m3 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 Time, s Subsurface concentration, mol/m3 Figure 1.4: Evolution in time of the subsurface concentration with the mixed kinetic model (left) and that obtained with the diffusion-controlled model for Henry’s isotherm (right), semi-log scale. 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 4x 10−6 Time, s Surface concentration, mol/m2 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 4x 10−6 Time, s Surface concentration, mol/m2 Figure 1.5: Evolution in time of the surface concentration Γ(t) with the mixed kinetic model (left) and that obtained with the diffusion-controlled model using Henry’s isotherm (right), semi-log scale. •Case vi: diffusion-controlled model with Henry’s isotherm, KH= 1.1×10−6m. Our aim is to compare the surface tension eγgiven by eγ(t) = eγ0−n R θ Γ(t), for each of the above cases, where eγ0= 0.072 N/m denotes the surface tension of pure water, θ= 293.71 K is the temperature, R= 8.31 J/(K mol) represents the gas constant and nis a constant which is equal to one for a non-ionic surfactant. Using 1.5. Numerical results 41 the discretization parameters h= 10−8m and k= 10−5s for cases ii-vi and k= 10−6s for case i, in Figure 1.6 the evolution in time of the surface tension obtained for each of the above six cases is represented (semi-log scale). We point out that these numerical calculations are in good agreement with the experimental and theoretical values of the surface tensions of the hexanol solution reported in Figure 6 of [44] and Figure 27 of [7]. As it can be expected, the time needed to reach the stationary value decreases meanwhile the value of the adsorption rate constant, ka H, increases. This is because, if ka Hincreases, then the incorporation of surfactant molecules at the surface becomes faster and, the increasing of the surface concentration is closely related to the decreasing of the surface tension. Furthermore, as the value of the adsorption rate constant increases, the behavior predicted by the mixed kinetic-diffusion model approaches to the behavior predicted by the diffusion-controlled one. 10−4 10−3 10−2 10−1 100 0.062 0.064 0.066 0.068 0.07 0.072 Time, s Surface tension, N/m case i case ii case iii case iv case v case vi Figure 1.6: Surface tension graphs obtained for the six cases of adsorption and desorption constants, semi-log scale. Finally, the evolution in time of the surface concentration is shown in Figure 1.7 for the above six cases using the same discretization parameters utilized to obtain Figure 42 Mixed kinetic-diffusion model for the Henry isotherm 1.6. As we can observe in Figure 1.7, the surface concentration reaches its equilibrium value faster in the diffusion-controlled model than in the mixed kinetic-diffusion one. Moreover, as the adsorption rate constant ka Hdecreases, the adsorption process becomes slower and more time is needed to achieve the saturation at the surface. 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 4x 10−6 Time, s Surface concentration, mol/m2 case i case ii case iii case iv case v case vi Figure 1.7: Evolution in time of the surface concentration obtained for the six cases of adsorption and desorption constants, semi-log scale. 1.5.3 Third example: simulation of heptanol As a third example, we consider now a dilute solution of the commercial alcohol heptanol (see [44] for further details): cb= 0.1 mol/m3, D = 6.5×10−10m2/s, ka H= 7.04 ×10−4m/s, kd H= 190.27 s−1, l = 10−6m, T = 1 s,Γ0= 0 mol/m2. Moreover, the initial condition c0is defined as c0(x) = cbfor all x∈[0,10−6]. Using the discretization parameters h= 10−8m and k= 10−4s, the evolution in time 1.5. Numerical results 43 of the subsurface and the surface concentrations are shown in Figure 1.8 (left-hand side and right-hand side, respectively). We note that the subsurface concentration evolves to the constant bulk concentration cbin a fast way. 10−4 10−3 10−2 10−1 100 0.05 0.055 0.06 0.065 0.07 0.075 0.08 0.085 0.09 0.095 0.1 Time, s Subsurface concentration, mol/m3 10−4 10−3 10−2 10−1 100 0 0.5 1 1.5 2 2.5 3 3.5 4x 10−7 Time, s Surface concentration, mol/m2 Figure 1.8: Evolution in time of the subsurface concentration (left) and the surface concentration (right), semi-log scale. Finally, in Figure 1.9 we plot the evolution in time of the surface tension for several bulk concentrations (cb= 0.1 mol/m3,cb= 0.5 mol/m3and cb= 0.9 mol/m3). 10−4 10−3 10−2 10−1 100 0.063 0.064 0.065 0.066 0.067 0.068 0.069 0.07 0.071 0.072 0.073 Time, s Surface tension, N/m cb=0.1 cb=0.5 cb=0.9 Figure 1.9: Evolution in time of the surface tension for several heptanol bulk concentrations, semi-log scale. As it can be observed, the time needed to reach stationary values of the surface ten- 44 Mixed kinetic-diffusion model for the Henry isotherm sion depends both on the values for the adsorption rate constants (see Figure 1.6) and on the bulk concentration (see Figure 1.9). Moreover, increasing the bulk concentration, the equilibrium value of the surface tension decreases and so, the fall of the surface tension curve is greater as the concentration increases. Besides, as the bulk concentration decreases, the number of molecules in the solution becomes smaller. Therefore, the quantity of molecules achieving the surface decreases which implies that the rate of adsorption also decreases (see equation (1.6)) and consequently, the time needed for reaching the equilibrium surface tension increases. Chapter 2 Diffusion-controlled model with Langmuir isotherm In this chapter, we focus on the problem of modeling the surfactant behavior at the air-water interface considering a diffusion-controlled model. As it was said in the introduction of this manuscript, in this family of models, diffusion is the mechanism that governs the process since adsorption is assumed to be instantaneous. The adsorption dynamics is described here by the Langmuir isotherm, which has been used in a huge amount of literature (see, for example, [6, 7, 35]). This expression states a nonlinear relationship between the surface and subsurface concentrations and it is based on a lattice-type model (see [7, 16]) which assumes that the adsorption places on the lattice are equivalent, the probability of adsorption of the monomers at one empty space is independent of the occupied sites in its neighborhood and neither interactions nor intermolecular forces between the monomers in the lattice are considered. Mathematically, in this chapter, we deal with a non-standard parabolic problem. The reason why this problem is non-standard is because the boundary condition at the subsurface is coupled with the Langmuir isotherm, which makes the system to be nonlinear. For this problem, we prove the existence of weak solution by using the Rothe’s method, an intermediate problem (for which the existence of a unique weak solution is obtained applying Brouwer’s fixed-point theorem), a priori estimates and passing to 45 46 Diffusion-controlled model with Langmuir isotherm the limit. The uniqueness issue is solved using some arguments already introduced in [25], as the integration in time of the respective weak equations and the definition of adequate test functions. Moreover, a semi-discrete problem in time associated to an equivalent formulation of the weak problem is analyzed, proving some a priori estimates from which the linear convergence is achieved under additional regularity conditions. Then, fully discrete approximations, obtained by using the finite element method for the spatial discretization and a hybrid combination of both backward and forward Euler schemes, are presented. An error estimate result is proved from which the linear convergence is deduced under suitable regularity conditions. Finally, some numerical examples are shown to demonstrate the accuracy of this algorithm and the behavior of two commercially available surfactants. 2.1 The mathematical model Denoting by ˜c(t, x) the concentration of surfactant at time t∈[0, T] and point x∈[0, l] and by Γ(t) the time-dependent surface concentration, as it is usual, and taking into account the Fick’s law, we consider the diffusion partial differential equation: ∂˜c ∂t(t, x)−D∂2˜c ∂x2(t, x) = 0, t > 0, x ∈(0, l),(2.1) together with the boundary conditions (see [7, 35]): D∂˜c ∂x(t, 0) = dΓ dt (t), t > 0,(2.2) ˜c(t, l) = cb, t > 0,(2.3) and the initial conditions: ˜c(0, x) = ˜c0(x), x ∈(0, l),(2.4) Γ(0) = Γ0.(2.5) 2.1. The mathematical model 47 In equation (2.4), ˜c0(x) is a function defined in [0, l] which equals cbon x=l. We remind that the time-dependent surface concentration, Γ(t), is also an unknown of the system, so an additional condition is needed in order to close the problem. As we said previously, in this chapter, we consider the well-known and classical Langmuir isotherm (see [7]): Γ(t)=Γm KL˜c(t, 0) 1 + KL˜c(t, 0), t ≥0,(2.6) where Γmis the maximum surface concentration and KLis the Langmuir equilibrium adsorption constant. Here, we are interested in surfactant solutions below their cmc (critical micelle concentration), that is to say, we are interested in single-molecule transport (see [5]), therefore, the parameter Γm, which is a theoretical limit, cannot be reached. Moreover, it is usual in chemistry literature (see [7, 35]) to approximate the Langmuir isotherm by the Henry isotherm when the concentration is low or when KLc(t)<< 1. Then KH= ΓmKL. The Langmuir isotherm (2.6) was first deduced using kinetic arguments (see [7, 16]), by assuming that the rate of change of the surface concentration due to adsorption is equal to the rate of change of the surface concentration due to desorption. However, the Langmuir isotherm can also be deduced by molecular thermodynamic arguments for ideal non-localized adsorption. For the sake of clarity in the presentation of this chapter and without lost of generality, hereinafter we assume that the constants D, KLand Γmare equal to 1 and we define the nondecreasing Lipschitz function F:R→Ras follows F(z) =    z 1 + zif z≥0, 0 if z < 0. (2.7) Notice that a primitive to Fgiven by H(z) =    z−ln(1 + z) if z≥0, 0 if z < 0, (2.8) 54 Diffusion-controlled model with Langmuir isotherm We subtract the resulting two equations obtained for cs=c1 sand cs=c2 s, respectively, and take c1 s−c2 s∈Vas a test function, then Zl 0 (c1 s−c2 s)2dx +τZl 0∂(c1 s−c2 s) ∂x 2 dx +(F(γ0(c1 s) + cb)−F(γ0(c2 s) + cb))γ0(c1 s−c2 s)=0.(2.24) Since Fis nondecreasing, all terms in the left-hand side are nonnegative. Therefore, we can conclude from (2.24) that all its terms are equal to zero, and then c1 s=c2 sfor x∈(0, l). In order to prove (2.21), we take v=c+ s= max{cs,0} ∈ Vas a test function in (2.20) to get Zl 0 (c+ s)2dx+(F(γ0(cs)+cb)−F(γ0(cs−1)+cb))γ0(c+ s)+τZl 0∂c+ s ∂x 2 dx =Zl 0 cs−1c+ sdx. Notice that, if γ0(c+ s) = 0, then the second term of the previous equation disappears. On the contrary, if γ0(c+ s) is positive then γ0(cs) is positive. Moreover, since cs−1≤0 a.e. in (0, l) and cs−1∈V⊂ C([0, l]) (see [38]), it follows that γ0(cs−1)≤0. Then, due to the nondecreasing behavior of function Fwe know that F(γ0(cs)+cb)−F(γ0(cs−1)+ cb)≥0. Therefore, in both cases, the left-hand side of the previous equality is nonnegative, while the right-hand side is nonpositive and we can conclude that c+ s= 0 a.e. in (0, l). Thus cs≤0 a.e. in (0, l). Finally, we take v= (cs+C)−= max{0,−(cs+C)} ∈ H1(0, l). Notice that v(l) = max{0,−(cs(l) + C)}= max{0,−C}= 0, then v∈Vand it can be taken as a test function in equation (2.20) to obtain Zl 0 (cs−cs−1)(cs+C)−dx + (F(γ0(cs) + cb)−F(γ0(cs−1) + cb))γ0(cs+C)− −τZl 0∂(cs+C)− ∂x 2 dx = 0.(2.25) 2.3. Existence and uniqueness results 55 By using the hypothesis −C≤cs−1a.e. x∈(0, l) we have Zl 0 (cs−cs−1)(cs+C)−dx =Z[cs≤−C] (cs−cs−1)(cs+C)−dx ≤0. Moreover, if γ0(cs)<−C, then γ0(cs+C)−>0 and γ0(cs)< γ0(cs−1). Taking into account that Fis nondecreasing we get F(γ0(cs) + cb)≤F(γ0(cs−1) + cb). Hence, all terms in equation (2.25) are nonpositive and then (cs+C)−= 0 a.e. in (0, l) and, consequently, −C≤csa.e. in (0, l).  Now, regarding csas the solution to problem (2.20) in time t=swe define the following both piecewise constant and piecewise linear in time functions. Definition 2.1 Assuming that c0∈V, let csbe the solution to problem (2.20) at time t=s, s ∈N. Then, for (0, T] = SK s=1((s−1)τ, sτ], with τ=T/K and K∈N, we define a piecewise linear and a piecewise constant in time functions: ˜cτ, cτ: [0, T]→V by ˜cτ(t, x) := cs(x),(2.26) cτ(t, x) := (s−t τ)cs−1(x)+(t τ−s+ 1) cs(x),(2.27) for x∈(0, l)and (s−1)τ≤t < sτ,s= 1, . . . , K. Moreover, we define Fτ: [0, T]→R as follows Fτ(t) := (s−t τ)F(γ0(cs−1) + cb)+(t τ−s+ 1) F(γ0(cs) + cb),(2.28) for (s−1)τ≤t < sτ,s= 1, . . . , K. Remark 2.1 Note that ∂cτ ∂t (t, x) = cs(x)−cs−1(x) τ,(2.29) 56 Diffusion-controlled model with Langmuir isotherm dFτ dt (t) = F(γ0(cs) + cb)−F(γ0(cs−1) + cb) τ,(2.30) for x∈(0, l)and (s−1)τ < t < sτ,s= 1, . . . , K, and problem (2.20) can be written for a.e. t∈(0, T)as follows: Zl 0 ∂cτ ∂t v dx +d Fτ dt γ0(v) + Zl 0 ∂˜cτ ∂x ∂v ∂x dx = 0,∀v∈V. (2.31) Note also that cτ−˜cτ= ( t τ−s) (cs−cs−1)=(t τ−s)τ∂cτ ∂t ,(2.32) for x∈(0, l)and (s−1)τ < t < sτ,s= 1, . . . , K. Definition 2.2 Regarding the functions Fand Hdefined in (2.7) and (2.8), respectively, for s= 1, . . . , K, we define Ms:= Zl 0 c2 s 2dx +F(γ0(cs) + cb)(γ0(cs) + cb)−H(γ0(cs) + cb), and Ns:= cbF(γ0(cs) + cb). We have the following energy decay property. Lemma 2.3 Assuming that the hypothesis (H1) holds with C=cb, it follows that Ms+τZl 0∂cs ∂x 2 dx ≤Ms−1+cb, s = 1, . . . , K. (2.33) MK−NK≤ ··· ≤ Ms−Ns≤Ms−1−Ns−1≤ ··· ≤ M0−N0,(2.34) where cs∈V,s= 1, . . . , K, are the solutions to problem (2.20). Moreover, s X n=1 τZl 0∂cn ∂x 2 dx ≤M0+cb, s = 1, . . . , K, (2.35) 2.3. Existence and uniqueness results 57 and, Ms≤M0+cb, s = 1, . . . , K. (2.36) Proof. Taking v=csas test function in problem (2.20), we get, for s= 1, . . . , K, Zl 0 (cs−cs−1)csdx +F(γ0(cs) + cb)−F(γ0(cs−1) + cb)γ0(cs) + τZl 0∂cs ∂x 2 dx = 0. Furthermore, using the fact that x(x−y)≥(x2−y2)/2, for x, y ∈R, in the first term of the latter expression, we have, for s= 1, . . . , K, Zl 0 c2 s 2dx −Zl 0 c2 s−1 2dx +F(γ0(cs) + cb)−F(γ0(cs−1) + cb)(γ0(cs) + cb−cb) +τZl 0∂cs ∂x 2 dx ≤0.(2.37) Keeping in mind that (F(γ0(cs) + cb)−F(γ0(cs−1) + cb))(γ0(cs) + cb) = F(γ0(cs) + cb)(γ0(cs) + cb) −F(γ0(cs−1) + cb)(γ0(cs−1) + cb) +F(γ0(cs−1) + cb)(γ0(cs−1) + cb)−(γ0(cs) + cb),(2.38) and, since the primitive Hof F, defined in (2.8), is convex, we get (see [17]) H(γ0(cs−1) + cb)−H(γ0(cs) + cb)≤F(γ0(cs−1) + cb)γ0(cs−1)−γ0(cs).(2.39) Taking into account (2.38) and (2.39) in (2.37), we obtain, for s= 1, . . . , K, Zl 0 c2 s 2dx −Zl 0 c2 s−1 2dx +F(γ0(cs) + cb)(γ0(cs) + cb)−F(γ0(cs−1) + cb)(γ0(cs−1) + cb) +H(γ0(cs−1) + cb)−H(γ0(cs) + cb)−cbF(γ0(cs) + cb) + cbF(γ0(cs−1) + cb) +τZl 0∂cs ∂x 2 dx ≤0. 58 Diffusion-controlled model with Langmuir isotherm Therefore, it follows that, for s= 1, . . . , K, Ms−Ms−1−Ns+Ns−1+τZl 0∂cs ∂x 2 dx ≤0,(2.40) and we find that, for s= 1, . . . , K, Ms+τZl 0∂cs ∂x 2 dx ≤Ms−1−Ns−1+Ns.(2.41) We remark here that, since hypothesis (H1) holds with C=cb, we have −cb≤ γ0(cs−1)≤0 and then 0 ≤γ0(cs−1) + cb≤cb. Therefore, 0 ≤Ns−1≤cb. Analogously and applying Lemma 2.2 we get 0 ≤Ns≤cband, from (2.41), we conclude that Ms+τZl 0∂cs ∂x 2 dx ≤Ms−1+Ns≤Ms−1+cb, s = 1, . . . , K, (2.42) and then (2.33) holds. Moreover, from (2.40) and taking into account that its fifth term is nonnegative, we get Ms−Ns≤Ms−1−Ns−1, s = 1, . . . , K, and (2.34) holds. Also, from (2.40) we have Ms−Ns+τZl 0∂cs ∂x 2 dx ≤Ms−1−Ns−1, s = 1, . . . , K, and adding the term τ s−1 X n=1 Zl 0∂cn ∂x 2 dx in both sides of the latter inequality it follows that Ms−Ns+τ s X n=1 Zl 0∂cn ∂x 2 dx ≤M0−N0, s = 1, . . . , K. Finally, considering that Ns∈[0, cb], s = 0, . . . , K, we obtain, for s= 1, . . . , K, 2.3. Existence and uniqueness results 59 Ms+τ s X n=1 Zl 0∂cn ∂x 2 dx ≤M0−N0+Ns≤M0+Ns≤M0+cb.(2.43) Note that we can guarantee that Ms≥0 taking into account that its first term is nonnegative and using (2.17). Thus, from (2.43) we obtain (2.35) and (2.36).  We have the following a priori error estimates. Proposition 2.1 Assuming the hypothesis (H1) with C=cb, then functions ˜cτand cτ, defined in (2.26) and (2.27), respectively, are bounded in the space L2(0, T;H1(0, l)). Moreover, cτis bounded in H1(0, T;H)and Fτ, defined in (2.28), is bounded in H1(0, l) independently of τ. Furthermore, kcτ−˜cτk2 L2(0,T;H)≤C1τ2,(2.44) kγ0(cτ)−γ0(˜cτ)k2 L2(0,T)≤C2τ2,(2.45) where C1and C2are real, positive constants independent of τ. Proof. First, we prove that ˜cτis bounded in L2(0, T;H1(0, l)). Indeed, by definition we have k˜cτk2 L2(0,T;H)=ZT 0k˜cτ(t)k2 Hdt = K X s=1 Zsτ (s−1)τZl 0 (˜cτ(t, x))2dxdt = K X s=1 Zsτ (s−1)τZl 0 (cs(x))2dxdt. (2.46) Using property (2.17) and Lemma 2.3, we know that Zl 0 (cs(x))2 2dx ≤Ms≤M0+cb, s = 1, . . . , K, 60 Diffusion-controlled model with Langmuir isotherm and thus, Zl 0 (cs(x))2dx ≤2(M0+cb), s = 1, . . . , K. (2.47) Keeping in mind (2.46) and (2.47), it follows that k˜cτk2 L2(0,T;H)≤ K X s=1 Zsτ (s−1)τ (2M0+ 2cb)dt = (2M0+ 2cb)τK = (2M0+ 2cb)T. (2.48) Moreover, considering inequalities (2.35), (2.47) and (2.48), we have k˜cτk2 L2(0,T;H1(0,l)) =ZT 0k˜cτ(t)k2 H1(0,l)dt = K X s=1 Zsτ (s−1)τ Zl 0 (cs(x))2dx +Zl 0∂cs ∂x (x)2 dx!dt ≤(2M0+ 2cb)T+ K X s=1 τZl 0∂cs ∂x (x)2 dx ≤(2M0+ 2cb)T+M0+cb. Thus, we can conclude that ˜cτis bounded in L2(0, T;H1(0, l)) independently of τ. The following step is to show that cτis bounded in L2(0, T;H1(0, l)) as well. Indeed, by definition we get kcτk2 L2(0,T;H)=ZT 0kcτ(t)k2 Hdt =ZT 0k(s−t τ)cs−1+ ( t τ−s+ 1)csk2 Hdt. Regarding that f(x) = kxk2is a convex function and for (s−1)τ≤t≤sτ, s = 1, . . . , K, we get that 0 ≤s−t τ<1, then kcτk2 L2(0,T;H)≤ZT 0(s−t τ)kcs−1k2 H+ ( t τ−s+ 1)kcsk2 Hdt =ZT 0(s−t τ)Zl 0 (cs−1(x))2dx + ( t τ−s+ 1) Zl 0 (cs(x))2dxdt. Now, using inequality (2.47), we have 2.3. Existence and uniqueness results 61 kcτk2 L2(0,T;H)≤ K X s=1 Zsτ (s−1)τ(s−t τ)2(M0+cb)+(t τ−s+ 1)2(M0+cb)dt = K X s=1 Zsτ (s−1)τ 2(M0+cb)dt = 2(M0+cb) K X s=1 τ= 2(M0+cb)T. Using the same arguments, we also get     ∂cτ ∂x     2 L2(0,T;H) =ZT 0    ∂cτ ∂x (t)    2 H =ZT 0   (s−t τ)∂cs−1 ∂x + ( t τ−s+ 1)∂cs ∂x     2 H dt ≤ K X s=1 Zsτ (s−1)τ (s−t τ)Zl 0∂cs−1 ∂x (x)2 dx + ( t τ−s+ 1) Zl 0∂cs ∂x (x)2 dx!dt =τ 2Zl 0∂c0 ∂x (x)2 dx + K−1 X s=1 τ 2Zl 0∂cs ∂x (x)2 dx + K X s=1 τ 2Zl 0∂cs ∂x (x)2 dx, and, using (2.35) and keeping in mind that τ≤T, we obtain     ∂cτ ∂x     2 L2(0,T;H)≤τ 2kc0k2 V+1 2(M0+cb) + 1 2(M0+cb)≤T 2kc0k2 V+M0+cb. Now, in order to prove that cτis bounded in H1(0, T;H), it is enough to show that ∂cτ ∂t is bounded in L2(0, T;H) since the boundedness of cτin L2(0, T;H) has been already proved. Taking cs−cs−1∈Vas a test function in (2.31), we get, for a.e. t∈(0, T) and s= 1, . . . , K, Zl 0 ∂cτ ∂t (cs−cs−1)dx +d Fτ dt γ0(cs−cs−1) + Zl 0 ∂cs ∂x ∂(cs−cs−1) ∂x dx = 0. Then, considering (2.29) and (2.30), it follows that, for a.e. t∈(0, T) and s= 1, . . . , K, 62 Diffusion-controlled model with Langmuir isotherm Zl 0 (cs−cs−1)2 τdx +F(γ0(cs) + cb)−F(γ0(cs−1) + cb) τγ0(cs−cs−1) +Zl 0 ∂cs ∂x ∂(cs−cs−1) ∂x dx = 0. Using the fact that x(x−y)≥x2 2−y2 2, for x, y ∈R,in the third term of the previous equality, we have, for s= 1, . . . , K, Zl 0 (cs−cs−1)2 τdx +F(γ0(cs) + cb)−F(γ0(cs−1) + cb) τγ0(cs−cs−1) +Zl 0 1 2∂cs ∂x 2 dx ≤Zl 0 1 2∂cs−1 ∂x 2 dx. (2.49) Now, using (2.18), we obtain, for s= 1, . . . , K, Zl 0 (cs−cs−1)2 τdx +(F(γ0(cs) + cb)−F(γ0(cs−1) + cb))2 τ+Zl 0 1 2∂cs ∂x 2 dx ≤Zl 0 1 2∂cs−1 ∂x 2 dx. Adding the term s−1 X n=1 Zl 0 (cn−cn−1)2 τdx +(F(γ0(cn) + cb)−F(γ0(cn−1) + cb))2 τ, in both sides of the previous inequality, we find that, for s= 1, . . . , K, s X n=1 Zl 0 (cn−cn−1)2 τdx + s X n=1 F(γ0(cn) + cb)−F(γ0(cn−1) + cb)2 τ +Zl 0 1 2∂cs ∂x 2 dx ≤Zl 0 1 2∂c0 ∂x 2 dx. Then, since all terms of the left-hand side are nonnegative, it follows that, for s= 2.3. Existence and uniqueness results 63 1, . . . , K, s X n=1 Zl 0 (cn−cn−1)2 τdx ≤Zl 0 1 2∂c0 ∂x 2 dx, (2.50) s X n=1 (F(γ0(cn) + cb)−F(γ0(cn−1) + cb))2 τ≤Zl 0 1 2∂c0 ∂x 2 dx. (2.51) Therefore, using (2.29) and (2.50) we have     ∂cτ ∂t     2 L2(0,T;H) =ZT 0    ∂cτ ∂t (t)    2 H dt = K X s=1 Zsτ (s−1)τZl 0∂cτ ∂t (t, x)2 dx dt = K X s=1 Zsτ (s−1)τZl 0 (cs−cs−1)2 τ2dx dt = K X s=1 τZl 0 (cs−cs−1)2 τ2dx ≤Zl 0 1 2∂c0 ∂x 2 dx =kc0k2 V 2, and the result follows. Moreover, regarding Fτand keeping in mind (2.30), we obtain kFτk2 H1(0,T)=ZT 0|Fτ(t)|2dt +ZT 0 dFτ dt (t) 2 dt ≤ K X s=1 Zsτ (s−1)τ(s−t τ)|F(γ0(cs−1) + cb)|2+ ( t τ−s+ 1)|F(γ0(cs) + cb)|2dt + K X s=1 Zsτ (s−1)τ (F(γ0(cs) + cb)−F(γ0(cs−1) + cb))2 τ2dt. Taking into account that |F(z)| ≤ 1, for all z∈Rand applying (2.51), we get kFτk2 H1(0,T)≤ K X s=1 τ+ K X s=1 τ(F(γ0(cs) + cb)−F(γ0(cs−1) + cb))2 τ2 ≤T+Zl 0 1 2∂c0 ∂x 2 dx =T+kc0k2 V 2. Note also that 70 Diffusion-controlled model with Langmuir isotherm Note also that, using the integration by parts formula, we have Zt 0 dFτ dt (t)dt =Fτ(t)−Fτ(0) = Fτ(t)−F(γ0(c0) + cb),for a.e. t∈(0, T). Besides, using (2.63), passing to the limit in the previous expression and applying the integration by parts formula, we get, for a.e. t∈(0, T), F(γ0(c(t)) + cb)−F(γ0(c(0)) + cb) = Zt 0 dF(γ0(c(t)) + cb) dt dt =F(γ0(c(t)) + cb)−F(γ0(c0) + cb). The previous expression yields F(γ0(c(0)) + cb) = F(γ0(c0) + cb),(2.65) and, using (2.44) and (2.59), we deduce that ˜cτ→cin L2(QT). Then, for a subsequence it holds (see [3]) ˜cτ→ca.e. in QT.(2.66) By using hypothesis (H1), −cb≤cs(x)≤0 a.e. in (0, l) and then, by construction, −cb≤˜cτ≤0 also holds a.e. in QTand, keeping in mind (2.66), we get (2.55). Uniqueness. In order to prove the uniqueness of solution to Problem PL W, we proceed using several arguments already introduced in [25]. Anyway, for the sake of clarity of the presentation, we detail the main steps of the proof. Therefore, we consider ψ∈V and we define vτ,n(t, x) = ϕτ,n(t)ψ(x), 2.3. Existence and uniqueness results 71 where ϕτ,n(t) =              1 if t∈[0, τ], n(τ−t) + 1 if t∈[τ, τ +1 n], 0 if t∈[τ+1 n, T], for τ∈(0, T) and n∈N. Since vτ,n ∈ V, we can use it as a test function in equation (2.15) to get ZT 0h∂c ∂t(t), vτ,n(t)iV0×Vdt +ZT 0 ((c(t), vτ,n(t))) dt +ZT 0 d(F(γ0(c(t)) + cb)) dt γ0(vτ,n(t)) dt = 0.(2.67) Notice that vτ,n ∈H1(0, T;V) and therefore, using Theorem 11.5 in [8] and taking into account that vτ,n(T, x) = 0 for a.e. x∈(0, l), the first term of the previous expression reads, for all τ∈(0, T), ZT 0h∂c ∂t(t), vτ,n(t)iV0×Vdt =−ZT 0hc(t),∂vτ,n ∂t (t)iV×V0−(c(0), vτ,n(0))H =nZτ+1 n τ (c(t), ψ)Hdt −(c0, ψ)H,(2.68) Furthermore, using the integration by parts formula, considering ϕτ,n(T) = 0 in the third term of equations (2.67) and taking into account expression (2.65), we obtain, for all τ∈(0, T), ZT 0 d(F(γ0(c(t)) + cb)) dt γ0(vτ,n(t)) dt =−ZT 0 F(γ0(c(t)) + cb)γ0(ψ)d ϕτ,n dt (t)dt −F(γ0(c0) + cb)γ0(ψ)ϕτ,n(0) =Zτ+1 n τ n F(γ0(c(t)) + cb)γ0(ψ)dt −F(γ0(c0) + cb)γ0(ψ).(2.69) 72 Diffusion-controlled model with Langmuir isotherm Therefore, taking into account (2.68) and (2.69), equation (2.67) reads Zτ+1 n τ (c(t), ψ)Hn dt +ZT 0 ϕτ,n(t) ((c(t), ψ)) dt +Zτ+1 n τ n F(γ0(c(t)) + cb)γ0(ψ)dt = (c0, ψ)H+F(γ0(c0) + cb)γ0(ψ),∀τ∈(0, T).(2.70) Now, let c1and c2be two solutions to Problem PL W. Subtracting the resulting equations obtained from the previous expression for c=c1and c=c2, we get, for all τ∈(0, T), Zτ+1 n τ (c1(t)−c2(t), ψ)Hn dt +Zτ+1 n 0 ((c1(t)−c2(t), ψ))ϕτ,n(t)dt +Zτ+1 n τ (F(γ0(c1(t)) + cb)−F(γ0(c2(t)) + cb)) n γ0(ψ) = 0,∀ψ∈V. (2.71) Now, taking into account that c1, c2∈W2(0, T)⊂ C([0, T]; H) (see [39]) and using the mean value theorem, we have, for t?∈[τ, τ +1 n], Zτ+1 n τ (c1(t)−c2(t), ψ)Hn dt = (c1(t?)−c2(t?), ψ)H.(2.72) We notice that Zτ+1 n 0 ((c1(t)−c2(t), ψ)) ϕτ,n(t)dt =ZT 0 χ(0, τ +1 n) ((c1(t)−c2(t), ψ)) ϕτ,n(t)dt, (2.73) where χ(0, τ +1 n) denotes the characteristic function over the interval (0, τ +1 n). Now, we define a sequence of functions given by fn(t) := χ(0, τ +1 n) ((c1(t)−c2(t), ψ))ϕτ,n(t), n ∈N. We remark that fn∈L1(0, T) for each n∈N, and the family of functions fn, n ∈N satisfies that fn(t)−→ f(t),a.e. t∈(0, T), 2.3. Existence and uniqueness results 73 where f(t) = χ(0, τ)((c1(t)−c2(t), ψ)) and |fn(t)| ≤ g(t),a.e. t∈(0, T),(2.74) being g(t) = ((c1(t)−c2(t), ψ)). Then, applying the Lebesgue dominated convergence theorem, we can conclude that f∈L1(0, T) and Zτ+1 n 0 ((c1(t)−c2(t), ψ)) ϕτ,n(t)dt −→ Zτ 0 ((c1(t)−c2(t), ψ)) dt. (2.75) Moreover, considering that F(γ0(ci(t)) + cb)∈H1(0, T)⊂ C([0, T]),for i= 1,2, and using the mean value theorem, it follows that, for a given t?? ∈[τ, τ +1 n], Zτ+1 n τ (F(γ0(c1(t)) + cb)−F(γ0(c2(t)) + cb)) n ψ(0) dt = (F(γ0(c1(t??)) + cb)−F(γ0(c2(t??)) + cb)) n1 nψ(0).(2.76) Therefore, passing to the limit when n→ ∞ in (2.71) and taking into account (2.72), (2.75) and (2.76), it follows that (c1(τ)−c2(τ), ψ)H+Zτ 0 ((c1(t)−c2(t), ψ)) dt +(F(γ0(c1(τ)) + cb)−F(γ0(c2(τ)) + cb)) ψ(0) = 0,∀ψ∈V, a.e. τ∈(0, T).(2.77) Now, we fix τ∈(0, T) and we take ψ=c1(τ)−c2(τ) in (2.77) to obtain Zl 0 (c1(τ, x)−c2(τ, x))2dx +Zτ 0 ((c1(t)−c2(t), c1(τ)−c2(τ))) dt +(F(γ0(c1(τ)) + cb)−F(γ0(c2(τ)) + cb))(γ0(c1(τ)) −γ0(c2(τ))) = 0. 74 Diffusion-controlled model with Langmuir isotherm Since F is nondecreasing, the last term of the previous equality is nonnegative, and then, for a.e. τ∈(0, T), kc1(τ)−c2(τ)k2 H +Zτ 0Zl 0 (∂c1 ∂x (t, x)−∂c2 ∂x (t, x))(∂c1 ∂x (τ, x)−∂c2 ∂x (τ, x)) dx dt ≤0.(2.78) Taking into account ∂c1 ∂x −∂c2 ∂x ∈L2(0, T;H), we define the function β(τ) := Zτ 0 (∂c1 ∂x (s)−∂c2 ∂x (s))ds which belongs to W1,2(0, T;H) (see [40], page 104), being dβ dτ (τ) = ∂c1 ∂x (τ)−∂c2 ∂x (τ). Thus, we deduce (see Chapter III, Corollary 1.1, in [40]), 1 2 d dτ kβ(τ)k2 H= (dβ dτ (τ), β(τ))H =Zl 0Zτ 0 (∂c1 ∂x (s)−∂c2 ∂x (s))ds (∂c1 ∂x (τ)−∂c2 ∂x (τ))dx. Therefore, taking into account the Fubini Theorem (see Theorem IV.5 in [3]), we can change the order of the integrals and then replace the previous equality in estimate (2.78) to obtain, for a.e. τ∈(0, T), kc1(τ)−c2(τ)k2 H+1 2 d dτ kβ(τ)k2 H≤0. Integrating from 0 to T, we have ZT 0kc1(τ)−c2(τ)k2 Hdτ +1 2ZT 0 d dτ kβ(τ)k2 Hdτ ≤0, 2.4. Analysis of a semi-discrete problem 75 and therefore ZT 0kc1(τ)−c2(τ)k2 Hdτ +1 2kβ(T)k2 H≤0. Consequently, c1=c2a.e in QT. 2.4 Analysis of a semi-discrete problem In this section, we study the approximation in time of Problem PL W, proving some a priori estimates depending on the time discretization parameter. First, we rewrite Problem PL Win the following equivalent form, in terms of the derivative of function Ftaking into account that F(γ0(c)) ∈H1(0, T). Problem PL,eq W. For a given c0∈V, find a function c∈H1(0, T;H)∩L2(0, T;V) such that c(0) = c0and, for every v∈Vand a.e. t∈(0, T), we get (∂c ∂t(t), v)H+ ((c(t), v)) + h(γ0(c(t)) + cb,∂(γ0(c)) ∂t (t))γ0(v)=0,(2.79) where function h:R2→Ris defined as h(u, v) =      v (1 + u)2if u≥0, 0 elsewhere. Now, in order to obtain an approximation in time of Problem PL,eq W, we consider the uniform partition of the time interval [0, T] given in Section 1.4 of Chapter 1. Moreover, the same notation introduced there is also used hereinafter in this chapter. Therefore, applying a hybrid combination of both implicit and explicit Euler schemes, we get the following semi-discrete form of Problem PL,eq W. Problem Pk L. Find a sequence of functions ck={ck n}N n=0 ⊂Vsuch that ck n−ck n−1 k, vH + ((ck n, v)) + h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(v) = 0,∀v∈V, (2.80) 76 Diffusion-controlled model with Langmuir isotherm where ck 0=c0. After a straightforward application of Lax-Milgram theorem and proceeding as in the proof of Lemma 2.2, assuming hypothesis (H1) we can prove that Problem Pk Lhas a unique solution ck⊂Vsuch that −cb≤ck n(x)≤0 for a.e. xin (0, l), n = 0, . . . , N. (2.81) Hereafter, in this section, we will obtain some a priori estimates depending on the time discretization parameter assuming the following additional regularity condition: c∈ C1([0, T]; V).(2.82) Theorem 2.2 Assuming that hypothesis (H1) holds with C=cb. Let cand ckdenote the respective solutions to problems PL,eq Wand Pk L. Under the additional regularity condition (2.82), we have the following a priori error estimates: max 0≤n≤Nkcn−ck nk2 H+k N X j=1 kcj−ck jk2 V+ max 0≤n≤N|γ0(cn)−γ0(ck n)|2 ≤Ck N X j=1 n    ∂c ∂t(tj)−cj−cj−1 k    2 V +k2o, where Cdenotes a generic positive constant which may depend on the continuous solution cbut it is independent of the discretization parameter k, and whose value may change from line to line. Proof. Writing equation (2.79) at time t=tnand subtracting it to equation (2.80) we find that, for all v∈V, ∂c ∂t(tn)−ck n−ck n−1 k, vH + ((cn−ck n, v)) +h(γ0(cn) + cb,γ0(cn)−γ0(cn−1) k)−h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(v) +h(γ0(cn) + cb,∂(γ0(c)) ∂t (tn)−γ0(cn)−γ0(cn−1) k)γ0(v) = 0. 2.4. Analysis of a semi-discrete problem 77 Then, taking in the previous expression v=cn−ck nas a test function we have ∂c ∂t(tn)−ck n−ck n−1 k, cn−ck nH + ((cn−ck n, cn−ck n)) +h(γ0(cn) + cb,γ0(cn)−γ0(cn−1) k)−h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(cn−ck n) +h(γ0(cn) + cb,∂(γ0(c)) ∂t (tn)−γ0(cn)−γ0(cn−1) k)γ0(cn−ck n)=0. Taking into account that ∂c ∂t(tn)−ck n−ck n−1 k, cn−ck nH =∂c ∂t(tn)−cn−cn−1 k, cn−ck nH +cn−cn−1 k−ck n−ck n−1 k, cn−ck nH , cn−cn−1 k−ck n−ck n−1 k, cn−ck nH≥1 2kkcn−ck nk2 H−kcn−1−ck n−1k2 H, ((cn−ck n, cn−ck n)) = kcn−ck nk2 V, h(γ0(cn) + cb,∂(γ0(c)) ∂t (tn)−γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) ≤C    ∂c ∂t(tn)−cn−cn−1 k   Vkcn−ck nkV, where H¨older, Cauchy and trace inequalities and estimate (2.55) has been employed, we find that 1 2kkcn−ck nk2 H+kcn−ck nk2 V+hγ0(cn) + cb,γ0(cn)−γ0(cn−1) kγ0(cn−ck n) −hγ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) kγ0(cn−ck n) ≤1 2kkcn−1−ck n−1k2 H+C    ∂c ∂t(tn)−cn−cn−1 k   Vkcn−ck nkV −∂c ∂t(tn)−cn−cn−1 k, cn−ck nH . Using the Cauchy inequality with a small parameter, see (1.16), the fact that the norms k·kH1(0,l)and k·kVare equivalent and multiplying by 2k, it follows that 78 Diffusion-controlled model with Langmuir isotherm kcn−ck nk2 H+β kkcn−ck nk2 V+ 2k hγ0(cn) + cb,γ0(cn)−γ0(cn−1) kγ0(cn−ck n) −2k hγ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) kγ0(cn−ck n) ≤ kcn−1−ck n−1k2 H+Ck     ∂c ∂t(tn)−cn−cn−1 k    2 V +Ckkcn−ck nk2 H, where βis a positive constant, which is independent of the time discretization parameter k, but, as in the case of constant C, it may change from line to line. Now, keeping in mind that h(γ0(cn) + cb,γ0(cn)−γ0(cn−1) k)−h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(cn−ck n) =h(γ0(cn) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) −h(γ0(ck n−1) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) +h(γ0(ck n−1) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) −h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(cn−ck n), the Lipschitz behavior of function N(z) = 1 (1 + z)2for z∈[0, cb] and h(γ0(cn) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) −h(γ0(ck n−1) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) ≤Rn 1 (1 + γ0(cn) + cb)2−1 (1 + γ0(ck n−1) + cb)2|γ0(cn−ck n)| ≤CRn|γ0(cn)−γ0(ck n−1)||γ0(cn−ck n)| ≤CRn|γ0(cn)−γ0(cn−1)|+|γ0(cn−1)−γ0(ck n−1)||γ0(cn−ck n)| ≤CRn|γ0(cn)−γ0(ck n)|2+|γ0(cn−1)−γ0(ck n−1)|2+|γ0(cn)−γ0(cn−1)|2, h(γ0(ck n−1) + cb,γ0(cn)−γ0(cn−1) k)γ0(cn−ck n) −h(γ0(ck n−1) + cb,γ0(ck n)−γ0(ck n−1) k)γ0(cn−ck n) ≥β 2k|γ0(cn)−γ0(ck n)|2−|γ0(cn−1)−γ0(ck n−1)|2, 2.4. Analysis of a semi-discrete problem 79 where Cauchy’s inequality and estimate (2.81) have been used and Rnis the error Rn= γ0(cn)−γ0(cn−1) k≤CkckC1([0,T];V),(2.83) we get kcn−ck nk2 H+β kkcn−ck nk2 V+β|γ0(cn)−γ0(ck n)|2 ≤ kcn−1−ck n−1k2 H+Ck     ∂c ∂t(tn)−cn−cn−1 k    2 V +Ckkcn−ck nk2 H +β|γ0(cn−1)−γ0(ck n−1)|2+Ck|γ0(cn)−γ0(ck n)|2+Ck|γ0(cn−1)−γ0(ck n−1)|2 +Ck|γ0(cn)−γ0(cn−1)|2. Considering inequality (2.83) and using the following notation αL n:= kcn−ck nk2 H+β|γ0cn−γ0ck n|2, λL n:= βkcn−ck nk2 V, φL n:= C    ∂c ∂t(tn)−cn−cn−1 k    2 V +Ckcn−ck nk2 H+C k2+C|γ0(cn)−γ0(ck n)|2 +β|γ0(cn−1)−γ0(ck n−1)|2, we deduce that αL n+k λL n≤αL n−1+k φL n, and then, we get kcn−ck nk2 H+β k n X j=1 kcj−ck jk2 V+β|γ0(cn)−γ0(ck n)|2≤Ck n X j=1 nkcj−ck jk2 H +    ∂c ∂t(tj)−cj−cj−1 k    2 V +|γ0(cj)−γ0(ck j)|2+k2o, where the initial condition ck 0=c0has been used. Then, defining aL n:= kcn−ck nk2 H+β k n X j=1 kcj−ck jk2 V+β|γ0(cn)−γ0(ck n)|2, gL n:= k n X j=1 k2+    ∂c ∂t(tj)−cj−cj−1 k    2 V!, 86 Diffusion-controlled model with Langmuir isotherm in mind that 0 ≤RL(·)≤cb, we obtain the following estimates, for all v∈V, |(h(RL(cb+γ0(cn)),∂ γ0(c) ∂t (tn)) −h(RL(cb+γ0(chk n−1)),∂ γ0(c) ∂t (tn))) γ0(v)| ≤C|γ0(cn−chk n−1)||γ0(v)|,(2.97) |(h(RL(cb+γ0(chk n−1)),∂γ0(c) ∂t (tn)) −h(RL(cb+γ0(chk n−1)), γ0(δcn))) γ0(v)| ≤C|γ0(∂c ∂t(tn)−δcn)||γ0(v)|,(2.98) h(RL(cb+γ0(chk n−1)), γ0(δcn)) −h(RL(cb+γ0(chk n−1)), γ0(δchk n))γ0(cn−chk n) ≥β 2k|γ0(cn−chk n)|2−|γ0(cn−1−chk n−1)|2,(2.99) where βis a positive constant which is independent of the discretization parameters. We recall that δcn= (cn−cn−1)/k. Therefore, using (2.96), (2.97), (2.98) and (2.99), equation (2.95) leads to the following estimates 1 2kkcn−chk nk2 H+kcn−chk nk2 V+β 2k|γ0(cn−chk n)|2≤1 2kkcn−1−chk n−1k2 H +β 2k|γ0(cn−1−chk n−1)|2+∂c ∂t(tn)−δchk n, cn−vhH + ((cn−chk n, cn−vh)) +h(RL(cb+γ0(chk n−1)), γ0(δcn)) −h(RL(cb+γ0(chk n−1)), γ0(δchk n))γ0(cn−vh) +C|γ0(cn−chk n−1)||γ0(cn−chk n)|+Cγ0∂c ∂t(tn)−δcn|γ0(cn−chk n)| +C|γ0(cn−chk n−1)||γ0(cn−vh)|+Cγ0∂c ∂t(tn)−δcn|γ0(cn−vh)| +δcn−∂c ∂t(tn), cn−chk nH ,∀vh∈Vh.(2.100) Taking into account that ∂c ∂t(tn)−δchk n, cn−vhH =∂c ∂t(tn)−δcn, cn−vhH +δcn−δchk n, cn−vhH, |γ0(cn−chk n−1)|≤|γ0(cn−cn−1)|+|γ0(cn−1−chk n−1)|, 2.5. Fully discrete approximations: a priori error estimates 87 and using several times H¨older and both Cauchy and Cauchy with ε, see (1.16), inequalities, and considering the property (1.31) with ℘= 2, we have, for all vh∈Vh, 1 2kkcn−chk nk2 H+kcn−chk nk2 V+β 2k|γ0(cn−chk n)|2≤1 2kkcn−1−chk n−1k2 H +β 2k|γ0(cn−1−chk n−1)|2+C  ∂c ∂t(tn)−δcn   2 H+Ckcn−vhk2 V +h(RL(cb+γ0(chk n−1)), γ0(δcn)) −h(RL(cb+γ0(chk n−1)), γ0(δchk n))γ0(cn−vh) +(δcn−δchk n, cn−vh)H+εkcn−chk nk2 V+C|γ0(cn−cn−1)|2+C|γ0(cn−1−chk n−1)|2 +C|γ0(cn−chk n)|2+Cγ0∂c ∂t(tn)−δcn 2+C|γ0(cn−vh)|2+Ckcn−chk nk2 H. Thus, by induction we find that kcn−chk nk2 H+k n X j=1 kcj−chk jk2 V+|γ0(cn−chk n)|2≤Ckc0−ch 0k2 H +C|γ0(c0−ch 0)|2+Ck n X j=1 n  ∂c ∂t(tj)−δcj   2 H+kcj−vh jk2 V +h(RL(cb+γ0(chk j−1)), γ0(δcj)) −h(RL(cb+γ0(chk j−1)), γ0(δchk j))γ0(cj−vh j) +(δcj−δchk j, cj−vh j)H+|γ0(cj−cj−1)|2+|γ0(cj−chk j)|2+γ0∂c ∂t(tj)−δcj 2 +|γ0(cj−vh j)|2+kcj−chk jk2 Ho,∀vh={vh j} ⊂ Vh.(2.101) Now, denoting by sj=1 (1 + RL(cb+γ0chk j−1))2, keeping in mind that functions N(z) = 1 (1 + z)2for z∈[0, cb] and RLare Lipschitz and using both (2.87) and (2.90), we deduce that, for all vh={vh j} ⊂ Vh, k n X j=1 h(RL(cb+γ0(chk j−1)), γ0(δcj)) −h(RL(cb+γ0(chk j−1)), γ0(δchk j))γ0(cj−vh j) = n X j=1 γ0(cj−chk j)−γ0(cj−1−chk j−1)sjγ0(cj−vh j) = n−1 X j=1 γ0(cj−chk j)sjγ0(cj−vh j)−sj+1 γ0(cj+1 −vh j+1) +γ0(cn−chk n)γ0(cn−vh n)sn+γ0(ch 0−c0)γ0(c1−vh 1)s1, 88 Diffusion-controlled model with Langmuir isotherm |sjγ0(cj−vh j)−sj+1γ0(cj+1 −vh j+1)| ≤sjγ0(cj−vh j−(cj+1 −vh j+1))+|γ0(cj+1 −vh j+1)(sj−sj+1)|, |sj−sj+1| ≤ C|γ0(chk j−chk j−1)| ≤ C√kkch 0kV. Therefore, taking into account both Cauchy and Cauchy with ε > 0 (see (1.16)) inequalities, we notice that k n X j=1 h(RL(cb+γ0(chk j−1)), γ0(δcj)) −h(RL(cb+γ0(chk j−1)), γ0(δchk j))γ0(cj−vh j) ≤C n−1 X j=1 (k|γ0(cj−chk j)|2+1 k|γ0(cj−vh j−(cj+1 −vh j+1))|2+kch 0k2 V|γ0(cj+1 −vh j+1)|2) +C ε|γ0(cn−chk n)|2+C|γ0(cn−vh n)|2+C|γ0(c0−ch 0)|2+C|γ0(c1−vh 1)|2. Using the previous estimate and estimate (1.39), see Chapter 1, expression (2.101) reads kcn−chk nk2 H+k n X j=1 kcj−chk jk2 V+|γ0(cn−chk n)|2≤Ckc0−ch 0k2 H +C|γ0(c0−ch 0)|2+Ck n X j=1 n  ∂c ∂t(tj)−δcj   2 H+kcj−vh jk2 V+γ0∂c ∂t(tj)−δcj 2 +|γ0(cj−chk j)|2+|γ0(cj−vh j)|2+kcj−chk jk2 H+|γ0(cj−cj−1)|2o +Ckcn−vh nk2 H+Ckc1−vh 1k2 H+C|γ0(cn−vh n)|2+C|γ0(c1−vh 1)|2 +C k n−1 X j=1 kcj−vh j−(cj+1 −vh j+1)k2 H+|γ0(cj−vh j−(cj+1 −vh j+1))|2 +C n−1 X j=1 |γ0(cj+1 −vh j+1)|2,∀vh={vh j} ⊂ Vh. Now, defining bL n:= kcn−chk nk2 H+k n X j=1 kcj−chk jk2 V+|γ0(cn−chk n)|2, 2.5. Fully discrete approximations: a priori error estimates 89 dL n:= kc0−ch 0k2 H+|γ0(c0−ch 0)|2+k n X j=1 n  ∂c ∂t(tj)−δcj   2 H+kcj−vh jk2 V +γ0∂c ∂t(tj)−δcj 2+|γ0(cj−vh j)|2+|γ0(cj−cj−1)|2o+kcn−vh nk2 H +kc1−vh 1k2 H+|γ0(cn−vh n)|2+|γ0(c1−vh 1)|2+C n−1 X j=1 |γ0(cj+1 −vh j+1)|2 +1 k n−1 X j=1 kcj−vh j−(cj+1 −vh j+1)k2 H+|γ0(cj−vh j−(cj+1 −vh j+1))|2, it follows that bL n≤C k n X j=1 bL j+C dL n, n = 1, . . . , N. Finally, applying the discrete version of Gronwall’s inequality presented in Section 1.4 of Chapter 1 (see, for example, [18]), the result follows.  Remark 2.3 We note that estimates (2.92) could be also obtained without the estimates on function h, keeping in mind that h(RL(cb+γ0(chk n−1)), γ0(δcn−δchk n))γ0(cn−chk n) =hRL(cb+γ0(chk n−1)),1 2kk2γ0(δcn−δchk n)2+γ0(cn−chk n)2−γ0(cn−1−chk n−1)2 ≥C 2kk2γ0(δcn−δchk n)2+γ0(cn−chk n)2−γ0(cn−1−chk n−1)2, and the estimate h(RL(cb+γ0(chk n−1)), γ0(δcn−δchk n))γ0(cn−vh) ≤Cεkγ0(δcn−δchk n)2+C kγ0(cn−vh)2, where ε > 0is assumed small enough. Estimates (2.92) are the basis for the convergence analysis. As an example, recall that the finite element space Vhis given in (1.41), and let us assume further regularity conditions on the solution to the continuous problem: c∈H1(0, T;H2(0, l)),∂2c ∂t2∈L2(0, T;V).(2.102) 90 Diffusion-controlled model with Langmuir isotherm Denoting by πh:C([0, l]) →Vhthe standard finite element interpolation operator (see [9]) and considering ch 0=πhc0, we are able to prove the following. Corollary 2.2 Let the assumptions of Theorem 2.3 and the additional regularity conditions (2.102) hold. Then the linear convergence of the algorithm is obtained; i.e. there exists a positive constant C > 0, independent of hand k, such that max 0≤n≤Nkcn−chk nkH+ max 0≤n≤N|γ0(cn−chk n)| ≤ C(h+k). Proof. Let us take vh j=πhcj, j = 1, . . . , N. Since c∈ C([0, T]; H2(0, l)) because H1(0, T)⊂ C([0, T]) we obtain (see [9]), k N X n=1 kcn−πhcnk2 V+|γ0(cn−πhcn)|2+kc0−ch 0k2 H+|γ0(c0−ch 0)|2 + max 0≤n≤Nkcn−πhcnk2 H+ max 0≤n≤N|γ0(cn−πhcn)|2≤C h2kck2 C([0,T];H2(0,l)). Keeping in mind the regularity ∂2c ∂t2∈L2(0, T;V) and expression (2.85), considering the trace inequality, the fact that k·kVand k·kH1(0,l)are equivalent norms in Vand the inequality (2.85), we have k N X n=1   ∂c ∂t(tn)−δcn   2 H+γ0∂c ∂t(tn)−δcn) 2≤C k N X n=1   ∂c ∂t(tn)−δcn   2 V ≤Ck2  ∂2c ∂ t2   2 L2(0,T;V). Now, following the ideas applied to estimate the damage error terms (see, for instance, [4]) we bound the terms 1 k N−1 X n=0 kcn−vh n−(cn+1 −vh n+1)k2 H+|γ0(cn−vh n−(cn+1 −vh n+1))|2. First, note that both cnand cn+1 belong to H2(0, l) and then, taking into account the 2.6. Numerical results 91 linearity of the interpolation operator, we get (see [9]), kcn+1 −cn−πh(cn+1 −cn)k2 H+|γ0(cn+1 −cn−πh(cn+1 −cn))|2≤C h4kcn+1 −cnk2 H2(0,l). We point out now that the second term in the latter expression, as well as the last term in estimates (2.92), are zero taking into account that γ0(cn) = γ0(πhcn), n= 0, . . . , N, and the linearity of the trace operator. On the other hand, using regularity condition (2.102) we deduce that cn+1 −cn=Ztn+1 tn ∂c ∂t(s)ds. Thus, we have kcn+1 −cnkH2(0,l)≤Ztn+1 tn  ∂c ∂t(s)  H2(0,l)ds ≤√kZtn+1 tn  ∂c ∂t(s)   2 H2(0,l)ds1/2 , and therefore, keeping in mind that h≤lit follows that 1 k N−1 X n=0 kcn−πhcn−(cn+1 −πhcn+1)k2 H≤C h4 N−1 X j=1 Ztn+1 tn  ∂c ∂t(s)   2 H2(0,l)ds ≤C h2  ∂c ∂t   2 L2(0,T;H2(0,l)). Combining all these estimates, the linear convergence is obtained.  2.6 Numerical results In this section, we first describe the numerical scheme implemented in MATLAB in order to obtain the numerical approximations of Problem Phk Land then, we present some numerical results to exhibit its accuracy in an academic example and its behavior in the simulation of two commercially available surfactants. Considering the finite element space defined in (1.41), for n= 1,2, . . . , N and given chk n−1∈Vh, the discrete concentration at time t=tnof surfactant, chk n, is then obtained 92 Diffusion-controlled model with Langmuir isotherm from equation (2.88); namely, it solves the problem: (chk n, vh)H+k((chk n, vh)) + h(RL(cb+γ0(chk n−1)), γ0(chk n)) γ0(vh) = (chk n−1, vh)H+h(RL(cb+γ0(chk n−1)), γ0(chk n−1)) γ0(vh),∀vh∈Vh. The algorithm implemented to solve this problem is described below: 1. Initial time step. At the beginning both chk 0and Γ0are given. We calculate µhk 0=1 (1 + RL(cb+γ0(chk 0)))2. 2. (n)th time step. The surfactant concentration at time tn−1,chk n−1, and the value µhk n−1are known. Then, at time tn,chk n,µhk nand Γhk nare obtained using the following algorithm: (a) We calculate chk nby solving the following linear problem: Zl 0 chk nvhdx +kZl 0 ∂chk n ∂x ∂vh ∂x dx +µhk n−1γ0(chk n)γ0(vh) =µhk n−1γ0(chk n−1)γ0(vh) + Zl 0 chk n−1vhdx, ∀vh∈Vh. (b) Now, µhk nis obtained by using the formula: µhk n=1 (1 + RL(cb+γ0(chk n)))2, (c) and the value of Γhk nis easily deduced: Γhk n=γ0(chk n) 1 + RL(cb+γ0(chk n)). This numerical scheme has been implemented on a 3.2 Ghz PC using MATLAB, and a typical run (h=k= 0.01) takes about 0.6 seconds of CPU time. 2.6. Numerical results 93 2.6.1 First example: numerical convergence As a first example, we consider the following test problem: ∂˜c ∂t(t, x)−5∂2˜c ∂x2(t, x) = 0, x ∈(0,1), t ∈(0,0.1), 5∂˜c ∂x(t, 0) = h(˜c(t, 0),∂˜c ∂t(t, 0)), t ∈(0,0.1), ˜c(t, 1) = 1, t ∈(0,0.1), ˜c(0, x) = ˜c0(x), with the initial condition ˜c0(x) = min{1,1000 x}. This problem corresponds to problem (2.1), (2.3)-(2.4) and (2.9) with the following data: l= 1, T = 0.1, cb= 1, D = 5,Γm= 1, KL= 1,Γ0= 0. Taking the solution obtained with parameters h= 1/16384 and k= 10−6as the “exact solution”, c, the numerical errors, which are given by max 0≤n≤Nkcn−chk nkH+ max 0≤n≤N|γ0(cn−chk n)|, are presented in Table 2.1 for several values of the discretization parameters hand k. As it can be seen, the numerical error tends to zero as both hand kdo. Moreover, the graph of the error with respect to the parameter h+kis shown in Figure 2.2, where the linear convergence, stated in Corollary 2.2, seems to be achieved. 2.6.2 Second example: simulation of propanol As a second problem, we consider a solution of propanol, using the following data from reference [6], namely: cb= 333 mol/m3, D = 5.2×10−10 m2/s, KL= 5.5×10−3m3/mol, Γm= 7.1×10−6mol2/m2, l = 10−4m, T = 10−4s,Γ0= 0 mol/m2. 94 Diffusion-controlled model with Langmuir isotherm h↓k→0.01 0.005 0.002 0.001 0.0005 1/8 0.437554 0.365515 0.312296 0.291710 0.280590 1/16 0.324326 0.247134 0.188115 0.164276 0.150916 1/32 0.267746 0.187936 0.125745 0.099936 0.085071 1/64 0.239626 0.158547 0.094766 0.067912 0.052196 1/128 0.225634 0.143942 0.079381 0.052003 0.035847 1/256 0.218658 0.136667 0.071724 0.044087 0.027711 1/512 0.215176 0.133037 0.067905 0.040136 0.023654 1/1024 0.213518 0.131309 0.066087 0.038261 0.021724 1/2048 0.213497 0.131288 0.066065 0.038238 0.021699 1/4096 0.213487 0.131277 0.066054 0.038227 0.021688 Table 2.1: Numerical errors (×101) for several time and spatial discretization parameters. 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0.2 0.4 0.6 0.8 1 1.2 1.4 h+k Numerical errors (x 10) Asymptotic convergence Figure 2.2: Example 1: linear convergence. Moreover, the initial condition ˜c0is here defined as ˜c0(x) =      0 if x= 0, 333 if x∈(0,10−4]. Using the time discretization parameter k= 10−9s and a non-uniform spatial mesh, refined as we approach to the point x= 0 and with the smallest element length 10−11 m, the evolution in time of both surface and subsurface concentrations are shown in Figure 2.6. Numerical results 95 2.3. As it can be seen, the subsurface concentration tends to the bulk concentration as time evolves, while the surface concentration increases but it converges to a value below Γmand determined by the Langmuir isotherm (2.6). 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x 10−4 0 50 100 150 200 250 300 350 Time, s Subsurface concentration, mol/m 3 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x 10−4 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5x 10−6 Time, s Surface concentration, mol/m 2 Figure 2.3: Evolution in time of subsurface and surface concentrations, respectively. The surface equation of state, relating the surface tension eγwith the subsurface concentration c(t, 0), is given by eγ(t) = eγ0−n R θ Γmln(1 + KLc(t, 0)),(2.103) where we take eγ0= 0.0725 N/m, θ= 293 K and we recall that R= 8.31 J/(K mol) and n= 1. In Figure 2.4 the evolution in time of the surface tension obtained with our algorithm is compared with the results provided in [6]. We observe that both results are in good agreement. 2.6.3 Third example: simulation of sodium dodecylsulfate In this last example, we consider a solution of sodium dodecylsufate (SDS) and we use the following data, obtained from [6]: D= 0.1×10−10 m2/s, KL= 0.11 m3/mol, l = 10−4m, Γm= 10 ×10−6mol2/m2, T = 10 s,Γ0= 0 mol/m2. 102 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation zero and so, a homogeneous Dirichlet boundary condition is imposed on the right end of the spatial interval. We turn now to the variational formulation of problem (3.1)-(3.5) and (3.8). So, assuming regularity, multiplying equation (3.1) by a smooth function zdefined in [0, l] such that z(l) = 0; integrating in (0, l) and using the integration by parts formula, we obtain Zl 0 ∂c ∂t(t, x)z(x)dx +DZl 0 ∂c ∂x(t, x)∂z ∂x(x)dx +D∂c ∂x(t, 0) z(0) = 0, for a.e. t∈(0, T). Using the equations (3.2) and (3.8), we find Zl 0 ∂c ∂t(t, x)z(x)dx +DZl 0 ∂c ∂x(t, x)∂z ∂x(x)dx +f(c(t, 0),Γ(t)) z(0) = kd LΓ(t)z(0), for a.e. t∈(0, T). From the latter, using (3.4), (3.5) and (3.8) and taking into account the notations introduced in Chapter 1, we obtain the following weak formulation of the problem: Problem PLH W. For given c0∈Hand Γ0∈R, find c∈W2(0, T) and Γ ∈H1(0, T) such that h∂c ∂t(t), viV0×V+D((c(t), v)) + f(γ0(c(t)),Γ(t))γ0(v) = kd LΓ(t)γ0(v), for a.e. t∈(0, T),∀v∈V, dΓ dt (t) = f(γ0(c(t)),Γ(t)) −kd LΓ(t),for a.e. t∈(0, T), c(0) = c0,Γ(0) = Γ0. Besides, in the sequel, we need the truncation operator R:R→Rgiven by R(s) =              0 if s < 0, sif 0 ≤s≤(1 −σ)Γm, (1 −σ)Γmif s > (1 −σ)Γm, (3.9) 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 103 where 0 ≤σ < 1 is a given small constant. Indeed, we use the previously defined function fand the truncation operator Rto introduce the following truncated version of the Langmuir–Hinshelwood equation, which becomes a modification of equation (3.8) and reads dΓ dt (t) = f(c(t, 0), R(Γ(t))) −kd LΓ(t), t > 0.(3.10) We notice that, in spite of being a mathematical tool, the truncated operator Rhas also physical sense since it acts on the surface concentration, its real values being nonnegative but below Γm, since, as mentioned, we are interested in surfactant concentrations below their cmc. Similarly as we did previously, from (3.1)-(3.5) and (3.10), we arrive to the following truncated problem associated to Problem PLH W: Problem PLH R. For given c0∈Hand Γ0∈R, find c∈W2(0, T) and Γ ∈H1(0, T) such that h∂c ∂t(t), viV0×V+D((c(t), v)) + f(γ0(c(t)), R(Γ(t)))γ0(v) = kd LΓ(t)γ0(v), for a.e. t∈(0, T),∀v∈V, dΓ dt (t) = f(γ0(c(t)), R(Γ(t))) −kd LΓ(t),for a.e. t∈(0, T), c(0) = c0,Γ(0) = Γ0. We remark that the initial conditions in Problems PLH Wand PLH Rmake sense since W2(0, T)⊂ C([0, T]; H) and H1(0, T)⊂ C([0, T]). Here and in what follows, by C([0, T]) we denote the space of continuous functions from [0, T] to Rwith the maximum norm kvkC([0,T]) = max{|v(t)|;t∈[0, T]}. In the following two sections we study the existence and uniqueness of solution to Problems PLH Rand PLH Wand we also analyze the relation between these two problems and their solutions. 104 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation 3.1.2 Existence and uniqueness results for Problem PLH R In this section we formulate and prove the main existence and uniqueness result for Problem PLH R. Theorem 3.1 Assume that D,kd L,ka Land Γmare positive constants, and Γ0∈R, c0∈H. Then Problem PLH Rhas a unique solution (c, Γ) ∈W2(0, T)×H1(0, T). The proof of Theorem 3.1 is carried out in several steps and it is based on the study of two intermediate problems, followed by the application of the Schauder fixed-point theorem. So, before demonstrating this result, we introduce all the tools needed to prove it. To simplify the presentation in this section, and without loss of generality, we can suppose that D=kd L=ka L= Γm= 1, and therefore the nonlinear term in Problem PLH Ris of the form f(r, s) = r(1 −s) for r,s∈R. Let abe a given positive constant, which represents an arbitrary time. Intermediate parabolic problem. Let η∈ C([0, a]) and consider the following problem: Problem Pη 1. Given c0∈H, find cη∈W2(0, a) such that h∂cη ∂t (t), viV0×V+ ((cη(t), v)) + γ0(cη(t))(1 −R(η(t)))γ0(v) = η(t)γ0(v), for a.e. t∈(0, a),∀v∈V, cη(0) = c0. For Problem Pη 1we have the following result regarding that Ctr denotes the L(V, R)- norm of the trace operator, see (1.9). Lemma 3.1 Problem Pη 1has a unique solution cη∈W2(0, a). Moreover, we have kγ0(cη)kL2(0,a)≤C2 tr √akηkC([0,a]) +Ctr kc0kH.(3.11) 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 105 Proof. The existence and uniqueness of solution is based on the classical result on evolution problems. We define a family of bilinear forms aη: (0, a)×V×V→Rby aη(t;u, v) = ((u, v)) + γ0(u) (1 −R(η(t)))γ0(v), t ∈(0, a),∀u, v ∈V, and the function gη: (0, a)→V0by hgη(t), viV0×V=η(t)γ0(v), t ∈(0, a),∀v∈V. Under this notation, Problem Pη 1has the form h∂cη ∂t (t), viV0×V+aη(t;cη(t), v) = hgη(t), viV0×V,for a.e. t∈(0, a),∀v∈V, cη(0) = c0. Exploiting the properties that aη(·, u, v) is measurable for all u,v∈V,aη(t;·,·) is continuous on V×Vfor a.e. t∈(0, a) (i.e. |aη(t;u, v)| ≤ Mkukkvkfor all u,v∈V, a.e. t∈(0, a) with M > 0), aη(t;·,·) is coercive for a.e. t∈(0, a) (i.e. aη(t;u, u)≥αkuk2 for all u∈V, a.e. t∈(0, a) with α > 0) and gη∈ V0, we apply Theorem 3.4 in [42] and we conclude that there exists a unique solution cη∈W2(0, a) to Problem Pη 1. Next, we establish estimate (3.11). Taking cη(t) as a test function in Problem Pη 1, we get 1 2 d dtkcη(t)k2 H+kcη(t)k2 V+ (1 −R(η(t))) (γ0(cη(t)))2=η(t)γ0(cη(t)), for a.e. t∈(0, a). Taking into account the fact that 0 ≤1−R(η(t)) for t∈(0, a) and using the following version of the Cauchy inequality with ε r s ≤1 2ε2r2+ε2 2s2,∀r, s ∈R, ε > 0,(3.12) 106 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation we have 1 2 d dtkcη(t)k2 H+kcη(t)k2 V≤η(t)γ0(cη(t)) ≤1 2ε2(η(t))2+ε2 2(γ0(cη(t)))2, for a.e. t∈(0, a) with an arbitrary and positive ε. Using the trace inequality with constant Ctr (see (1.9) in Chapter 1) and integrating from 0 to t, we obtain kcη(t)k2 H+2 C2 tr Zt 0 (γ0(cη(s)))2ds ≤1 ε2Zt 0 (η(s))2ds +ε2Zt 0 (γ0(cη(s)))2ds +kc0k2 H, for all t∈[0, a]. Choosing ε= 1/Ctr, we have 1 C2 tr Zt 0 (γ0(cη(s)))2ds ≤C2 tr Zt 0 (η(s))2ds +kc0k2 Hfor all t∈[0, a]. Since η∈ C([0, a]), it follows that kγ0(cη)k2 L2(0,t)≤C4 tr tkηk2 C([0,t]) +C2 trkc0k2 H,for all t∈[0, a], hence applying the property (1.31) with ℘= 1/2, we conclude that (3.11) holds.  Intermediate ordinary differential equation. Let η∈ C([0, a]) and let cη∈ W2(0, a) be the unique solution to Problem Pη 1corresponding to η. Consider the following problem. Problem Pη 2. Given Γ0∈R, find Γη∈H1(0, a) such that dΓη dt (t) = γ0(cη(t))(1 −R(η(t))) −Γη(t),for a.e. t∈(0, a), Γη(0) = Γ0. For Problem Pη 2we have the following existence and uniqueness result. Lemma 3.2 Problem Pη 2has a unique solution Γη∈H1(0, a)given by Γη(t) = Γ0e−t+e−tZt 0 γ0(cη(s))(1 −R(η(s))) esds, for all t∈[0, a].(3.13) 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 107 Moreover, we have kΓηkC([0,a]) ≤ |Γ0|+√akγ0(cη)kL2(0,a).(3.14) Proof. Let us define the function F(t, r) = γ0(cη(t))(1 −R(η(t))) −r, for a.e. t∈(0, a),∀r∈R. It is clear that F(t, ·) is Lipschitz continuous for a.e. t∈(0, a) and F(·, r)∈L2(0, a) for all r∈R. The existence and uniqueness of the solution to Problem Pη 2follows from the classical theorem of Cauchy–Lipschitz which can be found in Theorem 2.1 of [43]. A short computation entails the formula (3.13). From (3.13), using the estimate 1−R(η(t)) ≤1 for t∈[0, a] and applying the H¨older inequality, we obtain |Γη(t)|≤|Γ0|+√tkγ0(cη)kL2(0,t),for all t∈[0, a]. Since Γη∈ C([0, a]), the last estimate implies (3.14).  Next, we define the operator Λ1:C([0, a]) →W2(0, a) which to η∈ C([0, a]) assigns the unique solution cη∈W2(0, a) to Problem Pη 1. Moreover, we introduce the operator Λ2:C([0, a]) ×W2(0, a)→H1(0, a) which to η∈ C([0, a]) and cη∈W2(0, a) assigns the unique solution Γη∈H1(0, a) to Problem Pη 2. From Lemmata 3.1 and 3.2, it follows that the operators Λ1and Λ2are well defined. Now we define the operator Λ: C([0, a]) →H1(0, a)⊂ C([0, a]) by Λ(η) = Λ2(η, Λ1(η)) for η∈ C([0, a]).(3.15) We turn to investigate the properties of Λ. Given r > 0, we introduce the following notation Ba(r) = {u∈ C([0, a]); kukC([0,a]) ≤r} and two constants 108 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation T∗=1 2C2 tr ,(3.16) r∗= 2|Γ0|+√2kc0kH,(3.17) where Ctr is the trace constant. In the next step, we show that there exists a ball in C([0, T∗]) which is invariant under the operator Λ. Lemma 3.3 For the operator Λdefined by (3.15), we have Λ(BT∗(r∗)) ⊂BT∗(r∗). Proof. It is enough to prove that if η∈BT∗(r∗), then Λ(η)∈BT∗(r∗). Let η∈BT∗(r∗), i.e. η∈ C([0, T∗]) and kηkC([0,T∗]) ≤r∗. From the estimates (3.11), (3.14) with a=T∗ and taking into account (3.16) and (3.17), we obtain kΓηkC([0,T∗]) ≤ |Γ0|+√T∗kγ0(cη)kL2(0,T ∗) ≤ |Γ0|+√T∗(C2 tr√T∗kηkC([0,T ∗]) +Ctrkc0kH) ≤ |Γ0|+C2 trT∗(2|Γ0|+√2kc0kH) + Ctr√T∗kc0kH=r∗. This means that kΓηkC([0,T∗]) ≤r∗, i.e. Λ(η)∈BT∗(r∗), which proves the lemma.  The next two lemmata show that the operator Λ is compact. Lemma 3.4 The operator Λ: BT∗(r∗)→BT∗(r∗)is Lipschitz continuous with respect to the C([0, T∗])- topology. Proof. We prove that Λ is Lipschitz continuous from BT∗(r∗) into itself. Indeed, for η1,η2∈BT∗(r∗), let cη1,cη2∈W2(0, T∗) be the unique solutions to Problem Pη 1 corresponding to η1and η2, respectively. Moreover, let Γη1, Γη2∈H1(0, T∗) be the unique solutions to Problem Pη 2related to η1,cη1and η2,cη2, respectively. Thus we have 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 109 h∂cη1 ∂t (t), viV0×V+ ((cη1(t), v)) + γ0(cη1(t))(1 −R(η1(t)))γ0(v) = η1(t)γ0(v), h∂cη2 ∂t (t), viV0×V+ ((cη2(t), v)) + γ0(cη2(t))(1 −R(η2(t)))γ0(v) = η2(t)γ0(v), for a.e. t∈(0, T∗),∀v∈V, cη1(0) = cη2(0) = c0. Subtracting the above equations and taking v=cη1(t)−cη2(t)∈V, a.e. t∈(0, T∗) as a test function, we get 1 2 d dtkcη1(t)−cη2(t)k2 H+kcη1(t)−cη2(t)k2 V +γ0(cη1(t))1−R(η1(t))−γ0(cη2(t))1−R(η2(t))γ0(cη1(t)−cη2(t)) = (η1(t)−η2(t)) γ0(cη1(t)−cη2(t)),for a.e. t∈(0, T∗). Hence, adding and subtracting the term γ0(cη1(t))(1 −R(η2(t)))γ0(cη1(t)−cη2(t)) in the left hand-side of the previous expression and taking into account the linearity of the trace operator, it follows that 1 2 d dtkcη1(t)−cη2(t)k2 H+kcη1(t)−cη2(t)k2 V+ (1 −R(η2(t))) (γ0(cη1(t)−cη2(t)))2 = (η1(t)−η2(t)) γ0(cη1(t)−cη2(t)) −γ0(cη1(t))R(η2(t)) −R(η1(t))γ0(cη1(t)−cη2(t)),for a.e. t∈(0, T∗). Using the inequality 1 −R(η2(t)) ≥0 and integrating the last equation from 0 to t, we have 1 2kcη1(t)−cη2(t)k2 H+Zt 0kcη1(s)−cη2(s)k2 Vds ≤Zt 0|η1(s)−η2(s)||γ0(cη1(s)−cη2(s))|ds +Zt 0|γ0(cη1(s))||R(η2(s)) −R(η1(s))||γ0(cη1(s)−cη2(s))|ds, 110 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation for all t∈[0, T∗]. Next, taking into account the trace inequality and the fact that the truncation operator Ris 1-Lipschitz continuous, from the Cauchy inequality with positive ε, see expression (3.12), it follows 1 C2 tr Zt 0|γ0(cη1(s)−cη2(s))|2ds ≤1 2ε2Zt 0|η1(s)−η2(s)|2ds +ε2Zt 0|γ0(cη1(s)−cη2(s))|2ds +1 2ε2Zt 0|γ0(cη1(s))|2|η1(s)−η2(s)|2ds, for all t∈[0, T∗]. Choosing ε2=1 2C2 tr , we obtain 1 2C2 tr kγ0(cη1−cη2)k2 L2(0,t)≤C2 tr tkη1−η2k2 C([0,t]) +C2 tr kγ0(cη1)k2 L2(0,t)kη1−η2k2 C([0,t]), for all t∈[0, T∗]. Hence kγ0(cη1−cη2)k2 L2(0,T∗)≤2C4 tr kη1−η2k2 C([0,T∗])T∗+kγ0(cη1)k2 L2(0,T∗).(3.18) Now, subtracting the two equations obtained from (3.13) for η=η1,cη=cη1and for η=η2,cη=cη2, respectively, we find Γη1(t)−Γη2(t) = e−tZt 0γ0(cη1(s))1−R(η1(s))−γ0(cη2(s))1−R(η2(s))esds, for all t∈[0, T∗]. Adding and subtracting the term γ0(cη1(t))(1 −R(η2(t))) under the integral in the last equation, we get |Γη1(t)−Γη2(t)| ≤ Zt 0|γ0(cη1(s))||R(η2(s)) −R(η1(s))|ds +Zt 0|γ0(cη1(s)−cη2(s))||1−R(η2(s))|ds, for all t∈[0, T∗]. Exploiting the fact that 0 ≤R(η(·)) ≤1 and, again, the property that Ris 1-Lipschitz continuous, by the H¨older inequality we have 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 111 |Γη1(t)−Γη2(t)| ≤ √tkγ0(cη1−cη2)kL2(0,t)+√tkγ0(cη1)kL2(0,t)kη1−η2kC([0,t]), for all t∈[0, T∗]. Consequently kΓη1−Γη2kC([0,T∗]) ≤√T∗kγ0(cη1−cη2)kL2(0,T ∗) +√T∗kγ0(cη1)kL2(0,T∗)kη1−η2kC([0,T ∗]).(3.19) Using (3.18) in (3.19), we arrive at the following estimate kΓη1−Γη2kC([0,T∗]) ≤√2T∗C2 trkη1−η2kC([0,T ∗]) √T∗+kγ0(cη1)kL2(0,T∗) +√T∗kγ0(cη1)kL2(0,T∗)kη1−η2kC([0,T ∗]). Finally, using the estimate (3.11), we deduce that there exists a constant C > 0 depending on T∗,r∗and on the problem data such that kΓη1−Γη2kC([0,T∗]) ≤Ckη1−η2kC([0,T ∗]). Hence, it follows that the operator Λ is Lipschitz continuous from BT∗(r∗) into itself.  Lemma 3.5 The operator Λ: BT∗(r∗)→BT∗(r∗)is compact with respect to the C([0, T∗])- topology. Proof. From Lemma 3.4 we know that Λ is continuous. In order to prove the compactness of the operator Λ, let Bbe a bounded subset of BT∗(r∗). We show that Λ(B) is relatively compact in C([0, T∗]). Indeed, it is clear from Lemma 3.3 that all functions from Λ(B) are norm-bounded by r∗, thus Λ(B) is equibounded. We prove that the set Λ(B) is equicontinuous, that is, for every ε > 0 there exists ¯ δ > 0 such that for all t1, 118 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation Indeed, the fact that Γ(t)≥0for all t∈[0, T]follows from the proof of Lemma 3.6. Moreover, proceeding as in that proof, it is easy to show that Γ(t)≤(1 −σ)Γmfor all t∈[0, T]. In fact, we suppose that there exists t∗∈[0, T]such that Γ(t∗)>(1 −σ)Γm and we denote by t∗∗ = max{t∈[0, t∗]|Γ(t) = (1 −σ)Γm}. Then, Γ(t)>(1 −σ)Γm for all t∈(t∗∗, t∗). Then, using the ordinary differential equation of Problem PLH R, we obtain dΓ dt (t) = ka Lγ0(c(t))1−(1 −σ)Γm Γm−kd LΓ(t),for a.e. t∈(t∗∗, t∗). Integrating from t∗∗ to t, we have Γ(t) = Γ(t∗∗) + Zt t∗∗ σ ka Lγ0(c(s)) −kd LΓ(s)ds <(1 −σ)Γm+σ ka Lkγ0(c)kL∞(0,T)(t−t∗∗)−kd L(1 −σ)Γm(t−t∗∗) <(1 −σ)Γm,for all t∈[t∗∗, t∗], and then, we arrive to a contradiction. In the following two lemmata we show that if γ0(c)∈L∞(0, T), then the similar properties as in Lemma 3.6 hold also for the solutions to Problem PLH W. Lemma 3.7 If (c, Γ) ∈W+ 2(0, T)×H1(0, T)solves Problem PLH Wwith Γ0∈[0,Γm] and γ0(c)∈L∞(0, T), then Γ(t)∈[0,Γm]for all t∈[0, T]. Proof. First, we show that Γ(t)≤Γmfor all t∈[0, T]. Analogously to the proof of Lemma 3.6, we proceed by contradiction. Assume that Γ(t∗)>Γmfor some t∗∈[0, T] and define t∗∗ = max{t∈[0, t∗]|Γ(t)=Γm}. We have 0 ≤t∗∗ < t∗. Obviously, for t∈(t∗∗, t∗) it follows that Γ(t)>Γm. Moreover, for t∈(t∗∗, t∗), we have Γ(t) = Γ(t∗∗) + Zt t∗∗ ka L Γm γ0(c(s)) (Γm−Γ(s)) −kd LΓ(s)ds ≤Γm+ka L Γmkγ0(c)kL∞(0,T)Zt t∗∗ |Γm−Γ(s)|ds −kd LΓm(t−t∗∗). 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 119 If kγ0ckL∞(0,T)= 0, then Γ(t)≤Γmfor t∈(t∗∗, t∗). Otherwise, we choose ε= kd LΓ2 m ka Lkγ0(c)kL∞(0,T ). Since Γ ∈ C([0, T]), there exists ¯ δ > 0 such that if |t−t∗∗| ≤ ¯ δ, then |Γ(t∗∗)−Γ(t)| ≤ ε. For t∈(t∗∗,min{t∗, t∗∗ +¯ δ}), we have Γ(t)≤Γm+ka L Γmkγ0(c)kL∞(0,T)(t−t∗∗)kd LΓ2 m ka Lkγ0(c)kL∞(0,T)−kd LΓm(t−t∗∗) = Γm. So we obtain a contradiction. Next, we prove that Γ(t)≥0 for all t∈[0, T]. Suppose, by contradiction, that Γ(t∗)<0 for some t∗∈[0, T] and define t∗∗ = max{t∈[0, t∗]| Γ(t)=0}. We have 0 ≤t∗∗ < t∗and Γ(t)<0 for all t∈(t∗∗, t∗). Moreover, for all t∈[t∗∗, t∗] we have Γ(t) = Γ(t∗∗) + Zt t∗∗ ka L Γm γ0(c(s)) (Γm−Γ(s)) −kd LΓ(s)ds ≥0, which is a contradiction. The proof of the lemma is complete.  Lemma 3.8 Assume that σ= 0 in the definition (3.9) of the truncation operator and Γ0∈[0,Γm]. Let (c1,Γ1)∈W+ 2(0, T)×H1(0, T)with γ0(c1)∈L∞(0, T)and (c2,Γ2)∈W+ 2(0, T)×H1(0, T)with γ0(c2)∈L∞(0, T)be two solutions to Problem PLH W. Then c1(t) = c2(t)and Γ1(t) = Γ2(t)for a.e. t∈(0, T). Proof. It follows from Lemma 3.7 that Γ1(t)∈[0,Γm] and Γ2(t)∈[0,Γm] for all t∈[0, T]. Hence, (c1,Γ1) and (c2,Γ2) are two solutions to Problem PLH R. Now, the conclusion is a consequence of Theorem 3.1 and the proof of the lemma is complete.  3.1.4 Fully discrete approximation: a priori error estimate In this section, we consider a fully discrete approximation of Problem PLH Rwhich is done following the same two steps explained in Chapter 1. Moreover, we also use the same notation introduced there. Without loss of generality, we suppose in this section that D=kd L=ka L= Γm= 1. Using a hybrid combination of both backward and forward Euler schemes, we consider the following fully discrete approximations of Problem PLH R. 120 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation Problem Phk LH. Find chk ={chk n}N n=0 ⊂Vhand Γhk ={Γhk n}N n=0 ⊂Rsuch that chk 0=ch 0,Γhk 0= Γ0,(3.25) and, for n= 1, . . . , N and for all vh∈Vh, it holds (δchk n, vh)H+ ((chk n, vh)) + f(γ0(chk n), R(Γhk n−1)) γ0(vh)=Γhk n−1γ0(vh),(3.26) δΓhk n=f(γ0(chk n), R(Γhk n)) −Γhk n,(3.27) where ch 0∈Vhis an appropriate approximation of the initial condition c0. Under the assumptions of Theorem 3.1, using the Lax-Milgram lemma, we easily deduce the existence of a unique discrete solution to Problem Phk LH. In the sequel, we derive an error estimate for the differences cn−chk nand Γn−Γhk n. Under the following additional regularity of the solution to Problem PLH R c∈ C([0, T]; V)∩C1([0, T]; H) and Γ ∈ C1([0, T]),(3.28) we get the result presented below. We note that assuming regularity (3.28), we have γ0(c)∈ C([0, T]). Theorem 3.2 Assume the hypotheses of Theorem 3.1 and the regularity conditions (3.28) hold. Then there exists a positive constant C, independent of the discretization parameters hand k, such that the following error estimate is satisfied, for all {vh n}N n=1 ⊂ Vh, max 0≤n≤Nkcn−chk nk2 H+k N X n=1 kcn−chk nk2 V+σ(γ0(cn−chk n))2+ max 0≤n≤N|Γn−Γhk n|2 ≤Ckc0−ch 0k2 H+C k N X n=1 kcn−vh nk2 V+k2+  ∂c ∂t(tn)−δcn   2 H+ max 0≤n≤NI2 n +Cmax 0≤n≤Nkcn−vh nk2 H+C N−1 X n=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H,(3.29) 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 121 where Indenotes the integration error given by In=Ztn 0 f(γ0(c(s)), R(Γ(s))) −Γ(s)ds −k n X j=1 [f(γ0(cj), R(Γj)) −Γj].(3.30) Proof. Taking v=cn−vh∈Vin the parabolic equation in Problem PLH Rat time t=tn, we find that ∂c ∂t(tn), cn−vhH + ((cn, cn−vh)) + f(γ0(cn), R(Γn))γ0(cn−vh) = Γnγ0(cn−vh),(3.31) for n= 1, 2, . . . , N, and using equation (3.26) we have, for all vh∈Vh, (δchk n, cn−chk n)H+ ((chk n, cn−chk n)) + f(γ0(chk n), R(Γhk n−1)) γ0(cn−chk n) = (δchk n, cn−vh)H+ ((chk n, cn−vh)) + Γhk n−1γ0(cn−chk n) −Γhk n−1γ0(cn−vh) + f(γ0(chk n), R(Γhk n−1)) γ0(cn−vh).(3.32) From (3.31) and (3.32), we obtain, for all vh∈Vh, ∂c ∂t(tn)−δchk n, cn−chk nH +kcn−chk nk2 V +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−chk n) =∂c ∂t(tn)−δchk n, cn−vhH + ((cn−chk n, cn−vh)) +(Γn−Γhk n−1)γ0(cn−chk n)−(Γn−Γhk n−1)γ0(cn−vh) +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−vh), and therefore, for all vh∈Vh, (δcn−δchk n, cn−chk n)H+kcn−chk nk2 V +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−chk n) 122 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation = (δcn−δchk n, cn−vh)H+ ((cn−chk n, cn−vh)) +(Γn−Γhk n−1)γ0(cn−chk n)−(Γn−Γhk n−1)γ0(cn−vh) −∂c ∂t(tn)−δcn, cn−chk nH +∂c ∂t(tn)−δcn, cn−vhH +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−vh), where we recall that δcn= (cn−cn−1)/k. Moreover, reminding the following property of the divided differences (δan−δbn, an−bn)H≥1 2kkan−bnk2 H−1 2kkan−1−bn−1k2 H, the previous equation reads 1 2kkcn−chk nk2 H+kcn−chk nk2 V +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−chk n) ≤1 2kkcn−1−chk n−1k2 H+ (δcn−δchk n, cn−vh)H+ ((cn−chk n, cn−vh)) +(Γn−Γhk n−1)γ0(cn−chk n)−(Γn−Γhk n−1)γ0(cn−vh) −∂c ∂t(tn)−δcn, cn−chk nH +∂c ∂t(tn)−δcn, cn−vhH +f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−vh), for all vh∈Vh. Now, using the equality f(γ0(cn), R(Γn)) −f(γ0(chk n), R(Γhk n−1))γ0(v) =f(γ0(cn), R(Γn)) −f(γ0(cn), R(Γhk n−1))γ0(v) +f(γ0(cn), R(Γhk n−1)) −f(γ0(chk n), R(Γhk n−1))γ0(v),for all v∈V, and taking into account the following estimates 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 123 f(γ0(cn), R(Γhk n−1)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−chk n)≥σ(γ0(cn−chk n))2, f(γ0(cn), R(Γn)) −f(γ0(cn), R(Γhk n−1))γ0(v)≤ |γ0(cn)||Γn−Γhk n−1||γ0(v)|, f(γ0(cn), R(Γhk n−1)) −f(γ0(chk n), R(Γhk n−1))γ0(cn−vh) ≤ |γ0(cn−chk n)||γ0(cn−vh)|, for all v∈V, it follows that 1 2kkcn−chk nk2 H+kcn−chk nk2 V+σ(γ0(cn−chk n))2 ≤1 2kkcn−1−chk n−1k2 H+ (δcn−δchk n, cn−vh)H+kcn−chk nkVkcn−vhkV +|Γn−Γhk n−1||γ0(cn−chk n)|+|Γn−Γhk n−1||γ0(cn−vh)| +  ∂c ∂t(tn)−δcn  Hkcn−chk nkH+  ∂c ∂t(tn)−δcn  Hkcn−vhkH +|γ0(cn−chk n)||γ0(cn−vh)|+|γ0(cn)||Γn−Γhk n−1||γ0(cn−chk n)| +|γ0(cn)||Γn−Γhk n−1||γ0(cn−vh)|,for all vh∈Vh. Using now the Cauchy inequality with a small parameter and the property (1.31) with ℘= 2, considering the regularity condition (3.28) and the fact that k·kH1(0,l)and k·kV are equivalent norms in Vand keeping in mind that |Γn−Γhk n−1|≤|Γn−Γn−1|+|Γn−1−Γhk n−1| ≤ C k kΓkC1([0,T ]) +|Γn−1−Γhk n−1|,(3.33) we find 1 2kkcn−chk nk2 H+αkcn−chk nk2 V+α σ (γ0(cn−chk n))2 ≤1 2kkcn−1−chk n−1k2 H+C(δcn−δchk n, cn−vh)H+kcn−vhk2 V +|Γn−1−Γhk n−1|2+|γ0(cn−vh)|2+  ∂c ∂t(tn)−δcn   2 H+kcn−chk nk2 H+k2, 124 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation for all vh∈Vh, where the positive constants αand C, which are small and large enough, respectively, are independent of the discretization parameters hand k. Thus, from the previous estimate, by an induction argument with respect to n, we obtain kcn−chk nk2 H+α k n X j=1 kcj−chk jk2 V+α k σ n X j=1 (γ0(cj−chk j))2 ≤ kc0−ch 0k2 H+C k n X j=1 (δcj−δchk j, cj−vh j)H+kcj−vh jk2 V+|Γj−1−Γhk j−1|2 +|γ0(cj−vh j)|2+  ∂c ∂t(tj)−δcj   2 H+kcj−chk jk2 H+k2,(3.34) for all {vh j}N j=1 ⊂Vh. Now, we turn to obtain some estimates on the numerical errors for the surface concentration. We integrate the ordinary differential equation in Problem PLH Rand we get Γn= Γ0+Ztn 0f(γ0(c(s)), R(Γ(s))) −Γ(s)ds. From (3.27) and using (3.25), we have Γhk n= Γhk n−1+k f(γ0(chk n), R(Γhk n)) −kΓhk n= Γ0+k n X j=1 f(γ0(chk j), R(Γhk j)) −Γhk j. Using the last two equations, we find that |Γn−Γhk n| ≤ In+k n X j=1 h|f(γ0(cj), R(Γj)) −f(γ0(chk j), R(Γhk j))|+|Γj−Γhk j|i,(3.35) where we recall that the integration error Inis defined by expression (3.30). Considering the relation |f(γ0(cj), R(Γj)) −f(γ0(chk j), R(Γhk j))| ≤ |f(γ0(cj), R(Γj)) −f(γ0(cj), R(Γhk j))| +|f(γ0(cj), R(Γhk j)) −f(γ0(chk j), R(Γhk j))| ≤ |γ0(cj)||Γj−Γhk j|+|γ0(cj−chk j)|, 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 125 where the fact that the operator R is 1-Lipschitz continuous and the inequality 1 − R(·)≤1 have been used; and applying several times inequality (1.31) with ℘= 2, taking into account the regularity c∈ C([0, T]; V) and Cauchy-Schwarz inequality and recalling that k N =T, from expression (3.35), we find that |Γn−Γhk n|2≤C I2 n+C k n X j=1 |Γj−Γhk j|2+C k n X j=1 |γ0(cj−chk j)|2.(3.36) Moreover, keeping in mind that k n X j=1 (δcj−δchk j, cj−vh j)H= (cn−chk n, cn−vh n)+(ch 0−c0, c1−vh 1) + n−1 X j=1 (cj−chk j, cj−vh j−(cj+1 −vh j+1))H ≤1 2kc0−ch 0k2 H+1 2kc1−vh 1k2 H+1 2kcn−chk nk2 H+1 2kcn−vh nk2 H + n−1 X j=1 kkcj−chk jk2 H+1 4kkcj−vh j−(cj+1 −vh j+1)k2 H, combining the estimates (3.34) and (3.35), and using (3.36) and the initial condition Γ0= Γhk 0, we have kcn−chk nk2 H+α k n X j=1 kcj−chk jk2 V+α k σ n X j=1 (γ0(cj−chk j))2+|Γn−Γhk n|2 ≤Ckc0−ch 0k2 H+C k n X j=1 kcj−vh jk2 V+k2+|Γj−Γhk j|2 +|γ0(cj−vh j)|2+  ∂c ∂t(tj)−δcj   2 H+kcj−chk jk2 H+C I2 n+Ckc1−vh 1k2 H +Ckcn−vh nk2 H+C n−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H. Using now the notation 126 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation ean:= kcn−chk nk2 H+α k n X j=1 kcj−chk jk2 V+α k σ n X j=1 (γ0(cj−chk j))2+|Γn−Γhk n|2, where ea0=kc0−ch 0k2 Hand considering the trace inequality, see (1.9), we find that ean≤Cea0+C k n X j=1 eaj+C k n X j=1 kcj−vh jk2 V+k2+  ∂c ∂t(tj)−δcj   2 H+C I2 n +Ckc1−vh 1k2 H+Ckcn−vh nk2 H+C n−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H. Denoting by egn:= ea0+k n X j=1 kcj−vh jk2 V+k2+  ∂c ∂t(tj)−δcj  +I2 n +kc1−vh 1k2 H+kcn−vh nk2 H+ n−1 X j=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H, we can conclude that ean≤Cegn+C k n X j=1 eaj. Finally, applying a discrete version of the Gronwall inequality, like in Section 1.4 of Chapter 1 (see, for instance, [4]), the result is achieved.  Estimate (3.29) is the basis for the convergence analysis. As an example, we state the following corollary under the assumption that the finite element space Vhis given by (1.41) and under further regularity conditions on the solution to the continuous problem: c∈ C([0, T]; H2(0, l)) ∩H1(0, T;V)∩H2(0, T;H).(3.37) Corollary 3.1 Assume the hypotheses of Theorem 3.2 and the regularity condition (3.37) hold. Then, the convergence of the algorithm in Problem Phk LH is linear, i.e. 3.1. Mixed kinetic adsorption: Langmuir-Hinshelwood equation 127 there exists a constant β > 0, independent of hand k, such that max 0≤n≤Nkcn−chk nkH+ max 0≤n≤N|Γn−Γhk n| ≤ β(h+k). Proof. Let πh:C([0, l]) →Vhdenote the standard finite element projection operator, and let us take vh j=πhcjfor j= 1, . . . , N. Moreover, assume that the discrete initial condition is given by ch 0=πhc0. Since c∈ C([0, T]; H2(0, l)), we obtain (see [10]) max 0≤n≤Nkcn−πhcnkV≤β h kckC([0,T];H2(0,l)), max 0≤n≤Nkcn−πhcnkH≤β h2kckC([0,T];H2(0,l)). Keeping in mind (1.43), using H¨older inequality and from the regularity hypothesis c∈H2(0, T;H), we get k N X n=1   ∂c ∂t(tn)−δcn   2 H≤k N X n=1 1 k2Ztn tn−1Ztn t  ∂2c ∂t2(s)  Hds dt2 ≤1 k N X n=1 Ztn tn−1 k1/2  ∂2c ∂t2  L2(tn−1,tn;H)dt2 ≤ N X n=1 k2  ∂2c ∂t2   2 L2(tn−1,tn;H) =k2  ∂2c ∂t2   2 L2(0,T;H). From the definition of the integration error In, we obtain In≤ n X j=1 Ztj tj−1f(γ0(c(s)), R(Γ(s))) −f(γ0(cj), R(Γj))ds + n X j=1 Ztj tj−1 (Γj−Γ(s)) ds.(3.38) Now, using the regularity condition (3.28) and the mean value theorem, on one hand we have n X j=1 Ztj tj−1|Γj−Γ(s)|ds ≤ kΓkC1([0,T]) n X j=1 Ztj tj−1|tj−s|ds ≤T k kΓkC1([0,T ]).(3.39) 134 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation in this section can be proved following the same techniques we use in the previous section, and so we omit all the proofs. We have decided to perform a detailed mathematical analysis of the model concerning the Langmuir-Hinshelwood equation and not its modification because the former is actually the classical model that appears in several chemical literature and, on the contrary, the model presented in this section is only a modification of the classical one proposed by two authors. 3.2.1 Model setting and its weak formulation Now, we are interested in the problem consisting of the system of equations (3.1)-(3.5) and (3.43). In order to simplify the notation and taking into account the function f, defined in (3.7), expression (3.43) can be written as follows: dΓ dt (t) = f(c(t, 0),Γ(t))e−BΓ(t) Γm−kd LΓ(t)e−BΓ(t) Γm, t > 0.(3.44) Therefore, we are concerned in analyzing problem (3.1)-(3.5) coupled with (3.44). As we did in the analysis of problem (3.1)-(3.5) and (3.8), here, to simplify the calculations, we also assume that cbequals zero. Assume that cis a smooth function which solves the problem we are considering. Multiplying equation (3.1) by smooth function zdefined in [0, l] such that z(l) = 0, integrating in (0, l), using the integration by parts formula and equations (3.3) and (3.44), we get Zl 0 ∂c ∂t(t, x)z(x)dx +DZl 0 ∂c ∂x(t, x)∂z ∂x(x)dx +f(c(t, 0),Γ(t))e−BΓ(t) Γmz(0) =kd LΓ(t)e−BΓ(t) Γmz(0), for a.e. t∈(0, T). Then, using (3.4), (3.5) and (3.44), we get the following weak formulation of the problem. Problem PmLH W. Given c0∈Hand Γ0∈R, find c∈W2(0, T) and Γ ∈H1(0, T) such 3.2. Mixed kinetic adsorption: modified Langmuir-Hinshelwood equation 135 that h∂c ∂t(t), viV0×V+D((c(t), v)) + f(γ0(c(t)),Γ(t))e−BΓ(t) Γmγ0(v) = kd LΓ(t)e−BΓ(t) Γmγ0(v), for a.e. t∈(0, T),∀v∈V, dΓ dt (t) = f(γ0(c(t)),Γ(t))e−BΓ(t) Γm−kd LΓ(t)e−BΓ(t) Γm,for a.e. t∈(0, T), c(0) = c0,Γ(0) = Γ0. Analogously as in the previous section, we need the truncation operator given in (3.9), and we work with the following truncated version of the modified LangmuirHinshelwood equation: dΓ dt (t) = f(c(t, 0), R(Γ(t)))e−BR(Γ(t)) Γm−kd LR(Γ(t))e−BR(Γ(t)) Γm, t > 0.(3.45) Proceeding as before, we define the following weak formulation of the truncated problem, associated to Problem PmLH W. Problem PmLH R. For given c0∈Hand Γ0∈R, find c∈W2(0, T) and Γ ∈H1(0, T) such that h∂c ∂t(t), viV0×V+D((c(t), v)) + f(γ0(c(t)), R(Γ(t)))e−BR(Γ(t)) Γmγ0(v) =kd LΓ(t)e−BR(Γ(t)) Γmγ0(v),for a.e. t∈(0, T),∀v∈V, dΓ dt (t) = f(γ0(c(t)), R(Γ(t)))e−BR(Γ(t)) Γm−kd LR(Γ(t))e−BR(Γ(t)) Γm,for a.e. t∈(0, T), c(0) = c0,Γ(0) = Γ0. 3.2.2 An existence and uniqueness result for Problem PmLH R In this section, we introduce an existence and uniqueness result for Problem PmLH R. Its proof is obtained proceeding as in the proof of Theorem 3.1 and so, here we only give a brief scheme of this proof indicating the main steps. 136 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation Theorem 3.3 Let D, ka L, kd Land Γmbe positive constants, B, Γ0∈Rand c0∈H. Then Problem PmLH Rhas a unique solution (c, Γ) ∈W2(0, T)×H1(0, T). The proof is done in two steps: we first split the truncated problem into two intermediate problems for which the existence and uniqueness of solution is obtained. Then, the application of Schauder fixed-point theorem leads to the desired result. In order to simplify the notation, in this section we assume that D=ka L=kd L= Γm= 1. Now, let abe a given positive constant representing an arbitrary time. Intermediate parabolic problem. Let η∈ C([0, a]) and consider the following parabolic problem: Problem e Pη 1. Given c0∈H, find cη∈W2(0, a) h∂cη ∂t (t), viV0×V+ ((cη(t), v)) + γ0(cη(t))(1 −R(η(t))) e−BR(η(t))γ0(v) =η(t)e−BR(η(t)) γ0(v),for a.e. t∈(0, T),∀v∈V, cη(0) = c0. Lemma 3.9 There exists a unique solution cη∈W2(0, a)to Problem e Pη 1. Moreover kγ0(cη)kL2(0,a)≤¯ C C2 tr√akηkC([0,a]) +Ctrkc0kH,(3.46) where Ctr is the trace constant given by (1.9) and ¯ C:= max{1, e−B}. Intermediate ordinary differential equation. Let η∈ C([0, a]) and cη∈W2(0, a) be the unique solution to Problem e Pη 1corresponding to η. We formulate the following: Problem e Pη 2.Given Γ0∈R, find Γη∈H1(0, a) such that dΓη dt (t) = γ0(cη(t))(1 −R(η(t)))e−BR(η(t)) −R(η(t))e−BR(η(t)),for a.e. t∈(0, a), Γη(0) = Γ0. The following lemma establishes the existence and uniqueness of solution to Problem e Pη 2. 3.2. Mixed kinetic adsorption: modified Langmuir-Hinshelwood equation 137 Lemma 3.10 There exists a unique solution Γη∈H1(0, a)to Problem e Pη 2given by Γη(t)=Γ0+Zt 0γ0(cη(s))(1 −R(η(s)))e−BR(η(s)) −R(η(s))e−BR(η(s))ds, (3.47) for all t∈[0, a]. Moreover, the following estimate holds kΓηkC([0,a]) ≤ |Γ0|+¯ C√akγ0(cη)kL2(0,a)+¯ C a. (3.48) We define the operator e Λ1:C([0, a]) →W2(0, a) as follows, η→e Λ1(η) = cη, where cηis the unique solution to Problem e Pη 1. Moreover, we consider the operator e Λ2:C([0, a]) ×W2(0, a)→H1(0, a) given by (η, cη)→e Λ2(η, cη) = Γη, being Γηthe unique solution to Problem e Pη 2. We remark that Lemmata 3.9 and 3.10 guarantee both operator e Λ1and operator e Λ2are well defined. Furthermore, we introduce the operator e Λ: C([0, a]) →H1(0, a)⊂ C([0, a]) by e Λ(η) = e Λ2(η, e Λ1(η)) for η∈ C([0, a]).(3.49) For e T=1 2¯ C2C2 tr and er= 2|Γ0|+√2kc0kH+1 ¯ C C2 tr we obtain the following properties of the operator e Λ. Lemma 3.11 The operator e Λmaps the ball Be T(er)into itself. Moreover, it is compact with respect to the C([0,e T])-topology. We remark that the existence of solution to Problem PmLH Rfollows from Lemma 3.11 and the Schauder fixed-point theorem. Moreover, the uniqueness of solution can be demonstrated arguing as in the proof of Theorem 3.1. 138 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation 3.2.3 An existence and uniqueness result for Problem PmLH W In this section, we present three results that establish relations between the solutions of Problems PmLH Rand PmLH W. The following two lemmata study the existence of solution to Problem PmLH Wunder the assumption that c∈W+ 2(0, T). The third one deals with the uniqueness of solution to Problem PmLH W. Some proofs are omitted because they follow proceeding as in the proofs of the corresponding lemmata in Section 3.1.3. Lemma 3.12 Assume that σ= 0 in the definition of the truncation operator (3.9) and that the initial condition Γ0∈[0,Γm]. If (c, Γ) ∈W+ 2(0, T)×H1(0, T)is a solution to Problem PmLH R, then Γ(t)∈[0,Γm]for all t∈[0, T], and, consequently, (c, Γ) is a solution to Problem PmLH W. The following lemma states that a solution to Problem PmLH Wverifying that γ0(c)∈ L∞(0, T) is also a solution to Problem PmLH R. Lemma 3.13 If (c, Γ) ∈W+ 2(0, T)×H1(0, T)solves Problem PmLH Wwith Γ0∈[0,Γm] and γ0(c)∈L∞(0, T), then Γ(t)∈[0,Γm]for all t∈[0, T]. Proof. The fact that Γ(t)≥0 directly follows from the proof of Lemma 3.7 in the previous section. So, here we only show that Γ(t)≤Γmfor all t∈[0, T]. Assume now that there exists t∗∈[0, T] such that Γ(t∗)>Γm. We define t∗∗ = max{t∈[0, t∗]| Γ(t) = Γm}and we note that this definition makes sense since Γ is a continuous function from [0, T] to R. Obviously, for all t∈(t∗∗, t∗), we have Γ(t)>Γm. Integrating the ordinary differential equation in Problem PmLH Wfrom t∗∗ to tand taking into account that γ0(c)∈L∞(0, T), we obtain Γ(t)≤Γm+ka L Γmkγ0(c)kL∞(0,T)Zt t∗∗ (Γm−Γ(s))e−BΓ(s) Γmds −kd LZt t∗∗ Γ(s)e−BΓ(s) Γmds, for all t∈(t∗∗, t∗).(3.50) Now, two different cases are considered depending on the sign of B(the case B= 0 is done in the previous section). If B > 0, then the following estimates hold for all 3.2. Mixed kinetic adsorption: modified Langmuir-Hinshelwood equation 139 t∈(t∗∗, t∗), e−BΓ(t) Γm≤e−B,−e−BΓ(t) Γm≤ −e−BkΓkC([0,T ]) Γm. Therefore, from (3.50), we have for all t∈(t∗∗, t∗) Γ(t)≤Γm+ka L Γmkγ0(c)kL∞(0,T)e−BZt t∗∗ (Γm−Γ(s))ds −kd LΓme−BkΓkC([0,T ]) Γm(t−t∗∗). If kγ0(c)kL∞(0,T)= 0, then Γ(t)≤Γmfor t∈(t∗∗, t∗). Otherwise, choosing ε=kd LΓ2 me−BkΓkC([0,T ]) Γm ka Lkγ0(c)kL∞(0,T)e−B, there exists ¯ δ > 0 such that if |t−t∗∗| ≤ ¯ δ, then |Γ(t∗∗)−Γ(t)| ≤ ε. For t∈ (t∗∗,min{t∗, t∗∗ +¯ δ}), we get Γ(t)≤Γm+ka L Γmkγ0(c)kL∞(0,T)e−B(t−t∗∗)kd LΓ2 me−BkΓkC([0,T ]) Γm ka Lkγ0(c)kL∞(0,T)e−B −kd LΓme−BkΓkC([0,T ]) Γm(t−t∗∗)=Γm. So we have a contradiction. The case B < 0 is analogous.  Lemma 3.14 Assume that σ= 0 in the definition (3.9) of the truncation operator and Γ0∈[0,Γm]. Let (c1,Γ1)∈W+ 2(0, T)×H1(0, T)with γ0(c1)∈L∞(0, T)and (c2,Γ2)∈ W+(0, T)×H1(0, T)with γ0(c2)∈L∞(0, T)be two solutions to Problem PmLH W. Then c1(t) = c2(t)and Γ1(t)=Γ2(t)for a.e. t∈(0, T). 3.2.4 Fully discrete approximation: a priori error estimate In this section, we introduce a fully discrete approximation to Problem PmLH Rand, in doing so, the same notations as the ones taken in Section 1.4 are used here. We also assume, without loss of generality, that D=kd L=ka L= Γm= 1. Moreover, by using the finite element method to obtain the spatial discretization and a hybrid combination of both backward and forward Euler schemes to discretize the time derivatives, we deal 140 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation with the following fully discrete approximation of Problem PmLH R: Problem Phk mLH. Find chk ={chk n}N n=0 ⊂Vhand Γhk ={Γhk n}N n=0 ⊂Rsuch that chk 0=ch 0,Γhk 0= Γ0,(3.51) and, for n= 1, . . . , N, and all vh∈Vh, it holds (δchk n, vh)H+ ((chk n, vh)) + f(γ0(chk n), R(Γhk n−1))e−BR(Γhk n−1)γ0(vh) = Γhk n−1e−BR(Γhk n−1)γ0(vh),(3.52) δΓhk n=f(γ0(chk n), R(Γhk n−1))e−BR(Γhk n−1)−R(Γhk n−1)e−BR(Γhk n−1),(3.53) where ch 0∈Vhis an appropriate approximation of the initial condition c0. Under the assumptions of Theorem 3.3, Lax-Milgram lemma guarantees the existence and uniqueness of solution to Problem Phk mLH. The results presented in what follows focus on deriving an error estimate for the differences cn−chk nand Γn−Γhk n. Their proofs are omitted since they follow the same ideas and techniques proposed in Section 3.1.5. Assuming the additional regularity conditions of the solution to Problem PmLH Rgiven in (3.28), we get the following result. Theorem 3.4 Assume the hypotheses of Theorem 3.3 and the regularity condition (3.28) hold. Then, there exists a constant β > 0, independent of the discretization parameters hand k, such that the following error estimate is satisfied, for all {vh n}N n=1 ⊂ Vh, max 0≤n≤Nkcn−chk nk2 H+k N X n=1 kcn−chk nk2 V+σ(γ0(cn−chk n))2+ max 0≤n≤N|Γn−Γhk n|2 ≤βkc0−ch 0k2 H+β k N X n=1 kcn−vh nk2 V+k2+  ∂c ∂t(tn)−δcn   2 H+ max 0≤n≤NI2 n +βmax 0≤n≤Nkcn−vh nk2 H+β N−1 X n=1 1 kkcj−vh j−(cj+1 −vh j+1)k2 H,(3.54) 3.2. Mixed kinetic adsorption: modified Langmuir-Hinshelwood equation 141 where Inis the integration error defined now as follows, In=Ztn 0f(γ0(c(s)), R(Γ(s)))e−B(RΓ(s)) −R(Γ(s))e−B(RΓ(s))ds −k n X j=1 f(γ0(cj), R(Γj))e−BR(Γj)−R(Γj)e−BR(Γj). As an example of the convergence given by estimate (3.54), we state the following corollary under the assumption that the finite element space Vhis given by (1.41) and under further regularity condition on the solution to the continuous problem given by (3.37). Corollary 3.2 Assume the hypotheses of Theorem 3.4 and the regularity condition (3.37) hold. Then, the convergence of the algorithm in Problem Phk mLH is linear, i.e. there exists a constant β > 0, independent of hand k, such that max 0≤n≤Nkcn−chk nkH+ max 0≤n≤N|Γn−Γhk n| ≤ β(h+k). 3.2.5 Numerical results The numerical scheme, implemented in MATLAB, for approximating Problem Phk mLH is presented in this section. Moreover, some numerical simulations are introduced in order to show the behavior of this model. Given the finite element space defined in (1.41) and given chk n−1∈Vhand Γhk n−1∈R for n= 1,2, . . . , N, we calculate the discrete concentration at time t=tnof surfactant, denoted by chk n, by using equation (3.52), but here with generic constants; that is to say, we find the solution of the following linear problem (chk n, vh)H+D k((chk n, vh)) + k f(γ0(chk n), R(Γhk n−1))e−BR(Γhk n−1) Γmγ0(vh) = (chk n−1, vh)H+k kd LΓhk n−1e−BR(Γhk n−1) Γmγ0(vh),∀vh∈Vh. Once chk nis known, the discrete surface concentration Γhk nis calculated from equation 142 Mixed kinetic-diffusion model with the Langmuir-Hinshelwood equation (3.53) by the expression Γhk n= Γhk n−1+k f(γ0(chk n), R(Γhk n−1))e−BR(Γhk n−1)−k kd LR(Γhk n−1)e−BR(Γhk n−1). Besides, we describe below the implemented algorithm which solves this problem. 1. Initial time step. At the beginning both chk 0and Γ0are given. 2. (n)th time step. The bulk and surface concentrations at time tn−1,chk n−1and Γhk n−1, respectively, are known. Then, at time tn,chk nand Γhk nare obtained using the following algorithm: (a) We calculate chk nby solving the following linear problem: Zl 0 chk nvhdx +Dk Zl 0 ∂chk n ∂x ∂vh ∂x dx +k ka Lγ0(chk n)1−R(Γhk n−1) Γme−BR(Γhk n−1) Γmγ0(vh) =Zl 0 chk n−1vhdx +k kd LΓhk n−1e−BR(Γhk n−1) Γmγ0(vh),∀vh∈Vh. (b) Then, the value of Γhk nis determined by the formula: Γhk n= Γhk n−1+k f(γ0(chk n), R(Γhk n−1))e−BR(Γhk n−1) Γm−k kd LR(Γhk n−1)e−BR(Γhk n−1) Γm. This algorithm has been implemented on a 3.2 GHz PC using MATLAB, and a typical run (h=k= 0.01) takes about 0.047 seconds of CPU time. First example: numerical convergence We consider the following test problem: ∂c ∂t(t, x)−5∂2c ∂x2(t, x)=0, t ∈(0,0.1), x ∈(0,1), 5∂c ∂x(t, 0) = f(c(t, 0), R(Γ(t)))e−BRΓ(t)−kd LR(Γ(t))e−BR(Γ(t)), t ∈(0,0.1), c(t, 1) = 1, t ∈(0,0.1), c(0, x) = c0(x), x ∈(0,1), 3.2. Mixed kinetic adsorption: modified Langmuir-Hinshelwood equation 143 with the initial condition c0(x) = 1. This problem corresponds to Problem PmLH Rwith the following data: l= 1, T = 0.1, cb= 1, D = 5, ka L= 1, kd L= 0.25, B=−1,Γm= 1,Γ0= 0. Choosing the solution obtained with parameters h= 1/16384 and k= 10−6as the “exact solution”, c, the numerical errors given by max 1≤n≤Nkcn−chk nkH+|Γn−Γhk n| are presented in Table 3.2 for several values of the discretization parameters hand k. In Figure 3.4, the error with respect to the value of parameter h+kis plotted. It can be seen that the linear convergence is achieved as Corollary 3.2 states. h↓k→0.01 0.005 0.002 0.001 0.0005 1/8 20.251566 9.730384 3.196324 1.066322 0.438502 1/16 21.179591 10.672453 4.136028 1.917607 0.802959 1/32 21.41244 10.909478 4.375164 2.156568 1.039592 1/64 21.4707 10.968818 4.435143 2.216719 1.099735 1/128 21.48527 10.983659 4.450149 2.231776 1.114816 1/256 21.48891 10.987369 4.453902 2.235543 1.118588 1/512 21.489825 10.988297 4.454839 2.236485 1.119532 1/1024 21.490052 10.988529 4.455074 2.236721 1.119768 1/2048 21.490109 10.988587 4.455133 2.2367791 1.119826 1/4096 21.490123 10.988601 4.455148 2.2367939 1.119843 Table 3.2: Numerical errors (×104) for several time and spatial discretization parameters.