Upscaling of reactive flows in domains with moving oscillating boundaries

We consider the flow and transport of chemically reactive substances (precursors) in a channel over substrates having complex geometry. In particular, these substrates are in the form of trenches forming oscillating boundaries. The precursors react at the boundaries and get deposited. The deposited layers lead to changes in the geometry and are explicitly taken into account. Consequently, the system forms a free boundary problem. Using formal asymptotic techniques, we obtain the upscaled equations for the system where these equations are defined on a domain with flat boundaries. This provides a huge gain in computational time. Numerical experiments show the effectiveness of the upscaling process.


Introduction
This work deals with the upscaling of processes which are defined in a domain having oscillating (rough) boundaries. The aim is to approximate the model by one defined in a simpler domain with flat (non-oscillating) boundaries. The work presented here is motivated by reactive flow models where the reactive substances are transported by a carrier fluid (gas or liquid).
The reactions take place at the rough boundaries, and the reaction products are immobile species which get attached to or detached from the surface. Examples of these processes are crystal precipitation and dissolution processes [21] and chemical vapor deposition processes [10]. The formation of reaction products at the boundary leads to changes in the thickness of the solid layer with time. There are two options: one is to ignore the geometry change and take a fixed geometry model and describe the deposition/dissolution process by an equation defined at the boundary and the other option is to take the changes in the geometry explicitly into account in which the thickness of the solid layer provides information about the amount of deposition/dissolution. We follow the second approach here. This implies that the geometry may change with time and since the change in the geometry is unknown, we have a free (evolving) boundary problem. The problem is to find an upscaled equation defined in a simple domain with flat boundaries for this free boundary problem posed in a complex domain with oscillatory boundaries.
A motivation for this work is the process of making 3-D batteries through the process of chemical vapor deposition (CVD) where a part of the domain consists of trenches (rough surface) instead of being flat to enhance the surface area hence increasing the capacity of the battery [22,23,24]. The sizes of these trenches may vary and are in the order of micrometers (typical size ∼ 10 µm) and the substrate size is in the order of centimeters (typical size ∼ 30 cm). In the trenches, the thickness of different layers of deposition (of electrodes/anode/cathode) may not be negligible compared to the size of trenches. In this situation, it is important to take into account the thickness of the deposited layer and consider a model for describing the deposition process for the time-dependent geometry.
We use upscaling arguments, in particular a formal (matched) asymptotic approach, to provide a rational derivation of this approximation. The presence of free boundary makes the rigorous arguments difficult and here we restrict to formal asymptotic arguments. One important aspect of this upscaling approach is that the system of equations defined inside the domain remains unchanged and only the boundary conditions are modified. Naturally, such an upscaling leads to computational efficiency as the fine scales need not be resolved. Further, for small enough oscillations, one takes advantage of the scale separation and obtains better approximation. Hence, the finer these oscillations are, the better approximation can be achieved. Also, the upscaled domain is fixed in time and the movement of the boundary of the original domain leads to a time-dependent boundary condition for the upscaled equations.
In this work we provide the upscaling for both the transport and the flow problems defined in a domain with moving and oscillating boundaries. As we apply only formal arguments, non-linear reactions rates have been considered. Further, in the work presented here, for the sake of exposition, we have considered only one species, however, an extension to multi-species and multi-component is straightforward. Also, more general diffusion terms may be dealt with the similar techniques considered here.
There exists extensive amounts of work dealing with both the formal and the rigorous upscaling of rough boundaries. Most of the work however refers to the situation when the geometry remains fixed. Relevant to the work presented here are [10] for the transport problem and [13,14,15] for flow problem respectively. In the former, the formal approach is used for diffusion equation while the latter papers deal with rigorous treatment to the upscaling of Stokes equation in the fixed geometry situation. The difference with [10] is that the cited work deals with the fixed geometry situation and also the flow is absent. Although our techniques are inspired from [10], however, there are some differences in the way of our derivation. In an another related work, [8], rigorous treatment for the upscaling in the absence of flow and linear reaction rates when the geometry changes are explicitly taken into account. Our results are consistent with the results in the mentioned literature. Because of the formal approach employed here, we are able to take into account both the flow and the transport problem together and also much more general reaction rates are considered. For other related work on upscaling of reactive flows/transport involving oscillating boundaries we refer to [1,2,3,5,9,16,19]. These kind of problems find applications in making 3-D solid state batteries, flow in rough pipes, semiconductor industries, bio-film growth in porous medium and describing the evaporation processes in plants or through soil [4,7,11,18,20,25,26,27].
The article is organized as follows. We start with the description of the geometry and the formulation of the problem in Section 2. In Sections 3 and 4 we treat the flow problem and transport problem. In Section 5, we provide numerical computations where we compare the upscaled equations with the original equations.

