Crouzeix-Raviart MsFEM with Bubble Functions for Diffusion and Advection-Diffusion in Perforated Media

The adaptation of Crouzeix - Raviart finite element in the context of multiscale finite element method (MsFEM) is studied and implemented on diffusion and advection-diffusion problems in perforated media. It is known that the approximation of boundary condition on coarse element edges when computing the multiscale basis functions critically influences the eventual accuracy of any MsFEM approaches. The weakly enforced continuity of Crouzeix - Raviart function space across element edges leads to a natural boundary condition for the multiscale basis functions which relaxes the sensitivity of our method to complex patterns of perforations. Another ingredient to our method is the application of bubble functions which is shown to be instrumental in maintaining high accuracy amid dense perforations. Additionally, the application of penalization method makes it possible to avoid complex unstructured domain and allows extensive use of simpler Cartesian meshes.


Introduction
Many important problems in modern engineering context have multiple-scale solutions e.g., transport in truly heterogeneous media like composite materials or in perforated media, or turbulence in high Reynolds number flows are some of the examples. Complete numerical analysis of these problems are difficult simply because they exhaust computational resources. In recent years, the world sees the advent of computational architectures such as parallel and GPU programming; both are shown to be advantageous to tackle resource demanding problems. Nevertheless, the size of the discrete problems remains big. In some engineering contexts, it is sometimes sufficient to predict macroscopic properties of multiscale systems. Hence it is desirable to develop an efficient computational algorithm to solve multiscale problems without being confined to solving fine scale solutions. Several methods sprung from this purpose namely, Generalized finite element methods (Babuška et al., 2003), wavelet-based numerical homogenization method (Dorobantu and Engquist, 1998), variational multiscale method (Nolen et al., 2008),various methods derived from homogenization theory (Bourgeat, 1984), equationfree computations (Kevrekidis et al., 2003), heterogeneous multiscale method (Weinan and Engquist, 2003) and many others. In the context of diffusion in perforated media, some studies have been done both theoretically and numerically in (Cioranescu et al., 2006), (Cioranescu and Murat, 1997), (Henning and Ohlberger, 2009), (Hornung, 1997), and (Lions, 1980). For the case of advection-diffusion a method derived from heterogeneous multiscale method addressing oscillatory coefficients is studied in (Deng et al., 2009).
In this paper, we present the development of a dedicated solver for solving multiscale problems in perforated media most efficiently. We confine ourselves in dealing with only stationary diffusion and advection-diffusion problems as means to pave the way toward solving more complicated problems like Stokes. We begin by adapting the concept of multiscale finite element method (MsFEM) originally reported in (Hou and Wu, 1997). The MsFEM method relies on the expansion of the solution on special basis functions which are pre-calculated by means of local simulations on a fine mesh and which model the microstructure of the problem. By contrast to sub-grid modeling methodologies, the multiscale basis functions are calculated from the actual geometry of the domain and do not depend on an (often arbitrary) analytical model of the microstructure. A study on the application of MsFEM in porous media has been done in (Efendiev and Hou, 2007), and although it could have bold significance in geo-or biosciences, they can be applied also in different contexts, e.g., pollutant dispersion in urban area (Carballal Perdiz, 2011) or on similar problems which are extremely dependent on the geometry of perforations but their full account leads to very time consuming simulations. Textbook materials on the basics of MsFEM can be found in (Efendiev and Hou, 2009).
It is understood that when constructing the multiscale basis function, the treatments of boundary condition on coarse elements greatly influence the accuracy of the method of interest. For example, in the original work of Hou and Wu, the oversampling method was introduced to provide the best approximation of the boundary condition of the multiscale basis functions which is also of high importance when dealing with non-periodic perforations. Oversampling here means that the local problem in the course element are solved on a domain larger than the element itself, but only the interior information is communicated to the coarse scale equation. This reduces the effect of wrong boundary conditions and bad sampling sizes. The ways in which the sampled domain is extended lead to various oversampling methods, see (Efendiev and Hou, 2009), (Chu et al., 2008), (Henning and Peterseim, 2012), (Efendiev et al., 2013). The non-conforming nature of Crouzeix-Raviart element, see (Crouzeix and Raviart, 1973), is shown to provide great 'flexibility' especially when non-periodically perforated media is considered. In the construction of Crouzeix-Raviart multiscale basis functions, the conformity between coarse elements are not enforced in a strong sense, but rather in a weak sense i.e., the method requires merely the average of the "jump" of the function to vanish at coarse element edges. When very dense perforations are introduced, which often makes it virtually impossible to avoid intersections between coarse element edges and perforations, the benefit of using Crouzeix-Raviart MsFEM is significant for it allows the multiscale basis functions to have natural boundary conditions on element edges making it insensitive to complex patterns of perforations. Moreover, the integrated application of penalization method enables one to carry the simulations onto simple Cartesian meshes. Note that some methods derived from homogenization theory may provide robust and accurate results provided that the underlying multiscale structure or subgrid effects satisfies the necessary constraints which is not the case for problems with non-periodic perforations. In this paper several computational results with non-periodic perforations are given to highlight the feasibility of our method in such circumstances. Another important ingredient to our method is the multiscale finite element space enrichment with bubble functions. Again, when very dense perforations are considered, it is both crucial and difficult to capture correct approximations between perforations for which the application of bubble functions is offered as the remedy. We illustrate these problems in our paper to highlight the contribution of bubble function in improving the accuracy of our MsFEM. Our work continues the application of Crouzeix-Raviart MsFEM done on prototypical elliptic problems (Le Bris et al., 2013a) and on diffusion problems with homogeneous boundary condition (Le Bris et al., 2013b). Improvements are done to the earlier work by introducing bubble functions and to the latter by extending the application to advection-diffusion problems with non-homogeneous boundary conditions. The paper is organized as the following. In chapter 3 we begun with the formulation of the problem and the construction of our MsFEM. Here we explain the construction of Crouzeix-Raviart MsFEM functions space with bubble functions and the multiscale basis functions. In chapter 4 the application of non-homogeneous boundary conditions is explained. In chapter 5 we describe the application of penalization method. Demonstrations of our MsFEM in terms of computational simulations and its analysis can be found in chapter 6 followed by some concluding remarks.

