EINDHOVEN UNIVERSITY OF TECHNOLOGY

EINDHOVEN UNIVERSITY OF TECHNOLOGY Department of Mathematics and Computer Science

CASA-Report 09-25 August 2009 Analysis of a two-scale system for gas-liquid reactions with non-linear Henry-type transfer by A. Muntean, M. Neuss-Radu

Centre for Analysis, Scientific computing and Applications Department of Mathematics and Computer Science Eindhoven University of Technology P.O. Box 513 5600 MB Eindhoven, The Netherlands ISSN: 0926-4507

ANALYSIS OF A TWO-SCALE SYSTEM FOR GAS-LIQUID REACTIONS WITH NON-LINEAR HENRY-TYPE TRANSFER ADRIAN MUNTEAN∗ AND MARIA NEUSS-RADU† Abstract. In this paper, we consider a coupled two-scale nonlinear reaction-diffusion system modelling gas-liquid reactions. The novel feature of the model is the nonlinear transmission condition between the microscopic and macroscopic concentrations, given by a nonlinear Henry-type transfer function. The solution is approximated by using a Galerkin method adapted to the multiscale form of the system. This approach leads to existence and uniqueness of the solution, and can also be used for numerical computations for a larger class of nonlinear multiscale problems. Key words. non-linear reaction-diffusion systems, multiscale Galerkin approximation, structured porous media, gas-liquid reactions, Henry’s law AMS subject classifications. 35K57, 35B27, 65N30, 80A32, 76R50

1. Introduction. Gas-liquid reactions occur in a wealth of physicochemical processes in chemical engineering [2, 3] or geochemistry [15], e.g. A minimal reactiondiffusion scenario for such reactions is the following: A chemical species A1 penetrates an unsaturated porous material thorough the air-filled parts of the pores and dissolves in water along the interfaces between water and air. Once arrived in water, the species A1 transforms into A2 and diffuses then towards places occupied by another yet dissolved diffusing species A3 . As soon as A2 and A3 meet each other, they react to produce water and various products (typically salts). This reaction mechanism can be described as follows k

A1 A2 + A3 −→ H2 O + products.

(1.1)

For instance, the natural carbonation of stone follows the mechanism (1.1), where A1 := CO2 (g), A2 := CO2 (aq), and A3 := Ca(OH)2 (aq), while the product of reaction is in this case CaCO3 (aq). This reaction mechanism is also encountered in catalysis, see [5] and references therein, or in civil engineering, e.g. in chemical corrosion scenarios as reported in [10, 11, 12]. This reaction-diffusion mechanism can be translated into a system of partial differential equations on the microscopic (pore) scale, together with transmission conditions at the interface separating the gas region from the region occupied by the liquid. Typical laws describing the mass transfer at these air-water interfaces are Henry-type laws. When some information about the pore geometry is available, e.g. industrially produced porous media have periodic geometry, these microscopic models can be homogenized to obtain multiscale descriptions of the underlying processes. These consist of a macroscopic system formulated on the entire domain, coupled with a microscopic system, formulated on the standard cell associated to the microstructure. The coupling between these two systems is given on one hand by a sink/source term appearing in the macroscopic equation and involving an integral operator over the ∗ CASA,

Center for Analysis, Scientific computing and Applications, Department of Mathematics, Eindhoven Institute of Technology, PO Box 513, 5600 MB, Eindhoven, The Netherlands, ([email protected]). † IWR, Heidelberg Graduate School of Mathematical and Computational Methods for the Sciences, University of Heidelberg, Im Neuenheimer Feld 294, D-69120 Heidelberg, Germany, ([email protected]). 1

2

A. MUNTEAN AND M. NEUSS-RADU

microscopic solution. On the other hand side by the boundary conditions for the microscopic solution, which is a (linear) Henry-type law involving the macroscopic concentration. For the derivation of such two-scale models (also called double porosity models) via homogenization techniques see e.g. [1, 7]. In our paper, we start from such a two-scale description of gas-liquid reactions in which the coupling between the microscopic and macroscopic problem is given by a nonlinear Henry-type law. Our study is partly motivated by some remarks from [4] mentioning the occurrence of nonlinear mass-transfer effects at air-water interfaces, partly by the fact that there are still a number of incompletely understood fundamental issues concerning gas-liquid reactions, and hence, a greater flexibility in the choice of the micro-macro coupling my help to identify the precise mechanisms. We assume that the microstructure consists of solid matrix and pores partially filled with water and partially filled with dry air; for details see Section 2.1. We assume the microstructure to be constant and wet, i.e. we do not account for any variations of the microstructure’s boundaries. The wetness of the porous material is needed to host the chemical reaction (1.1). We assume the wet parts of the pore to be static so that they do not influence the microscopic transport. The local geometry of the porous media we are interested in is given by a standard cell. The aim of the paper is to show that the two-scale model is well-posed, and to give an approach which can be used for the numerical computation of the solution. To this end, we use a Galerkin method to approximate the solutions of the two-scale system. We remark that this approach is of general interest for handling multiscale problems numerically. To show the convergence of the finite-dimensional approximates, we prove estimates which guarantee the compactness of the approximate solutions. Here, the challenging feature is to control the dependence of the cell solutions on the macroscopic variable. This is difficult because, on one hand we need strong L2 -convergence of the solutions and their traces to pass to the limit in the nonlinear terms. On the other hand, the macroscopic variable enters the cell problem just as a parameter, and therefore cannot be handled by standard methods. We remark that the techniques employed here are very much inspired by the analysis performed in [13]. There a multiscale Galerkin approach was developed to investigate transport and nonlinear reaction in domains separated by membranes. Numerical tools for multiscale problems for linear diffusion processes were developed also in [8]. There, the aim was to construct a sparse finite element discretization for the high-dimensional multiscale problem. The paper is organized as follows. We start with the statement of the problem and the assumptions on the data. The positivity and boundedness of weak solutions are proved in Sections 3 and 4. Next we show in Section 5 that there exists at most one weak solution of the considered multiscale system. In Section 6, we define the Galerkin approximations and show existence and uniqueness of the resulting systems of ordinary differential equations. The most delicate part of the paper is showing estimates which guarantee the compactness of the approximate solutions. Such estimates are proven assuming stability of the projections on finite dimensional subspaces with respect to appropriate norms. Finally, the convergence of the Galerkin approximates to the unique solution of the multiscale problem is shown.