Modeling equations
We describe the setting of the problem and start with the precise description of the geometry considered here. To define the model, we consider the flow carrying reactive substances (precursors). The flow passes over the substrate and the precursor is deposited on the substrate making it a transport process with reactions at the boundary. To take into account the thickness of the deposited layer, we provide a model that describes the deposition process for the time-dependent geometry. The model consists of two components, one the flow component and the other the transport component. For the flow, we take for simplicity of computations, the Stokes equation and the transport part is described by a linear convection-diffusion equation. The boundary conditions for both the transport and flow equations are posed on a time-dependent boundary and we also prescribe a law for describing the movement of the geometry. Figure 1: Schematic for the domain Ω ε : channel with rough boundary Γ ε . Initially, the boundary Γ ε ⊂ ∂Ω ε consists of a periodic function εh( x ε ). Note that the geometry may change in time due to reactions taking place.

Geometry
Let Ω ε denote the domain having rough boundaries denoted by Γ ε . We define a system of equations describing the flow and transport processes in Ω ε with reactive boundary conditions defined on Γ ε . The question that we are concerned with is defining the appropriate system of equations in the simpler domain Ω with flat (non-oscillating) boundary such that the solutions in a simpler domain Ω approximates those in Ω ε . We proceed to give the precise definitions of these geometries. A schematic of the geometry is presented in Figure 1.
First, the homogenized domain with flat boundaries is denoted by Ω := (0, 1) 2 . Let h be a given smooth function that is 1periodic. Let us define h ε (x, t) using the periodicity of h such that h ε (x, 0) := εh( x ε ).
We use the above to define Ω ε where h ε (x, 0) is as defined above and hence is a smooth function. Let ∂Ω ε denote the boundary of Ω ε . The oscillating boundary Γ ε (⊂ ∂Ω ε ) is defined as: We remark that the h ε is not a priori known except for t = 0 and depends on the unknown solution of the transport equation. Hence, this defines a free-boundary problem. The inlet and the outlet boundaries are correspondingly defined as This particular scaling for the boundary ensures that the Lebesgue measure of Γ ε remains bounded and of order O(1). Notice that we emphasize the time-dependence of the definition of Ω ε , Γ ε as opposed to the homogenized domain Ω which is time-independent.