Multiscale Methodology
In our paper, we tackle our problems with the non-classical Multiscale Finite Element Method (MsFEM) developed by Hou and Wu (1997). In this section we revisit the general framework in which a class of multiscale methods is commonly described. In general, one is interested in solving the following problem, namely in which X is a functional space with very high dimension m such that m= For instance, X could be the finite element space on a very fine mesh that captures smallest characteristic length in d spatial dimension.
The maxim of multiscale methods is to reduce the difficulties in dealing with a large set of unknowns inherent to Eq. 2.1 by transforming it into smaller dimensional problems. In doing so, one builds a functional space X 0 with smaller dimension than X, the elements of which constitutes an "optimal" reduced representation of the local solutions of problem 2.1. The dimension of X 0 , m 0 is ideally very small compared to m. To link X and X 0 one could construct the projection operator P : X → X 0 , and the reconstruction operator R : X 0 → X satisfying P R = Id, but RP = Id. A multiscale method renders the original problem 2.1 into smaller sized problem In the context of MsFEM, this method consists of two principal steps namely (i) determination of a basis that spans space X 0 , and (ii) the Calculations of 2.2. The reconstruction operator is used in the end to recover the multiscale approximation of u, that is v = Ru 0 . The accuracy of such approximation depends heavily on step (i). Several techniques can be used to construct the X 0 basis, and in terms of MsFEM, the finite element method is used to compute such "multiscale" basis function. Then, the same finite element method is used to compute the coarse scale solutions.