2. Setting of the problem.

3

Gas-liquid reactions with nonlinear transfer

2.1. Geometry of the domain. Let Ω and Y be connected domains in R3 with Lipschitz continuous boundaries. We denote by λk the k-dimensional Lebesgue measure (k ∈ {2, 3}), and assume that λ3 (Ω) 6= 0 and λ3 (Y) 6= 0. Here, Ω is the macroscopic domain, while Y denotes the standard pore associated with the microstructure within Ω. We have that Y = Y ∪ Y g, where Y and Y g represent the wet region and the air-filled part of the standard pore respectively. Let Y and Y g be connected. The boundary of Y is denoted by Γ, and consists of two parts Γ = ΓR ∪ ΓN , where ΓR ∩ ΓN = ∅, and λ2y (ΓR ) 6= 0. Note that ΓR is the gas/liquid interface along which the mass transfer occurs, and λ2y denotes the surface measure on ∂Y . Furthermore, we denote by θ := λ3 (Y g ) the porosity of the medium. A possible geometry for our standard pore is illustrated in Figure 2.1.

Fig. 2.1. Standard pore Y. Typical shapes for Y, Y g ⊂ Y. ΓR is the interface between the gaseous part of the pore Y g and the wet region Y along which the mass transfer occurs.

2.2. Setting of the equations. Let us denote by S the time interval S =]0, T [ for a given T > 0. Let U , u and v denote the mass concentrations of the species A1 , A2 , and A3 respectively, see (1.1). The mass-balance for the vector (U, u, v) is described by the following two-scale system: Z θ∂t U (t, x) − D∆U (t, x) = − b(U (t, x) − u(t, x, y))dλ2y in S × Ω, (2.1) ΓR

∂t u(t, x, y) − d1 ∆y u(t, x, y) = −kη(u(t, x, y), v(t, x, y)) ∂t v(t, x, y) − d2 ∆y v(t, x, y) = −αkη(u(t, x, y), v(t, x, y))

in S × Ω × Y, (2.2) in S × Ω × Y, (2.3)

with macroscopic non-homogeneous Dirichlet boundary condition U (t, x) = U ext (t, x)

on S × ∂Ω,

(2.4)

4


and microscopic homogeneous Neumann boundary conditions ∇y u(t, x, y) · ny = 0 ∇y v(t, x, y) · ny = 0

on S × Ω × ΓN , on S × Ω × Γ.

(2.5) (2.6)

The coupling between the micro- and the macro-scale is made by the following nonlinear Henry-type transmission condition on ΓR −∇y u(t, x, y) · ny = −b(U (t, x) − u(t, x, y))

on S × Ω × ΓR .

(2.7)

The initial conditions U (0, x) = UI (x) u(0, x, y) = uI (x, y) v(0, x, y) = vI (x, y)

in Ω, in Ω × Y, in Ω × Y,

(2.8) (2.9) (2.10)

close the system of mass-balance equations. R Note that the sink/source term − ΓR b(U − u)dλ2y models the contribution in the effective equation (2.1) coming from mass transfer between air and water regions at microscopic level. The parameter k is the reaction constant for the competitive reaction between the species A2 and A3 , while α is the ratio of the molecular weights of these two species. 2.3. Assumptions on data and parameters. For the transport coefficients, we assume that (A1) D > 0, d1 > 0, d2 > 0. Concerning the micro-macro transfer and the reaction terms, we suppose (A2) The sink/source term b : R → R+ is globally Lipschitz, and b(z) = 0 if z ≤ 0. This implies that it exists a constant cˆ > 0 such that b(z) ≤ cˆz if z > 0. (A3) η : R × R → R+ is defined by η(r, s) := R(r)Q(s), where R, Q are globally Lipschitz continuous, with Lipschitz constants cR and cQ respectively. Furthermore, we assume that R(r) > 0 if r > 0 and R(r) = 0 if r ≤ 0, and similarly, Q(s) > 0 if s > 0 and Q(s) = 0 if s ≤ 0. Finally, k, α ∈ R, k > 0, and α > 0. For the initial and boundary functions, we assume 2 ∞ (A4) U ext ∈ H 1 (S, H 2 (Ω)) ∩ H 2 (S, L2 (Ω)) ∩ L∞ + (S × Ω), UI ∈ H (Ω) ∩ L+ (Ω), ext 1 2 2 ∞ UI − U (0, ·) ∈ H0 (Ω), uI , vI ∈ L (Ω, H (Y )) ∩ L+ (Ω × Y ). The classical choice for b in the literature on gas-solid reactions, see e.g. [3], is the linear one given by b : R → R+ , b(z) = cˆz for z > 0 and b(z) := 0 for z ≤ 0. However, there are applications, see e.g. [4], where extended Henry’s Law models are required. Our assumptions on b include self-limiting reactions, like e.g. MichaelisMenten kinetics. Typical reaction rates satisfying (A3) are power law reaction rates, sometimes also refereed to as generalized mass action laws ; see e.g. [2]. These laws have the form R(r) := rp and Q(s) := sq , where the exponents p ≥ 1 and q ≥ 1 are called partial orders of reaction. However, to fulfill the Lipschitz condition, for large values of the arguments the power laws have to replaced by Lipschitz functions. 2.4. Weak formulation. Our concept of weak solution is given in the following. Definition 2.1. A triplet of functions (U, u, v) with (U − U ext ) ∈ L2 (S, H01 (Ω)), ∂t U ∈ L2 (S × Ω), (u, v) ∈ L2 (S, L2 (Ω, H 1 (Y )))2 , (∂t u, ∂t v) ∈ L2 (S × Ω × Y )2 , is

5


called a weak solution of (2.1)–(2.10) if for a.e. t ∈ S the following identities Z Z Z Z d θU ϕ + D∇U ∇ϕ + b(U − u)ϕdλ2y dx = 0 dt Ω Ω Ω ΓR Z Z Z Z d uφ + d1 ∇y u∇y φ − b(U − u)φdλ2y dx dt Ω×Y Ω×Y Ω ΓR Z +k η(u, v)φ = 0 Z Z Z Ω×Y d vψ + d2 ∇y v∇y ψ + αk η(u, v)ψ = 0, dt Ω×Y Ω×Y Ω×Y

hold (2.11)

(2.12) (2.13)

for all (ϕ, φ, ψ) ∈ H01 (Ω) × L2 (Ω; H 1 (Y ))2 , and U (0) = UI in Ω,

u(0) = uI , v(0) = vI in Ω × Y.

3. Non-negativity of weak solutions. Let us first prove that weak solutions are non-negative. Lemma 3.1. Assume that hypotheses (A1)–(A4) hold, and that (U, u, v) is a weak solution of problem (2.1)–(2.10). Then, for a.e. (x, y) ∈ Ω × Y and all t ∈ S, we have U (t, x) ≥ 0,

u(t, x, y) ≥ 0,

v(t, x, y) ≥ 0.

(3.1)

Proof. We use here the notation u+ := max{0, u} and u− := max{0, −u}. Testing in (2.11)–(2.13) with (ϕ, φ, ψ) := (−U − , −u− , −v − ), we obtain Z Z Z Z d d θ|U − |2 + |u− |2 + D |∇U − |2 + d1 |∇y u− |2 dt Ω dt Ω×Y Ω Ω×Y Z Z Z −k η(u, v)u− + b(U − u)(u− − U − )dλ2y dx = 0 (3.2) Ω×Y Ω ΓR Z Z Z d |v − |2 d2 |∇v − |2 − αk η(u, v)v − = 0. (3.3) dt Ω×Y Ω×Y Ω×Y Note that by (A3) the last but one term of the r.h.s. of (3.2) and the last term of (3.3) vanish. We denote by H(·) the Heaviside function and estimate the last term of (3.2) as follows: Z Z Z Z b(U − u)(U − − u− )dλ2y dx ≤ b(U − u)U − dλ2y dx Ω ΓR Ω ΓR Z Z ≤ cˆ H(U − u)(U − u)U − dλ2y dx Ω ΓR Z Z = cˆ H(U − u)[U U − − u+ U + + u− U − ]dλ2y dx Ω ΓR Z Z Z 2 − 2 ≤ cˆλy (ΓR ) |U | + cˆ u− U − dλ2y dx Ω Ω ΓR Z Z 1 cˆ 2 − 2 ≤ cˆλy (ΓR )(1 + )||U ||L2 (Ω) + ||u− ||2H 1 (Y ) . (3.4) 2 2 Ω Y Now, we choose = dcˆ1 and apply Gronwall’s inequality in(3.2) and (3.3) to conclude the proof of the lemma.

6


4. Upper bounds for weak solutions. Next, we show that weak solutions are bounded. Lemma 4.1. If the hypotheses (A1)–(A4) hold, and (U, u, v) is a weak solution of problem (2.1)–(2.10), then for a.e. (x, y) ∈ Ω × Y and all t ∈ S, we have U (t, x) ≤ M1 ,

u(t, x, y) ≤ M2 ,

v(t, x, y) ≤ M3 ,

(4.1)

where M1 := max{||U ext ||L∞ (S×Ω) , ||UI ||L∞ (Ω) }, M2 := max{||uI ||L∞ (Ω×Y ) , M1 }, M3 := ||vI ||L∞ (Ω×Y ) . Proof. Choosing in (2.11)–(2.13) the test functions (ϕ, φ, ψ) := ((U − M1 )+ , (u − M2 )+ , (v − M3 )+ ), yields Z Z d d θ|(U − M1 )+ |2 + |(u − M2 )+ |2 + D |∇(U − M1 )+ |2 dt Ω dt Ω Z Z + 2 + d1 |∇y (u − M2 ) | + k η(u, v)(u − M2 )+ Ω×Y Ω×Y Z Z Z Z + 2 + b(U − u)(U − M1 ) dλy dx = b(U − u)(u − M2 )+ dλ2y dx (4.2) Ω

ΓR

Ω

ΓR

and d dt

Z

+ 2

Z

|(v − M3 ) | + d2 Ω×Y

Z + αk

|∇y (v − M3 )+ |2

Ω×Y

η(u, v)(v − M3 )+ = 0.

(4.3)

Ω×Y

Since the last two terms from the l.h.s of (4.2) and the last one from the l.h.s. of (4.3) are positive, the only term, which still needs to be estimated, is the term on the r.h.s of (4.2). We proceed as follows: Z Z Z Z + 2 b(U − u)(u − M2 ) dλy dx ≤ cˆ H(U − u)(U − u)(u − M2 )+ dλ2y dx Ω ΓR Ω ΓR Z Z Z Z ≤ cˆ H(U − u)(U − M1 )(u − M2 )+ dλ2y dx − cˆ H(U − u)|(u − M2 )+ |2 dλ2y dx Ω Γ Ω ΓR Z Z R Z Z + + 2 ≤ cˆ H(U − u)(U − M1 ) (u − M2 ) dλy dx − cˆ H(U − u)|(u − M2 )+ |2 dλ2y dx Ω Γ Ω ΓR Z ZR Z Z cˆ cˆ + 2 2 H(U − u)|(U − M1 ) | dλy dx − H(U − u)|(u − M2 )+ |2 dλ2y dx ≤ 2 Ω ΓR 2 Ω ΓR cˆ ≤ λ2y (ΓR )||U − M1 )+ ||2L2 (Ω) . (4.4) 2 The desired estimates on the solution now follow from (4.2), (4.3), and (4.4) by Gronwall’s inequality. Remark that, if the solution (U, u, v) is sufficiently smooth, and (A1)–(A3) hold, then one can use the technique from Lemma 3.1 of [5] to prove that the system (2.1)–(2.10) satisfies the classical maximum principle.