The model
We start with the modeling for the processes taking place in the geometry described above. For the sake of simplicity, we write all the equations in dimensionless form (see [21] for the dimensional model and nondimensionalization). As described above, because of the reactions taking place at the boundaries, the thickness of the layer attached to the boundary may change with time. The thickness of this layer is h ε (x, t) − h ε (x, 0). The growth of the layers lead to variation in the domain, and implicitly to a change in the flow, an effect that is taken into account.
For the reactants in the domain, the different processes are diffusion, transport by the fluid flow, and reactions taking place at the boundaries of the domain. The flow problem is modeled by the Stokes equations, where p ε is the pressure field, q ε is the flow field and µ is the dynamic viscosity. At the inlet, we take a parabolic velocity profile normal to the inlet, where Q > 0 is a normalisation factor related to the total flow in the x-direction and e 1 is the unit vector along the x-axis. For the outlet we prescribe the pressure p ε = 0. At other boundaries, we assume no-flow boundary conditions, in particular q ε = 0 at Γ ε .
For the transport equation, for the sake of presentation, we consider only one species here, with its concentration denoted by u ε . The convection-diffusion equation and the boundary conditions describing the transport process of the reactant are where u b is a given non-negative concentration. The diffusion coefficient is taken as 1. For simplicity, we assume homogeneous Neumann boundary conditions on the remaining boundary, that is, on The reactions take place only at the interface between the fluid and the solid, that is, Γ ε (t). The corresponding mathematical model involves the outer normal ν ε (t) to Γ ε (t), At Γ ε (t), mass conservation yields where ρ is the molar density of the reactant in the solid phase (the adsorbed substance, or the deposited layer) and v n is the outward normal velocity of the interface. Furthermore, the normal velocity of the interface v n is proportional to the reaction rate Here f is a given function assumed to be sufficiently smooth. Similar results can be derived formally even for non-Lipschitz rates like Freundlich isotherms or multi-valued reaction rates as encountered for dissolution processes [16].
Then the normal velocity v n becomes and (6) transforms into Further, we assume ρ to be of O(1/ε) so that the movement in the boundary is of order ε thereby retaining the roughness as the reactions change the geometry.

The flow problem
We use here formal matched asymptotics to derive the upscaled equations. The procedure consists of considering the solutions in different regions, one close to Γ ε (the inner solution) and far from it (the outer solution). We later derive the matching conditions at an intermediate scale.
Using the matching conditions, we obtain the upscaled boundary conditions.

The inner region
Close to Γ ε we assume q ε = q ε (x, y, t), p ε = p ε (x, y, t) and for the expansions we assume In the expansions we assume that the thickness of the deposition layer remains of the order O(ε) to retain roughness. Define: ξ := x ε , η := y ε . Following the periodicity in the geometry of oscillations, we assume periodicity of the inner variables in the ξ co-ordinate. This implies Denoting by e 1 and e 2 the unit vectors in the xand y-directions respectively, we start with the gradient terms.
This implies for i = 1, 2 Similarly, We proceed by inserting these expansions in the equations. For the continuity equation, we have Equating the terms of similar order, for ε −1 terms we get and for ε 0 terms we have We proceed with the Stokes equation. Plugging in the expansions for ∆q ε and ∇p ε , we obtain after cascading the different orders of ε With the expansion for h ε we have for the leading order term for the normal of Γ ε (t) For the boundary conditions, the no-slip boundary condition leads to Since there is no forcing term combined with the periodicity in ξ and homogeneous Dirichlet data at a part of the boundary, we conclude q 0 = 0.
For the stokes equation, we first write down separately ∆Q ε , ∇P ε for convenience. For i = 1, 2, we have which upon expansion and rearrangement gives Similarly, Plugging these expansions into the Stokes equation and separating the different orders of ε we obtain For the e 2 direction, for the Stokes equation the different terms of different orders of ε provide From (26) and the periodicity in ξ we obtain that Q 0 is independent of ξ and using this in (27) we conclude 1 (x, y).
Integrating (25) over ξ and using periodicity result in We integrate (28) over ξ from 0 to 1 and use the fact that Q (2) 0 and P 0 are independent of ξ and we get Note that Q are ξ-independent. This means that they satisfy the original Stokes model.