Basic Formulation
We consider an advection-diffusion problem laid in a bounded domain Ω ∈ R d within which a set B of perforations is included. From here on we assume that the ambient dimension is d = 2. The perforated domain with voids left by perforations is denoted Ω =Ω\B , where denotes the minimum width of perforations. The advection-diffusion problem is then to find u : Ω → R which is the solution to where f : Ω → R is a given function, g is a function fixed on boundary ∂Ω and w is a given velocity field. In this paper, we consider only the Dirichlet boundary condition on ∂B namely u |∂B = 0 thereby assuming that the perforation is opaque. Other kinds of tion for MsFEM basis function is considered, it is difficult to approximate correct coarse node solutions when one or more perforations coincide with any of the coarse element's boundaries. The approximation of the MsFEM basis function will be distorted and the whole solution will be affected. This problem often and can be relaxed by using oversampling methods. However, in practice, this brings an inherent inconvenience since the size and position of perforations are most of the times unpredictable which requires some problem dependent parameters to be introduced. Moreover, the computational cost also increases due to enlarged sampled domain. The Crouzeix-Raviart basis functions are non-conforming throughout the computational domain. The continuity of the functions are enforced weakly i.e., it requires no fixed values across the boundaries but rather vanishing "jump" averages on each edge. In order to explain the MsFEM space in the vein of Crouzeix-Raviart's finite element, we define a mesh T H in Ω which are discrete polygons with each diameter at most H and made up by n H coarse elements and n e coarse element edges. Denote H the set of all edges e of T H which includes edges on the domain boundary ∂Ω. It is assumed that the mesh does not include any hanging nodes and each edge is shared by two elements except those on ∂Ω which belongs only to one element. T H is assumed a regular mesh. By regular mesh, we mean for any mesh element T ∈ T H , there exists a smooth one-to-one mapping M :T → T whereT ⊂ R d is the element of reference, and that ∇M L ∞ ≤ DH, We introduce the functional space for Crouzeix-Raviart type MsFEM with bubble function enrichment where [[u]] denotes the jump of u over an edge. The MsFEM approximation to Eq.

Calculation of Crouzeix-Raviart and Bubble functions
The basis for V H contains functions associated to edges e and mesh elements T ∩(Ω\ B )∩Ω . The former has the notation Φ e and the latter Φ B . The edges composing T k are denoted Γ k i with i = 1,··· , N Γ whereas k = 1,··· ,n H . The Crouzeix-Raviart multiscale basis functions Φ k e i are then the unique solution in where λ k i are the Lagrange multipliers. In a weak form, the finite element approximation of the above problem is to find In this paper, Eq. 3.7 is solved with standard Galerkin approximation of Q1 FEM within a coarse elements T k . We choose a finite fine basis function {Ψ m } n f m=1 where n f is the number of nodes in a coarse element T k such that Φ e k = ∑ n f m=1 Φ em Ψ m . Hence within T k Eq. 3.7 can be discretized and written in matrix form as (3.10) and p = 1,2,··· , N Γ . The bubble function Φ k B can be obtained by solving for each element T k In a weak form, the finite element approximation of the above problem is to find The bubble functions Φ k B are also computed with standard Galerkin approximation of Q1 FEM within T k by choosing fine scale basis function Ψ such that Φ k B = ∑ n f m=1 Φ Bm Ψ m . Having acquired the Crouzeix-Raviart multiscale basis functions Φ k e and the bubble functions Φ k B , the coarse calculation Eq. 3.3 can be done with Galerkin approximation such that the approximated solution u H is described as y). (3.13) Here, u i are finite element approximation of u corresponding to coarse element edges i ∈ H whereas u k are finite element approximation of u corresponding to coarse element T k ∈T H . With this, the reconstruction of the quantity u H on the fine mesh can also be done. In this paper, although the general formulation is focused on advection-diffusion problem, various tests on diffusion-only cases are presented for showing the contributions of bubble functions and Crouzeix-Raviart multiscale basis functions. The MsFEM formulation for diffusion cases is largely similar to the advection-diffusion counterpart. Regarding the advection-diffusion problems, one could consider another approach for selecting the space for the basis functions namely by following the Petrov-Galerkin formulation which could be useful for problems with very high Péclet numbers, nevertheless this is not the emphasis of this paper. Rigorous studies on Crouzeix-Raviart MsFEM's numerical analysis and error estimates for highly oscillatory elliptic problems and for diffusion problems in perforated media can be found in (

Boundary condition
We propose to approximate the non-homogeneous dirichlet boundary condition in Eq.