7

5. Uniqueness of weak solutions. Before starting with the proof of the existence of weak solutions for the problem (2.1)–(2.10), let us show that if they exist then weak solutions are unique. Proposition 5.1. If (A1)–(A4) hold, then the weak solution to (2.1)–(2.10) is unique. Proof. Let (Ui , ui , vi ), i ∈ {1, 2}, be arbitrary weak solutions of the problem (2.1)– 2 (2.10), thus satisfying for all (ϕ, φ, ψ) ∈ H01 (Ω) × L2 (Ω; H 1 (Y )) the equalities Z Z Z Z d θUi ϕ + D∇Ui ∇ϕ + b(Ui − ui )ϕdλ2y dx = 0 dt Ω Ω Ω ΓR Z Z Z Z Z d ui φ + d1 ∇y u i ∇y φ − b(Ui − ui )φdλ2y dx + k η(ui , vi )φ = 0 dt Ω×Y Ω×Y Ω ΓR Ω×Y Z Z Z d η(ui , vi )ψ = 0. d2 ∇y vi ∇y φ + αk vi ψ + dt Ω×Y Ω×Y Ω×Y We subtract the weak formulation written for (U2 , u2 , v2 ) from that one written for (U1 , u1 , v1 ) and choose in the result as test function (ϕ, φ, ψ) := (U2 − U1 , u2 − u1 , v2 − 2 v1 ) ∈ H01 (Ω) × L2 (Ω; H 1 (Y )) . We obtain: d dt