Matching conditions
Here we couple the inner solutions and the outer solutions. We use the asymptotic expansion, for the flow vector, and Taylor's expansion at y = εη to obtain For the inner solution we have the asymptotic expansion Therefore, by matching, we conclude and further using the above, lim η→∞ ∂ η q 1 (x, ξ, η) − ∂ y Q 0 = 0.
Now we are in a position to provide the leading order approximation and the first order approximation. For the leading order approximation Q 0 we simply have Q 0 = 0 at y = 0.
Remark 3.1 To obtain the next order approximation, let us define Q e := Q 0 +εQ 1 . From the matching condition we have This provides us the effective boundary condition Q e = εα(t)∂ y Q e at y = 0.
Notice that the factor α depends on the geometry which changes with time. However, here time only acts as a parameter and the free boundary does not appear explicitly. The computation requires the information about the geometry and we need to provide the law for the movement of the boundary of the cell, which we describe in the next section.

The transport equation
We consider the convection-diffusion equation (3) with boundary conditions (5), (6) and the movement of the boundary defined in (8). As before, we decompose the solution in the two regions, one close to Γ ε (the inner solution) and away from it (the outer solution). Denote the outer and inner solutions respectively by u ε and U ε and assume for the variables • Inner solution (close to Γ ε :) u ε := u ε (x, x ε , y ε , t) • Outer solution (away from Γ ε ): U ε := U (x, x ε , y, t).
We make the following asymptotic ansatz:

The outer solution
We start with the outer solution. For the different terms in (3) we have (as before with ξ = x/ε; η = y/ε) Inserting the asymptotic expansion in the transport equation (3), we obtain for the lowest order term, that is, the ε −2 term: ∂ ξξ U 0 = 0 which combined with periodicity in ξ implies, The ε 0 term provides

The inner solution
Using the above asymptotic ansatz, we obtain Using the above asymptotic expansions in the equation (3) for the inner region and rearranging the terms and using q 0 = 0 We separate the terms according to the different orders of ε; for the ε −2 term, we have Proceeding further with ε −1 terms, and the corresponding ε 0 terms provide For the boundary condition (5) we get −ν 0 · ∇ ξ,η u 0 = 0, and Using (35), the boundary condition (38), the periodicity in ξ and that U 0 does not depend on ξ , we conclude For the leading order term of the normal of the oscillating boundary we have Notice that the normal is dependent on time because of the change in geometry being taken into account. For any time t we denote by Γ L and Γ(T ) respectively the lower and upper parts of the boundaries of the inner region, Integrating (36) over the inner region, since by the periodicity with respect to ξ, the contributions to the lateral terms cancel (corresponding to the sides ξ = 0 and ξ = 1) and we get We will see that the above equation provides the effective boundary condition.

The matching conditions
As before, we couple the inner solutions and the outer solutions by deriving the matching conditions. We have for U ε , the asymptotic expansion and Taylor's expansion at y = εη to obtain For the inner solution we have the asymptotic expansion Therefore, by matching, we conclude and further using (43) above, We use the matching conditions for deriving the effective boundary condition below.

The effective boundary condition
Using the matching condition (44), we obtain Use the boundary condition at Γ L to obtain The fact that u 0 = u 0 (x, t) and the matching condition (43) imply u 0 = U 0 . This leads to The equations (46) and (47) together with (41) provide us the effective boundary condition Note that Γ L is a function of time t and using parametrization of the boundary, we have For the movement of the boundary, using (8) for the leading order term h 0 and using the assumption that ρ is of order O(1/ε), we have Recalling that Ω = (0, 1) 2 and the definition of Γ L in (40), the upscaled equations become,