Application of penalization method
Solving Eq. (3.1) in Ω as it is often requires complex and ad-hoc grid generation methods. For highly non-periodic perforations, complicated unstructured mesh is likely what engineers would resort to. In order to confine our computations in a simple uniform Cartesian domain, we incorporate the penalization method to solve Eq. (3.1). Henceforth, we solve instead the following Here h is the width of a fine scale element used to capture highly oscillatory basis functions. We introduce the penalization coefficient σ β which forces the solution to vanish 6 Numerical results

Application of Bubble Functions
In this paper, we first give numerical examples that would exhibit a case with very dense presence of perforations throughout the domain, so as to highlight the contribution of bubble function enrichment to the basis function set. In the first example, we applied bubble enrichment to a standard linearly boundary-conditioned, nodal based MsFEM without oversampling. We set aside the application of Crouzeix  that the one without bubble function enrichment fails to exhibit the correct solution at the interiors of coarse mesh whereas the solution with bubble function enrichment exhibits more consistency with that of the reference. In Fig. 3, we plot the standard MsFEM basis function alongside a bubble function used in this test. It clearly illustrates that with the presence of perforations this dense, the contribution of a standard MsFEM basis function inside the coarse element is insignificant. The bubble function applied in the coarse element is shown to contribute greatly to the approximation of the solution.
In Fig. 4, we plot the relative L2 errors with respect to the size of coarse element H of the standard MsFEM with or without bubble function. We notice that the one with bubble function enrichment gives a decreasing relative error when H increases away from . This is of course not the behaviour exhibited by the more standard MsFEM. However, we notice that both methods shows increasing error the moment H is reduced to be lower that . This is due to the fact that at this region, the edges of coarse mesh start to coincide with the perforations causing incorrect solution at the interior of the MsFEM basis function. Oversampling methods had been implemented to overcome this problem

Application of Crouzeix-Raviart MsFEM
In this section, we test the Crouzeix-Raviart MsFEM with bubble functions and compare it with the standard linearly-boundary-conditioned MsFEM also with bubble functions. The test is designed to analyse the sensitivity of the methods subject to placement of perforations. The computational domain remains the same with that of the previous section. The size of each perforation is now set as = 0.025.
The methods underwent two tests: In the first test, the arrangement of the perforations is made such that none of the coarse mesh edges coincides with the them. We call this test the non-intersecting case. In the second test, the allocation of these perforations is shifted both in x and y direction until all coarse element edges coincide with perforations. We call this test the intersecting case. In this example, we tried to illustrate a possible worst case scenario where each and every element edges coincide with perforations at three different locations. In both cases, we implemented 8×8 coarse elements each consists of 128×128 fine elements. Within each of these coarse elements, the Crouzeix-Raviart multiscale basis function and the bubble function are calculated on the fine mesh by solving Eqs. 3.6 and 3.11. The reference solution is calculated using standard Q1 FEM on 1024×1024 elements.
First, in Figs. 5, the results of these two methods for non-intersecting case are compared with the reference solution. The results shows quantitatively good accuracies displayed by both methods. The Crouzeix-Raviart MsFEM with bubble functions records 0.11407 L2 relative error whereas the standard MsFEM with bubble function records 0.11738. However, in the second test, where all coarse element edges coincide with per-  To get a better understanding on why the two methods exhibit such different accuracies, we plot the basis functions of the Crouzeix-Raviart and the standard MsFEM in Figs. 7(a) and (b). Here one can see that the Crouzeix-Raviart basis function cope very well with perforations on the cell edges and provide natural boundary conditions around them without violating the applied constraints. By contrast, the basis function of the standard MsFEM with linear boundary condition fails to give a correct approximation in the penalized region. Again we note that several methods including oversampling methods have been introduced as remedies to this kinds of problem on standard MsFEM. Nevertheless, Crouzeix-Raviart MsFEM also has the benefit of not increasing the size of the sampled domain for constructing the MsFEM basis functions. Moreover, the exhibited natural boundary condition gives a good deal of flexibility in tackling delicate cases for it is prohibitively difficult to avoid intersections between perforations and coarse element boundaries especially when simple Cartesian mesh is implemented. In the later examples, the applicability of our method on non-periodic pattern of perforations is demonstrated. For the case of diffusion with homogeneous Dirichlet boundary condition, one can refer to the previous works in (Le Bris et al., 2013b) where detailed comparison of performances between Crouzeix-Raviart MsFEM and other types of MsFEM including those with oversampling methods can be found. In this paper, more detailed study on the convergence behaviour of our method will be emphasized more on advection-diffusion case.