Z

d θ|U2 − U1 | + dt Ω 2

Z Ω×Y

d |u2 − u1 | + dt 2

Z

|v2 − v1 |2 Z |∇y (u2 − u1 )|2 + d2 |∇y (v2 − v1 )|2 Ω×Y

Z Z +D |∇(U2 − U1 )|2 + d1 Ω×Y Ω×Y ZΩ +k (η(u2 , v2 ) − η(u1 , v1 )) [(u2 − u1 ) + α(v2 − v1 )] Ω×Y Z Z + [b(U2 − u2 ) − b(U1 − u1 )][(U2 − U1 ) − (u2 − u1 )]dλ2y dx = 0 Ω

(5.1)

ΓR

The last two terms in (5.1), say I1 and I2 , can be estimated as follows. Z I1 ≤ |k (η(u2 , v2 ) − η(u1 , v1 )) [(u2 − u1 ) + α(v2 − v1 )] | Z Ω×Y ≤k |Q(v2 )(R(u2 ) − R(u1 )) [(u2 − u1 ) + α(v2 − v1 )] | ZΩ×Y +k |R(u1 )(Q(v2 ) − Q(v1 )) [(u2 − u1 ) + α(v2 − v1 )] | Ω×Y

Using the boundedness of the weak solutions and the Lipschitz continuity of the functions R and Q, we obtain Z ≤ kQmax |R(u2 ) − R(u1 )||(u2 − u1 ) + α(v2 − v1 )| Z Ω×Y + kRmax |Q(v2 ) − Q(v1 )||(u2 − u1 ) + α(v2 − v1 )| Ω×Y Z 1 3 α2 ≤ k(Qmax cR + Rmax cQ ) max |u2 − u1 |2 + |v2 − v1 |2 , , +α+ 2 2 2 Ω×Y

8


where Rmax := maxr∈[0,M2 ] R(r), and Qmax := maxs∈[0,M3 ] Q(s). We estimate I2 as follows. Z Z I2 ≤ | [b(U2 − u2 ) − b(U1 − u1 )][(U2 − U1 ) − (u2 − u1 )]dλ2y dx| Ω ΓR Z Z Z Z 2 2 ≤ 2ˆ c |U2 − U1 | dλy dx + 2ˆ c |u2 − u1 |2 dλ2y dx. (5.2) Ω

ΓR

Ω

ΓR

To estimate the second term in (5.2) we use the interpolation-trace inequality (5.3), see Proposition 5.2 below, with θ = 12 . We obtain Z Z 2ˆ c |u2 − u1 |2 λ2y dx Ω ΓR Z ≤C ||∇y (u2 − u1 )||L2 (Y ) + ||u2 − u1 ||L2 (Y ) ||u2 − u1 ||L2 (Y ) , ZΩ Z 2 ≤ ||∇y (u2 − u1 )||L2 (Y ) + c ||u2 − u1 ||2L2 (Y ) Ω

Ω

Choosing now = d21 , and applying Gronwall’s inequality, we conclude the statement of the Proposition. Proposition 5.2 (Theorem 5.9 in [14]). Let Y be a bounded domain in Rn with piecewise smooth Robin boundary ΓR and let u ∈ W 1,p (Y ). The following inequality holds: θ (5.3) ||u||Lq (ΓR ) ≤ C ||∇y u||Lp (Y ) + ||u||Lγ (Y ) ||u||1−θ Lr (Y ) , where θ =

p qn−r(n−1) q p(n+r)−nr