Numerical simulations
To study the approximation of the upscaled equations to the original equations, we make the following choices for the geometry. Initially, Γ ε (t = 0) is defined as follows: and we conduct the numerical experiments for different ε. We choose a linear rate for the reaction, specifically, f (u, h) = ku, and choose ρ = 50 so that 1 ρ is of order ε. This ensures that the deposition profile retains the roughness of the boundary. Furthermore, for the equations, we choose the following parameters: For the computations for (1)-(3), we use finite element method with BDF time stepping for solving the equations. Since, the geometry changes with time, we use the Arbitrary Lagrangian Eulerian (ALE) method to solve the problem on a moving domain. This method is a generalization of Eulerian and Lagrangian descriptions of the free boundaries. We use ALE method as implemented in the COMSOL Multiphysics package [12], with Laplacian smoothing [6]. For more details on ALE method we refer to [6].
For the upscaled equations (50) with the boundary condition (52), we use finite difference with implicit in time for the diffusion and convection term and explicit time stepping for the reaction term. Note that for the upscaled geometry, we have an explicit solution for Stokes equation (48),(49) under the given conditions, namely, the parabolic inlet profile and no-slip boundary conditions on the lateral flat boundaries (51); Also, for (53) we use time-explicit discretization. For the sake of convenience, as the geometry is simple, we implement the scheme in Matlab [17].
For the numerical simulations, we consider the following situations.

Concentration at the boundary
We compute the full solution for ε = 0.1, 0.04, 0.02 and plot u ε and u at the boundary Γ ε and y = 0 boundary of Ω respectively at t = 0.5. The plot is shown in Figure 2. Due to the oscillations in the boundary, we have the boundary layer and the maximum error takes place at the boundary itself. We compute the concentration at the oscillating boundary for given ε and then plot it against the upscaled concentration u at y = 0. Because of the oscillating boundary, the concentration has an oscillating profile while the upscaled concentration has a monotonic profile; however, as ε decreases, u ε converges to u.
In Table 1 we give the values of the error as well as the convergence rate α computed from the preceding formula. The convergence rate is close to ε and this is expected as for linear reaction rates, H 1 convergence in Ω is of the order √ ε which suggests that we expect a better rate than √ ε for L 2 at the boundary.

Mass balance: full vs upscaled
Further, we compare the mass balance for the full model with the upscaled model. We definē Observe thatū provides information regarding the average concentration (and hence the flux) into one periodic unit for ε model. We compare this with the upscaled equation in Figure 3. We note that as ε decreases, the mass balance for the ε model tends to the upscaled one. The agreement is very good for small ε which indicates the quality of upscaling. Also, the mass balance suggests that in one periodic unit, the net amount of flux in the upscaled model approximates very well to the ε model which means that by defining a cell problem with this flux provides how to compute the deposition profile.

Deposition profile
The upscaled model also provides information about the deposition profile. The solution to (53) gives h 0 which after suitable rescaling is compared to h ε . The deposition profile for the upscaled equation approximates that of the ε model and the error is of ε order. The corresponding plot is shown in Figure 4.

Conclusions
We have used formal asymptotic arguments to derive the upscaled equations. The model defined in a geometry with oscillating boundary is approximated by a model defined in a simpler geometry with flat boundaries. Also, the geometry changes due to the reactions taking place at the oscillating boundary have been incorporated. In the upscaled model, the geometry changes lead to time-dependent boundary condition. For the flow problem, we derived upscaled equations for both the leading order and the next order approximations. For the leading order, the no-flow boundary condition at the oscillating boundary is approximated by homogeneous Dirichlet boundary condition; however, in the first order effect, the upscaled boundary condition is of slip-type with the slip parameter α implicitly dependent on time t because of the changes in geometry.
For the transport equation, the approximation model retains the equation inside the domain but with a modified boundary condition. This boundary condition takes into account the geometry of ε dependent model. For the numerical experiments, we compute the solutions for different ε and then compare the approximate solution. As ε decreases, we see that the upscaled model approximates the original model pretty well. For the concentration at the boundary, for the original problem, the profile is oscillating while for the upscaled model, the profile is monotonic. However, as ε decreases, the approximations for the concentrations at the boundary are excellent. Further, we have investigated the convergence of L 2 error for the concentrations at the boundary and we find that the order is ε. Also, the upscaling process provides a good approximation for the deposition profile.
The quality of upscaling as evidenced by the numerical experiments provides a convincing argument for using these upscaling techniques instead of solving the original problem in a complex geometry.