Advection-diffusion Problems on Perforated Domain
In this section we test our method on advection-diffusion problems on perforated domain with homogeneous boundary conditions. We implement Crouzeix-Raviart MsFEM with bubble function enrichment simply in the context of standard Galerkin approximation without any stabilizations. We reuse the computational set up done in 6.2 but with different source terms. The vector field w = (2y(1−x 2 ),−2x(1−y 2 )) set in a domain Ω = [−1,1] 2 . w determines a recirculating flow with streamlines {(x,y)|(1−x 2 )(1−y 2 ) = constant}. The source term f (x,y) is defined as follows  terns of perforations without having to resort to some perforation-dependent parameters nor to enlarge the sampled domain.

Non-periodically perforated domain
In this section, we test the applicability of our method on domain with non-periodic perforations. We consider two kinds of patterns of perforation as can be seen in figures 9. The Crouzeix-Raviart MsFEM with bubble function enrichment is implemented on Ω = [−1,1] 2 domain. The first case (case a) includes 400 perforations each with width of = 0.025 whereas on the second case (case b) we include 3600 perforations each with width of =0.005. We reuse the vector field w=(2y(1−x 2 ),−2x(1−y 2 )) which determines streamlines {(x,y)|(1−x 2 )(1−y 2 ) = constant}. The source term (6.1) is also applied. On both of these cases, we utilize a diffusion coefficient of A = 0.03. The result of the convergence tests done on these two cases is given on table 1. In figures (10) the contours of u solved on 8×8,16×16,32×32,64×64, and 128×128 elements are given alongside the reference solution on 1024×1024 solved using standard Q1-Q1 FEM. Although on 124×124 elements, the method already returns quite an identical result in comparison to the reference, the result solved on 32×32 elements is often deemed sufficient for many engineering purposes. This converging characteristic is also exhibited when solving case b, as evident from figures (11). Here the L2 relative deviations of the two cases are proportional to the values of H/ as expected.

Application of non-homogeneous boundary condition
Here we test the applicability of our method on solving advection-diffusion problems with non-homogeneous boundary condition. Again we set a computational domain on Ω = [−1,1] 2 where the vector field w = (2y(1−x 2 ),−2x(1−y 2 )) is set and no source term is included. Discontinuities in parts of the boundaries are introduced. At the top edge the value at the boundary is set as u ∂Ω =1 and u ∂Ω =0 everywhere else. Randomly placed 100 perforations are considered each with width of = 0.04 as shown in figure (12). In figures (13) the contours of u solved on 8×8,16×16,32×32,64×64, and 128×128 elements are given alongside the reference solution on 1024×1024 solved using standard Q1-Q1 FEM. In table 2, it is shown that the method returns grid converging results toward the reference solution as exhibited in previous tests with homogeneous boundary conditions.

Concluding remarks
In this paper, the feasibility of Crouzeix-Raviart MsFEM with bubble function enrichments for solving diffusion and advection-diffusion problems in perforated media through means of penalization methods have been demonstrated without much major constraints. The resulting method allows us to address multiscale problems with inconvenient patterns of perforations and still obtain accurate solutions between perforations. Although  Table 2: Deviation from reference solution (with non-homogeneous B.C. and = 0.04) in the given examples, the diffusion coefficient A are taken as constants, Crouzeix-Raviart MsFEM has been shown to be able to solve highly oscillatory problems (Le Bris et al., 2013a). Crouzeix-Raviart MsFEM with bubble function enrichment has shown good performance in comparison with more conventional MsFEMs especially in as far as insensitivity to size and placements of perforations is concerned. We also include the cases for non-periodic perforations where the robustness of our method is tested in more realistic circumstances.