∈]0, 1[, 1 ≤ γ < ∞,

1≤r
p n−p n−p

1 ≤ r, q < ∞ if p = n or if p > n.

6. Global existence of weak solutions. In this section, we prove existence of a global weak solution of problem (2.1)–(2.10) by using the Galerkin method. This method has to be adapted to the two-scale form of our system, and yields an ansatz for the numerical treatment of a more general class of multiscale problems. One important aspect is the choice of the bases which are used to define finite dimensional approximations of the solution. Here, the structure of the basis elements reflects the two-scale structure of the solution; the basis elements on the domain Ω × Y are chosen as tensor products of basis elements on the macroscopic domain Ω and on the standard cell Y . The most challenging aspect in the analysis of the discretized problems is to control the nonlinear sink/source integral term coupling the macroscopic and microscopic problems. More precisely, enough compactness for the finite-dimensional approximations, with respect to both, the microscopic and macroscopic variables, is needed in order to be able to pass to the limit in the discretized problems. This is not straightforward due to the fact that the macroscopic variable enters the cell problem just as a parameter.


9

In [13], a Galerkin approach for the approximation of a multiscale reactiondiffusion system with transmission conditions was developed. We will adapt the concepts used there to the structure of our actual problem. Numerical tools for elliptic linear multiscale problems were developed in [8]. There, the aim was to construct a sparse finite element discretization for the high-dimensional multiscale problem. 6.1. Galerkin approximation. Global existence for the discretized problem. We introduce the following Schauder bases: Let {ξi }i∈N be a basis of L2 (Ω), with ξj ∈ H01 (Ω), forming an orthonormal system (say o.n.s.) with respect to L2 (Ω)norm. Furthermore, let {ζjk }j,k∈N be a basis of L2 (Ω × Y ), with ζjk (x, y) = ξj (x)ηk (y), where {ηk }k∈N is a basis of L2 (Y ), with ηk ∈ H 1 (Y ), forming an o.n.s. with respect to L2 (Y )-norm. Let us also define the projection operators on finite dimensional subspaces PxN , PyN associated to the bases {ξj }j∈N , and {ηk , }k∈N respectively. For (ϕ, ψ) of the form ϕ(x) =

X

aj ξj (x),

j∈N

X

ψ(x, y) =

bjk ξj (x)ηk (y),

j,k∈N

we define (PxN ϕ)(x) =

N X

aj ξj (x),

(6.1)

j=1

(PxN ψ)(x, y) =

N X X

bjk σj (x)ηk (y)

(6.2)

bjk σj (x)ηk (y).

(6.3)

j=1 k∈N

(PyN ψ)(x, y) =

N XX j∈N k=1

The bases {σj }j∈N , and {ηk }k∈N are chosen such that the projection operators PxN , PyN are stable with respect to the L∞ -norm and H 2 -norm; i.e. for a given function the L∞ norm and H 2 -norm of the truncations by the projection operators can be estimated by the corresponding norms of the function. Now, we look for finite-dimensional approximations of order N ∈ N for the functions U0 := U − U ext , u, and v, of the following form U0N (t, x)

=

N X

αjN (t)ξj (x),

(6.4)

j=1

uN (t, x, y) =

N X

N βjk (t)ξj (x)ηk (y),

(6.5)

N γjk (t)ξj (x)ηk (y),

(6.6)

j,k=1

v N (t, x, y) =

N X j,k=1

10


N N where the coefficients αjN , βjk , γjk , j, k = 1, . . . , N are determined by the following relations:

Z

θ∂t U0N (t)ϕdx

Z

D∇U0N (t)∇ϕdx =

+

(6.7) − b (U0N + U ext − uN )(t) dλ2 (y) + θ∂t U ext (t) + D∆U ext (t) ϕdx Ω ΓR Z Z N ∂t u (t)φ dxdy + d1 ∇y uN (t)∇y φ dxdy = (6.8) Ω×Y Ω×Y Z Z Z b (U0N + U ext − uN )(t) φ dλ2 (y)dx − k η uN (t), v N (t) φ dydx Ω ΓR Ω×Y Z Z N N d2 ∇y v (t)∇y ψ dydx = (6.9) ∂t v (t)ψ dydx + Ω×Y Ω×Y Z − αk η uN (t), v N (t) ψ dydx Ω

Ω

Z Z

Ω×Y

for all ϕ ∈ span{ξj : j ∈ {1, . . . , N }}, and φ, ψ ∈ span{ζjk : j, k ∈ {1, . . . , N }}, and

αjN (0) :=

Z

(UI − U ext (0))ξj dx, ZΩ Z N βjk (0) := uI ζjk dxdy, ZΩ ZY N γjk (0) := vI ζjk dxdy. Ω

(6.10) (6.11) (6.12)

Y

Taking in (6.7)-(6.9) as test functions ϕ = ξj , φ = ξjk , and ψ = ξjk , for j, k = 1, . . . , N , we obtain the following system of ordinary differential equations for the coefficients N N αN = (αjN )j=1,...,N , β N = (βjk ) j=1,...,N , and γ N = (γjk ) j=1,...,N : k=1,...,N

∂t αN (t) +

N X

k=1,...,N

Ai αiN (t) = F (αN (t), β N (t)),

(6.13)

i=1

∂t β N (t) +

N X

Bil βilN (t) = F˜ (αN (t), β N (t)) + G(β N (t), γ N (t)),

(6.14)

Cil γilN (t) = αG(β N (t), γ N (t)),

(6.15)

i,l=1

∂t γ N (t) +

N X i,l=1

11


where for j, k, i, l = 1, . . . , N we have Z (Ai )j := D∇ξi (x)∇ξj (x)dx, Ω Z (Bil )jk := d1 ∇ζil (x, y)∇ζjk (x, y) dydx, Ω×Y Z (Cil )jk := d2 ∇ζil (x, y)∇ζjk (x, y) dydx, Ω×Y Z 1 ± ± ± ∇ξrm (x)∇ξj0 (x)dx, (Crm )j := D± 2 Ω± Z Z 1 Fj := − b (U0N + U ext − uN )(t) dλ2y + θ∂t U ext (t) + D∆U ext (t) ξj (x)dx, θ Ω ΓR Z Z F˜jk := b (U0N + U ext − uN )(t) ζjk (x, y) dλ2 (y)dx, Ω ΓR Z Gjk := k η uN (t), v N (t) ζjk dydx. Ω×y

Due to the assumptions (A2)–(A3) on b and η, the functions F, F˜ , and G are globally Lipschitz continuous, and the Cauchy problem (6.10)-(6.15) has a unique solution 2 2 (αN , β N , γ N ) in C 1 ([0, T ])N × C 1 ([0, T ])N × C 1 ([0, T ])N . We conclude this section by proving the global Lipschitz property of F˜ , the proof of the Lipschitz continuity of F and G is similar. Let (U0N , uN ) and (W0N , wN ) be of 2 the form (6.4), (6.5), with coefficients (α1N , β1N ), (α2N , β2N ) ∈ RN × RN . We have: F˜jk (αN , β N ) − F˜jk (α2N , β2N ) Z Z 1 1 = b(U0N + U ext − uN )(t) − b(W0N + U ext − wN )(t) ζjk dλ2y dx Ω ΓR Z Z ≤ cˆ |(U0N − W0N ) − (uN − wN )||ζjk |dλ2y dx Ω

ΓR

Z Z |

= cˆ Ω

≤ cˆ

N X

ΓR

N X

(α1N (t) − α2N (t))i ξi −

i=1

+ cˆ

Z Z Ω

i=1

Z Z Ω

i=1 `=1

≤ cˆ max{cijk }

|ξi (x)ζjk (x, y)|dλ2y dx

ΓR

|(β1N (t) − β2N (t))i` | N X

(β1N (t) − β2N (t))i` ξi` ||ζjk |dλ2y dx

i=1 `=1

|(α1N (t) − α2N (t))i |

N X N X

N N X X

|ζi` (x, y)ζjk (x, y)|dλ2y dx

ΓR

|(α1N (t) − α2N (t))i | + cˆ max{ci`jk }

i=1

N X N X i=1 `=1

where the coefficients cijk and ci`jk are given by Z Z cijk := |ξi (x)ζjk (x, y)|dλ2y dx Ω ΓR Z Z ci`jk := |ζi` (x, y)ζjk (x, y)|dλ2y dx Ω

ΓR

|(β1N (t) − β2N (t))i` |,

12


for i, l, j, k = 1, . . . , N . Thus, we obtain |F˜ (α1N , β1N ) − F˜ (α2N , β2N )| ≤ c(N ) |α1N − α2N | + |β1N − β2N | .

(6.16)

6.2. Uniform estimates for the discretized problems. In this section, we prove uniform estimates for the solutions to the finite-dimensional problems. Based on this estimates, in the next section, we are able to pass in (6.7)–(6.9) to the limit N → ∞. Theorem 6.1. Assume that the projection operators PxN , PyN , defined in (6.1)(6.3), are stable with respect to the L∞ -norm and H 2 -norm, and that (A1)–(A4) are satisfied. Then the following statements hold: (i) The finite-dimensional approximations U0N (t), uN (t), and v N (t) are positive and uniformly bounded. More precisely, we have for a.e. (x, y) ∈ Ω × Y , all t ∈ S, and all N ∈ N 0 ≤ U0N (t, x) ≤ m1 ,

0 ≤ uN (t, x, y) ≤ m2 ,

0 ≤ v N (t, x, y) ≤ m3 , (6.17)

where m1 := 2||U ext ||L∞ (S×Ω) + ||UI ||L∞ (Ω) , m2 := max{||uI ||L∞ (Ω×Y ) , m1 }, m3 := ||vI ||L∞ (Ω×Y ) . (ii) There exists a constant c > 0, independent of N , such that ||U0N ||L∞ (S,H 1 (Ω)) + ||∂t U0N ||L2 (S,L2 (Ω)) ≤ c,

(6.18)

||uN ||L∞ (S,L2 (Ω;H 1 (Y ))) + ||∂t uN ||L2 (S,L2 (Ω;L2 (Y ))) ≤ c,

(6.19)

N

N

||v ||L∞ (S,L2 (Ω;H 1 (Y ))) + ||∂t v ||L2 (S,L2 (Ω;L2 (Y ))) ≤ c,

(6.20) (6.21)

Proof. (i) We consider the function U N := U0N + U ext . Then (U N , uN , v N ) satisfies the equations Z Z θ∂t U N (t)ϕdx + D∇U N (t)∇ϕdx = (6.22) Ω Ω Z Z b (U N − uN )(t) ϕdλ2 (y)dx − Ω ΓR Z Z N ∂t u (t)φ dxdy + d1 ∇y uN (t)∇y φ dxdy = (6.23) Ω×Y ZΩ×Y Z Z b (U N − uN )(t) φ dλ2 (y)dx − k η uN (t), v N (t) φ dydx Ω ΓR Ω×Y Z Z N N ∂t v (t)ψ dydx + d2 ∇y v (t)∇y ψ dydx = (6.24) Ω×Y Ω×Y Z − αk η uN (t), v N (t) ψ dydx Ω×Y

for all ϕ ∈ span{ξj : j ∈ {1, . . . , N }}, and φ, ψ ∈ span{ζjk : j, k ∈ {1, . . . , N }}, and U N (0) = PxN (UI − U ext (0)) + U ext (0) in Ω, uN (0) = PxN (PyN uI ) in Ω × Y. v N (0) = PxN (PyN vI ) in Ω × Y.

13


Using the stability of the projection operators with respect to the L∞ -norm, the proof follows the lines of the proofs of Lemma 3.1 and Lemma 4.1. (ii) Let us first take (ϕ, φ, ψ) = (U0N , uN , v N ) as test function in (6.7)–(6.9). Using the positivity of the solution, we get 1 d N ||v (t)||2L2 (Ω×Y ) + d2 ||∇y v N (t)||2L2 (Ω×Y ) ≤ 0, 2 dt θ d 1 d N ||U0N (t)||2L2 (Ω) + ||u (t)||2L2 (Ω×Y ) 2 dt 2 dt + D||∇U0N ||2L2 (Ω) + d1 ||∇y uN ||2L2 (Ω×Y ) Z Z + b(U0N + U ext − uN )(U0N − uN )dλ2y dx Ω ΓR Z + (θ∂t U ext + ∆U ext )U0N ≤ 0.

(6.25) (6.26)

Ω

To estimate the following term, we the Lipschitz property of b and the interpolationtrace inequality (5.3) Z Z Ω

b(U0N + U ext − uN )(U0N − uN )dλ2y dx

(6.27)

ΓR

Z Z ≤ cˆ Ω

|U0N + U0ext − uN ||U0N − uN |dλ2y dx

Γ

Z ZR

|U0N |2 + |U ext |2 + |uN |2 dλ2y dx Ω ΓR Z ≤ C ||U0N ||2L2 (Ω) + ||U ext ||2L2 (Ω) + C ||uN ||H 1 (Y ) ||uN ||L2 (Y ) Ω ≤ C ||U0N ||2L2 (Ω) + ||U ext ||2L2 (Ω) + ||∇y uN ||2L2 (Ω×Y ) + C()||uN ||2L2 (Ω×Y ) . ≤C

Taking := d21 in (6.27), inserting (6.27) in (6.26), and using the regularity properties of U ext , we obtain d d d ||U N (t)||2L2 (Ω) + ||uN (t)||2L2 (Ω×Y ) + ||v N (t)||2L2 (Ω×Y ) dt 0 dt dt + ||∇U0N (t)||2L2 (Ω) + ||∇y uN ||2L2 (Ω×Y ) + ||∇y v N ||2L2 (Ω×Y )

(6.28)

≤ C + ||U0N ||2L2 (Ω) + ||uN ||2L2 (Ω×Y ) . Integrating with respect to time, and applying Gronwall’s inequality yields the estimates ||U0N ||L∞ (S,L2 (Ω)) + ||∇U0N ||L2 (S,L2 (Ω)) ≤ c, N

(6.30)

N

(6.31)

||u ||L∞ (S,L2 (Ω×Y ))) + ||∇y u ||L2 (S,L2 (Ω×Y )) ≤ c, N

(6.29)

N

||v ||L∞ (S,L2 (Ω×Y ))) + ||∇y v ||L2 (S,L2 (Ω×Y )) ≤ c.

To obtain L∞ -estimates with respect to time of the gradients, and the estimates for the time derivatives, we differentiate with respect to time the weak formulation

14


(6.7)–(6.9), and test with (ϕ, φ, ψ) = (∂t U0N , ∂t uN , ∂t v N ) and obtain: Z Z Z 1 d 1 d 1 d θ|∂t U0N (t)|2 + |∂t uN (t)|2 + |∂t v N (t)|2 (6.32) 2 dt Ω 2 dt Ω×Y 2 dt Ω×Y Z Z Z +D |∇∂t U0N |2 + d1 |∇y ∂t uN |2 + d2 |∇y ∂t v N |2 Ω Ω×Y Ω×Y Z Z 2 + b0 (U0N + U ext − uN ) ∂t U0N + ∂t U ext − ∂t uN ∂t U0N − ∂t uN dλy t t Ω ΓR Z Z =− θ∂tt U ext ∂t U0N − D∆∂t U ext ∂t U0N Ω Ω Z Z 0 N N N 2 N N N N 2 +k R (u )Q(v ) |ut | + αut vt + k R(u)Q0 (v) uN . t vt + α|vt | Ω×Y

Ω×Y

Integrating this expression with respect to time, using the Lipschitz properties of the nonlinear terms b and η, and the regularity properties of U ext , as well as the interpolation-trace inequality (5.3), we obtain for all t ∈ S Z

|∂t U0N (t)|2 +

Ω Z t

Z

|∂t uN (t)|2 +

Z

Ω×Y

Z

|∇∂t U0N |2 Z0 Ω×Y Z ≤ |∂t U0N (0)|2 +

0

Z tZ |∇y ∂t uN |2 + |∇y ∂t v N |2 Ω×Y 0 Ω×Y Z N 2 N |∂t u (0)| + |∂t v (0)|2

Ω×Y

Z tZ

+C 1+

0

(6.33)

Z tZ

+

Ω

|∂t v N (t)|2

Ω×Y

Ω×Y

|∂t U0N (t)|2

Ω

Z tZ

N 2

Z tZ

|∂t u | +

+ 0

Ω×Y

N 2

|∂t v | 0

.

Ω×Y

To proceed, we have to estimate the norm of (∂t U0N , ∂t uN , ∂t v N ) at t = 0. For this purpose, we evaluate the weak formulation (6.7)–(6.9) at t = 0, and test with (∂t U0N (0), ∂t uN (0), ∂t v N (0)). We obtain Z

θ|∂t U0N (0)|2

Ω

Z +

D∇U0N (0)∇∂t U0N (0)

(6.34)

Ω

Z Z

=− b(U0N (0) + U ext (0) − uN (0))∂t U0N (0)dλ2y dx Ω ΓR Z Z − θ∂t U ext (0)∂t U0N (0) − D∆U ext (0)∂t U0N (0) Ω Ω Z Z N 2 |∂t u (0)| + |∂t v N (0)|2 Ω×Y Ω×Y Z Z N N + d1 ∇y u (0)∇y ∂t u (0) + d2 ∇y v N (0)∇y ∂t v N (0) Ω×Y Ω×Y Z Z = b U0N (0) + U ext (0) − uN (0) ∂t uN (0)dλ2y dx Ω ΓR Z +k η(uN (0), v N (0))(∂t uN (0) + α∂t v N (0))

(6.35)

Ω×Y

Integrating by parts in the higher order terms, and using the fact that uN satisfies


15

the transmission condition (2.7) in a weak sens, we obtain Z

Z

Z |∂t uN (0)|2 + |∂t v N (0)|2 (6.36) Ω Ω×Y Ω×Y Z Z Z = D∆U0N (0)∂t U0N (0) + d1 ∆y uN (0)∂t uN (0) + d2 ∆y v N (0)∂t v N (0) Ω Ω×Y Ω×Y Z Z − b(U0N (0) + U ext (0) − uN (0))∂t U0N (0)dλ2y dx Ω ΓR Z Z ext N − θ∂t U (0)∂t U0 (0) − D∆U ext (0)∂t U0N (0) Ω Ω Z +k η(uN (0), v N (0))(∂t uN (0) + α∂t v N (0)) θ|∂t U0N (0)|2 +

Ω×Y

Now, the regularity properties of the initial and boundary data, together with the stability of the projection operators PxN and PyN with respect to the H 2 -norm, yield the desired bounds on the time derivatives at t = 0: Z

|∂t U0N (0)|2

Ω

Z

N

2

Z

|∂t u (0)| +

+ Ω×Y

|∂t v N (0)|2 ≤ c.

(6.37)

Ω×Y

Inserting now (6.37) into (6.33), and using Gronwall’s inequality, the estimates of the Theorem are proved. The estimates proved in Theorem 6.1 still don’t provide the compactness for the solutions (U0N , uN , v N ) needed to pass to the limit N → ∞ in the nonlinear terms of the variational formulation (6.7)–(6.9). In the next theorem, additional regularity of the solutions uN , v N with respect to the macroscopic variable x is proved. Theorem 6.2. ||∇x uN ||L∞ (S,L2 (Ω×Y ) + ||∇x v N ||L∞ (S,L2 (Ω×Y ) ≤ c N

N

||∇y ∇x u ||L2 (S,L2 (Ω×Y ) + ||∇y ∇x v ||L2 (S,L2 (Ω×Y ) ≤ c

(6.38) (6.39)

Proof. Let Ω0 be an arbitrary compact subset of Ω, and let h ∈]0, dist(Ω0 , ∂Ω)[. Denote by Uhi , uih , and vhi the difference quotients with respect to the variable xi , for i = 1, . . . , n, of U N , uN , and v N respectively. For example,

uih (t, x, y) :=

uN (t, x + hei , y) − uN (t, x, y) , h

for all t ∈ S, x ∈ Ω0 , and y ∈ Y . We consider the following variational formulation

16


satisfied by these difference quotients, tested with (ϕ, φ, ψ) := (Uhi , uih , vhi ): Z Z d i 2 θ|Uh | + D |∇Uhi |2 = (6.40) dt Ω Ω Z Z 1 b(U0N (t, x + hei ) + U ext (x + hei ) − uN (t, x + hei , y)) − h Ω ΓR − b(U0N (t, x) + U ext (x) − uN (t, x, y)) UhN dλ2y dx, Z Z d |uih |2 + d1 |∇y uih |2 = (6.41) dt Ω×Y Ω×Y Z Z 1 b(U0N (t, x + hei ) + U ext (x + hei ) − uN (t, x + hei , y)) h Ω ΓR − b(U0N (t, x) + U ext (x) − uN (t, x, y)) uih dλ2y dx Z 1 N +k η(u (t, x + hei , y), v N (t, x + hei , y)) − η(uN (t, x, y), v N (t, x, y)) uih , Ω×Y h Z Z d |vhi |2 + d2 |∇y vhi |2 = (6.42) dt Ω×Y Ω×Y Z 1 N αk η(u (t, x + hei , y), v N (t, x + hei , y)) − η(uN (t, x, y), v N (t, x, y)) vhi . Ω×Y h Firstly, we add equations (6.40) and (6.41), and estimate the coupling terms using the Lipschitz property of b, the regularity of U ext , and the interpolation-trace estimate (5.3) as follows. Z Z 1 b(U0N (t, x + hei ) + U ext (x + hei ) − uN (t, x + hei , y)) (6.43) − h Ω ΓR 2 − b(U0N (t, x) + U ext (x) − uN (t, x, y)) (Uhi − uN i )dλy dx Z Z ≤ cˆ |Uhi − uih |2 + |Uhext,i ||Uhi − uih |dλ2y Ω Γ Z Z R |Uhi |2 + |Uhext,i |2 + |uih |2 dλ2y dx ≤c Ω ΓR Z Z Z Z Z i 2 i 2 ≤c 1+ |Uh | + C() |uh | + |∇y uih |2 . Ω

Ω

Y

Ω

Y

Next, we estimate the reaction term in the equation for uih by using the Lipschitz properties of R and Q and the uniform boundedness of uN and v N as follows. Z k R(uN (t, x + hei , y))Q(v N (t, x + hei , y)) − R(uN (t, x, y))Q(v N (t, x, y)) uih Ω×Y h Z R(uN (t, x + hei , y)) − R(uN (t, x, y)) = Q(v N (t, x + hei , y)) h Ω×Y Q(v N (t, x + hei , y)) − Q(v N (t, x, y)) i N + R(u (t, x, y)) uh h i ≤ Qmax cR |uih |2 + Rmax cQ uN h vh ≤ 2 (Qmax cR + Rmax cQ ) (|uih |2 + |vhi |2 ).

(6.44)

The reaction term in the equation for vhi can be estimated analogously. Finally,

17


choosing :=

d1 2

in (6.43), and summarizing the above estimates, we get the inequality Z Z Z d d d i 2 i 2 |U | + |u | + |v i |2 dt Ω h dt Ω×Y h dt Ω×Y h Z Z Z + |∇Uhi |2 + |∇y uih |2 + |∇y vhi |2 Ω Ω×Y Ω×Y Z Z Z i 2 i 2 i 2 ≤c 1+ |Uh | + |uh | + |vh | (6.45) Ω

Ω×Y

Ω×Y

Integrating with respect to time in (6.45), and applying Gronwall’s inequality yields for all i = 1, . . . , n, and N ∈ N the estimates ||uih ||L∞ (S,L2 (Ω×Y ) + ||vhi ||L∞ (S,L2 (Ω×Y ) ≤ c ||∇y uih ||L2 (S,L2 (Ω×Y )

+

||∇y vhi ||L2 (S,L2 (Ω×Y )

≤ c.

(6.46) (6.47)

with a constant c independent on i, h, and N . Now applying the result on difference quotients from Lemma 7.24, in [6], the estimates (6.38)-(6.39) follow. 6.3. Convergence of the Galerkin approximates. In this section, we prove the convergence of the Galerkin approximations (U0N , uN , v N ) to the weak solution of the two-scale problem (2.1)–(2.10). Based on the estimates proved in Section 6.2, we first derive the following convergence properties of the sequence of finite-dimensional approximations. Theorem 6.3. There exists a subsequence, again denoted by (U0N , uN , v N ), and a 2 limit (U0 , u, v) ∈ L2 (S; H 1 (Ω)) × L2 (S; L2 (Ω; H 1 (Y ))) , with (∂t U0N , ∂t uN , ∂t v N ) ∈ 2 L2 (S × Ω) × L2 (S × Ω × Y ) , such that 2 (U0N , uN , v N ) → (U0 , u, v) weakly in L2 (S; H 1 (Ω)) × L2 (S; L2 (Ω; H 1 (Y ))) (6.48) (∂t U0N , ∂t uN , ∂t v N ) → (∂t U0 , ∂t u, ∂t v) weakly in L2 (U0N , uN , v N ) → (U0 , u, v) strongly in L2 uN |ΓR → u|ΓR strongly in L2 (S × Ω, L2 (ΓR ))

(6.49) (6.50) (6.51)

Proof. The estimates from Theorem 6.1 immediately imply (6.48) and (6.49). Since ||U0N ||L2 (S,H 1 (Ω)) + ||∂t U0N ||L2 (S,L2 (Ω)) ≤ c, Lions-Aubin’s compactness theorem, see [9], Theorem 1, page 58, implies that there exists a subset (again denoted by U0N ) such that U0N −→ U0 strongly in L2 (S × Ω). To get the strong convergences for the cell solutions uN , v N , we need the higher regularity with respect to the variable x, proved in Theorem 6.2. We remark that the estimates (6.38)-(6.39) imply that ||uN ||H 1 (Ω,H 1 (Y )) + ||v N ||H 1 (Ω,H 1 (Y )) ≤ c. Moreover, from Theorem 6.1, we have that ||∂t uN ||L2 (S×Ω×Y ) + ||∂t v N ||L2 (S×Ω×Y ) ≤ c.

18


Since the embedding H 1 (Ω, H 1 (Y )) ,→ L2 (Ω, H β (Y )) is compact for all 12 < β < 1, it follows again from Lions-Aubin’s compactness theorem that there exist subsequences (again denoted uN , v N ), such that (uN , v N ) −→ (u, v) strongly in L2 (S × L2 (Ω, H β (Y )), for all

1 2

(6.52)

< β < 1. Now, (6.52) together with the continuity of the trace operator H β (Y ) ,→ L2 (ΓR ),

for

1