A multiscale finite element method for Neumann problems in porous microstructures

In this paper we develop and analyze a Multiscale Finite Element Method (MsFEM) for problems in porous microstructures. By solving local problems throughout the domain we are able to construct a multiscale basis that can be computed in parallel and used on the coarse-grid. Since we are concerned with solving Neumann problems, the spaces of interest are conforming spaces as opposed to recent work for the Dirichlet problem in porous domains that utilizes a non-conforming framework. The periodic perforated homogenization of the problem is presented along with corrector and boundary correction estimates. These periodic estimates are then used to analyze the error in the method with respect to scale and coarse-grid size. An MsFEM error similar to the case of oscillatory coefficients is proven. 
A critical technical issue is the estimation of Poincare constants in perforated domains. This issue is also addressed for a few interesting examples. Finally, numerical examples are presented to confirm our error analysis. This is done in the setting of coarse-grids not intersecting and intersecting the microstructure in the setting of isolated perforations.

1. Introduction. The modeling and simulation of porous media has wide ranging applications from subsurface flows to simulation of charge and discharge of lithium-ion batteries. At the pore-scale the simulation is constrained by the complex geometry of the material microstructure. As the small pore scale features must be resolved, however, this makes solving such a problem by direct numerical simulation very costly. In this work, we propose a Multiscale Finite Element Method (MsFEM) based on the now classical works in [10,16]. In this method, we compute many local problems with linear boundary conditions to build a coarse-grid basis. This helps to encode fine-scale information into coarse-grid basis functions. The computation of the local problems can be done in parallel in an offline phase, then computation can be completed on a cheaper coarse-grid in the online phase.
Multiscale methods of this nature have been explored vigorously in recent years. The primary motivation being multiscale phenomenon arising from oscillatory or heterogenous coefficients. Many multiscale frameworks to attack such problems have arisen, examples include the Variational Multiscale methods [18,19] and the Heterogenous Multiscale Method (HMM) [1,8] and references therein, to name a few. In this work, we will develop a method in the MsFEM framework [10,16], and utilize the theory of homogenization to obtain error bounds with respect to coarse-grids and micro-scale parameters for the case of multiscale behavior arising from perforated domains and Neumann boundary conditions.
The problem of considering partial differential equations (PDEs) in domains with perforated and complicated microstructure has a long history. This is especially true in the area of effective media theory and homogenization c.f. [5,7,24,31], where an effective, non-perforated, PDE is derived and auxiliary cell problems are proposed. Several computational approaches to perforated problems have been applied and explored. Using an HMM, the authors in [14] solve for an unknown diffusion coefficient at each coarse-grid node by solving a local auxiliary cell problem. Then, computation of the coarse-grid problem is carried out on a non-perforated domain. Utilizing a non-conforming Crouzeix-Raviart approach, a perforated MsFEM is developed for the Laplace problems in perforated domains with Dirichlet data in [21]. This approach was extended to the case of Stokes flow in complicated microstructures in [30], however, with an orthogonal splitting approach close to Variational Multiscale Methods [18,19]. In the work [4], the authors develop a Local Orthogonal Decomposition (LOD) [13,23] for perforated Neumann Problems. By truncating multiscale corrections to coarse-scale basis functions, an efficient computational scheme was developed and analyzed.
In this work, we develop and analyze a MsFEM method for Neumann problems in perforated domains. This is similar to the equation considered in [4], however, here we shall utilize homogenization techniques to obtain error estimates based on characteristic scale and coarse-grid size. The advantage of the MsFEM approach is the localized support of the basis functions compared to the patch extensions of the method based on the LOD [4,13,23]. However, due to the limitations of homogenization theory to mostly periodic problems, the error estimates obtained in MsFEM are confined to periodic media. Applicability beyond this setting must be verified numerically. Due to technical considerations arising from the Poincaré constants in perforated domains c.f. [4,29], we consider here domains only created by isolated particles so that the Poincaré constant remains uniformly bounded with respect to pore and coarse-grid sizes. We will briefly discuss details of this technicality.
The paper is organized as follows. We first give a problem setup and a brief necessary background on periodic homogenization in perforated domains. With this foundation we are able to state our MsFEM algorithm for perforated Neumann problems. Due to technical considerations of Poincaré constants in perforated and complicated domains, we give a brief overview of their possible impact to error estimates in this setting. With the assumption that these constants are bounded with respect to scale and coarse-grid block size, we are able to derive the standard error expression for MsFEM. The auxiliary results needed for this estimate are given in the appendix for clarity of presentation. Finally, we test our error estimates by implementing the method on two geometries. The first geometry has the perforations contained entirely in the coarse-grid blocks, and the second where the perforations intersect the coarse-grid blocks and the nodal points reside inside the perforations. The second example is critical to demonstrate that the technical restriction of having the perforations inside the coarse-grids is only for simplicity and the method has possible wider applicability.
2. Neumann problem in porous microstructures. In this section, we briefly give the problem statement of Neumann problems in porous microstructures. Due to the multiscale nature of such microstructures and complex geometries, we use the tools of periodic homogenization to derive the effective equation and generate auxiliary cell problems that connect macro-and micro-scale information. From this homogenization framework we develop a MsFEM method where local basis functions are computed in parallel. These basis functions are then used to compute on a coarse-grid as multiscale information is contained in the basis functions. We then present error analysis that closely follows the work of [10,16], and references therein and the standard MsFEM analysis. In this context, we extend this idea to domains with perforations. Finally, we remark on the possible effects of the microstructure on the Poincaré constants in perforated domains [4,29].
We suppose we have a Lipschitz domain Ω ⊂ R n such that it is decomposed into a solid microstructure S ε and an open connected pore space Ω ε with a characteristic pore size of ε. The interior interface of such a microstructure is denoted as Γ ε and the outer global boundary we denote as ∂Ω ε \Γ ε . Given a f ∈ L 2 (Ω ε ), we consider the following Neumann problem in Ω ε . We wish to find a c ε such that Note that we could consider a similar Neumann problem with oscillatory diffusion coefficients in addition to the pore microstructure, however in terms of MsFEM such problems are well studied both numerically and theoretically cf. [9,16,17].
2.1. Periodic homogenization. As noted, solving the equations (1) is complicated by the fact that the geometry has many scales and can be very complex. A useful tool in circumventing this issue is homogenization. In particular we employ the use of two-scale asymptotic expansions to homogenize the system [31]. These techniques are closely related to the design and analysis of MsFEM. A key assumption in many of these methods is the assumption of periodic structure, however, this assumption may be relaxed in the application of MsFEM. We suppose that our medium has additional periodic structure and can be written as where Y is the reference cell, Y * is the open pore space, and Y Γ is the interface, or perforations boundary, in the cell. We briefly recall the periodic homogenization of (1), similar derivations can also be found in [5]. We expand using the two-scale asymptotic expansion as here y = x ε , the periodic variable in the unit cell. Applying the above expansion into (1), and gathering the ε −2 terms in (1a) and ε −1 terms in (1b), we obtain −∆ yy c 0 = 0 in Y * , −∇ y c 0 · n = 0 on Y Γ , hence, c 0 (x, y) = c 0 (x). Taking the next orders in ε we obtain the cell equations with c 1 being y-periodic. We denote the average over a domain by · Y * = 1 |Y | Y * · dy, and require c 1 Y * = 0 to eliminate the arbitrary constant. We may simplify the above cell problem by writing the two-scale function Here, Q(y) ∈ (H 1 # (Y * )) n , where # signifies periodic, satisfies the following cell problems with Q being y-periodic and Q Y * = 0. Here we write I to mean the n × n identity matrix or (I) ij = δ ij . Taking the final order in ε of (1), we arrive at Averaging over the y variable, the second term vanishes on the boundary of the interface. The homogenized problem can then be written as where (D * ) ij = δ ij + (∇ y Q) ij Y * and we write c T = c 0 Ω .

2.2.
A perforated multiscale FEM. In this section, we outline the method of MsFEM for perforated Neumann problems. Due to the Neumann condition, we are able to construct a conforming multiscale finite element method. This is in contrast to the case when the problem has a Dirichlet condition on the perforations and weak conditions can be effectively utilized [3,30]. Moreover, the coarse-grid blocks may intersect the perforations and so the basis functions can have holes on portions of the boundary. Our analysis however, will focus on the case where the holes are entirely contained in each coarse-grid block. This assumption is to avoid complicated technical details in the homogenization and error estimates [6,12]. We will present numerical examples in Section 3 with coarse-grids intersecting the microstructures. We will highlight the areas where this non-intersecting assumption is made and can be weakened. We begin first with some notation. Suppose we have a domain with microstructure Ω ε , not necessarily periodic, and a quasi-uniform coarse-grid T h (of characteristic grid size h) of the domain. In our error analysis, we will suppose that the coarse-grid does not intersect the perforations along the edges. We suppose ε < h, as the case where h < ε corresponds to a full direct numerical simulation and the coarse mesh will resolve the geometry. The standard finite element error of local problems is ignored in the analysis.
We let K ∈ T h be an element in the partition (may also be triangulation) without perforations and denote K = K ∩ Ω ε to be a perforated element. We denote the boundaries of the perforations Γ K and the boundary of the element ∂K. We denote the global set of vertices N h associated with T h . For x i ∈ N h , we build the corresponding multiscale basis function φ i such that for K ∈ T h we have Here we denoted the linear classical nodel shape function associated with the global node x i ∈ N h by φ i L . Note here we have dropped the K index from the multiscale basis.
Remark 2.1. If the coarse-grid does not intersect the perforations, then, ∂K\Γ K = ∂K, and this is the case we will consider for our analysis and proofs. However, the variational form for a local problem with intersections is to find for all ψ ∈ H 1 ( K), ψ = φ i L on ∂K\Γ K . In our algorithm, we suppose that the perforations leave portions of the linear Dirchlet condition intact so that the above right hand data gives a useful basis function.

We denote the space of solutions by
Note that this approximation space is conforming. We will approximate the solutions of (1) by c h ∈ V h ms . The corresponding variational form is given by where the boundary term over Γ ε vanishes due to condition Neumann condition. We take a similar approach as in [16,17]. We denote the Lagrange interpolation operator for the multiscale basis as I h ms : V → V h ms . We let I h ms (c 0 ) be the interpolant of c 0 (the solution of (7)) using the multiscale basis φ i . More specifically, the interpolant is Here we emphasize what domain the interpolant its expressed on as in the analysis it will be important to distinguish between perforated and unperforated quantities. For example, even if x i ∈ S ε , c 0 (x i ), a solution to (7), is well defined as this function is expressed on the unperforated domain Ω and is what we refer to as a coarse-grid quantity.

2.3.
Poincaré constants in perforated domains. Before we start the error analysis of our MsFEM, it is important to note that we do not track the Poincaré constants as we suppose that such a constant is uniformly bounded with respect to the microstructural parameters of scale and separation length. For our case of isolated particles of size and separation length ε, this uniform bound is known to hold c.f. [4,29]. Here we briefly discuss some of the cases that may occur. The analysis of Poincaré constants has a long history and too vast of literature to discuss here. We will follow the method and examples introduced in [29], and the references therein, for the application of weighted Poincaré inequalities. It is known that the shape of the domain and moreover, the microstructure may have an effect on the Poincaré constant. First recall that for K, with diam(K) = h, we have where C P (h, ε) is the Poincaré constant that may depend on the size of the diameter of the domain as well as the separation length of the microstructure and scale. We take the scale of the pore and the scale of the separations to both be of order ε.
A classical illustration of domain dependence of the Poincaré constant is the dumbbell shaped domain. Suppose we take K = [0, h] 2 and remove the two thin pieces Then, we let K = K\S ε . In Figure 1, we see S ε is black, and K as the white background. Using isoperimetric inequality methods to bound C P for the dumbbell domain, [25,26] yields the pessimistic bound where C is a benign constant independent of h or ε. By using methods derived for high-contrast problems and weighted Poincaré inequalities, the authors in [29] are able to obtain a more optimistic bound by a constructive method. For the dumbbell domain it is known that A particularly bad Poincaré constant occurs when the domain is particularly tortuous. Using methods from [29] adapted to the perforated case, the authors in [4] show that for a reticulated interlocking filamented structures, Figure 2, have the bound However, in the geometries we will consider in our analysis and in our numerics we largely consider isolated particles that do not create the ill effects that may be possible in the above domains. This was first demonstrated for high-contrast domains in [29] and in perforated domains in [4]. Given that the particles are of characteristic size and separation length ε in R 2 , along with technical assumptions for the constructive approach it can be shown in [4] that where C is a benign constant independent of h or ε. The above bound is the Poincaré bound we shall suppose throughout the rest of the paper.

2.4.
Error analysis of the MsFEM. We will now show the errors introduced in the method are similar to that of the classic elliptic MsFEM error estimates. Indeed, due to the structure of (1), the homogenization procedure and related analysis are very similar to the elliptic oscillatory case considered in [16,17]. Thus, the proofs here follow closely those works. In this analysis, we suppose that the perforations do not intersect the coarse-grid, i.e. all the microstructure is in the interior of the basis functions (8). This is to avoid technical considerations of periodic unfolding methods, however, we will discuss the possible extensions in this direction. We will ignore the fine-scale discretization error throughout and consider only coarse-grid error h. This is also the case when we investigate the numerics. Finally, we will also suppose the microstructure has periodic structure and, the domain and data are sufficiently smooth for the following estimates to hold. More specifically, we suppose regularity so that c ε ∈ H 1 (Ω ε ), c 0 ∈ C 2,1 (Ω), and Q ∈ (C 1 (Ȳ * )) n . We proceed as in [16], note that the basis functions may be expanded in K as This is simply the two-scale expansion similar to (3), applied to the multiscale basis. We denote the residual in this expansion by R ε and estimate its magnitude in Lemma 2.4. Lemma Again, as with the multiscale basis for x i ∈ N h , we define φ i 0 that satisfies the homogenized equation with φ i L the classical nodal shape function associated to x i , and φ i where Q is given by (6) after rescaling, and As we would like to multiply quantities that are on the periodic cell and the homogenized domain and then rescale to the perforated media, we note here that there is some ambiguity in the notation of Q(y) for y ∈ Y * and the same function rescaled to Ω ε . When there is ambiguity in the notation we shall denote in Ω ε . In the case of linear Lagrange triangular finite elements, for example, the functions φ i 0 are a basis for solving the homogenized equations (7) with approximation order h in H 1 (Ω). Note that linear Lagrange is just a choice of convenience and other MsFEM elements have been developed, for example in the case of non-conforming elements in [11], and recent work in the setting of higher order finite elements in [15]. Using the same notation as in [16], we denote the space V h 0 ⊂ H 1 0 (Ω) to be the space of first order finite elements (order h approximation) with zero global Dirichlet conditions. More specifically, we define V h 0 to be the space spanned by φ i 0 , satisfying (18). Further, we define the Lagrange interpolation operator I h 0 : V → V h 0 , to be the operator that gives the Lagrange interpolation in the V h 0 basis. Note that I h ms (c 0 ), given by (11), satisfies −∆I h ms (c 0 ) = 0 in K, −∇I h ms (c 0 ) · n = 0 on Γ K .
Hence, we may expand, as an ansatz as in (17) I h ms this is simply the two-scale expansion similar to (3), applied to the multiscale interpolant in K. We denote the residual in this expansion by R I,ε and estimate its magnitude in Lemma 2.5. Lemma Here we may expand the first term as where φ i 0 satisfy (18). The above quantity may be expressed on the unperforated domain. The next order corrector is given by Q ε ∇I h 0 (c 0 ) in K. Finally, the local boundary layer correction θ I,ε satisfies −∆θ I,ε = 0 in K, −∇θ I,ε · n = 0 on Γ K , Further we will need the so-called global boundary corrector −∇θ ε · n = 0 on Γ ε (25b) We will need the following technical lemma that we will restate and prove in the Appendix B related to the above global boundary corrector.
Lemma 2.4. Let c ε be a solution of (1), let c 0 be a solution to (7), and Q given by (6). We suppose the global boundary layer correction θ ε satisfies (25). Then, we have the following estimate In addition, we have the following estimate for Lemma 2.5. The interpolants satisfy the following corrector estimate Proof. We may compactly write the differential operator associated to (1) as L ε , and by (21), we have that L ε I h ms (c 0 ) = 0. Hence, we may write the formal two-scale ansatz expansion in K of I h ms (c 0 ) as I h ms (c 0 ) = I h 0 (c 0 ) + ε Q ε ∇I h 0 (c 0 ) − θ I,ε + · · · , we may apply the corrector estimate (26) to I h ms (c 0 ) so that we have This is proved in Appendix B in the general setting. Again, by briefly following the arguments of [16], we let c bl be a standard bilinear approximation of c 0 . We have by standard approximation theory and elliptic regularity of c 0 that and . Thus, using the above inequality and (28), summing over K, using the bound (29) we obtain (27).
We are now in a position to state our main theorem. Again, we note that it is for the most part a translation of the elliptic oscillatory case ( [16,17]) into the language of perforated homogenization and perforated multiscale finite elements. The main differences being the need to emphasize what are perforated fine grid quantities and what can be represented on the unperforated coarse-grid.
We will be assuming sufficient smoothness on the domains Ω ε and Ω and that the perforations do not intersect the global boundary ∂Ω or the edges and vertices of the discretization T h . Moreover, we suppose that the Poincaré constants in perforated domains do not interfere with the estimate and satisfy the uniform bound (16). These assumptions may be loosened by periodic unfolding ( [6]) or carefully tracking the Poincaré constants in perforated domains ( [4,29]), however, these may be considered in future works. First we proceed by noting that from the so-called Cea's Lemma we have the following error.
Theorem 2.6. Suppose that c ε is a solution to (1) and c h satisfies the variational form (10). Further, we suppose that the Poincaré constant C P satisfies a uniform bound (16). Then, we have the following error estimate, that there exists a C > 0 independent of h and ε such that Proof. From the classical Cea's Lemma and Galerkin orthogonality we have for C > 0, independent of h and ε, that Taking c I = I h ms (c 0 ) given by (11), and using Theorem 2.7, we arrive at our result.
The above error relies on the following estimate.
Theorem 2.7. Suppose that c ε is a solution to (1) and I h ms (c 0 ) given by (11). Further, we suppose that the Poincaré constant C P satisfies a uniform bound (16). Then, we have the following error estimate, that there exists a C > 0 independent of h and ε such that Proof. Using the expansions (3) and (22) for c ε and I h ms (c 0 ) respectively and the corresponding corrector estimates (26) and (27) we have We begin by estimating each of these term by term. Recall that from (23) I h 0 (c 0 ) is in V h 0 and is a finite element approximation to (7) with the basis spanned by φ i 0 . Thus, we have the estimate Here, we used that both c 0 and I h ms (c 0 ) are able to be represented as unperforated quantities. Note that assuming sufficient smoothness of the perforations, we have Q L ∞ ≤ C and ∇Q L ∞ ≤ C/ε. Using the expression c 1 = Q ε ∇c 0 , we have element-wise Taking the integral over the whole domain on the unperforated coarse-grid quantities, squaring, and summing over K, and using the approximability of V h 0 we have Finally, using the estimates (45) and (47) in the Appendix A for θ ε , θ I,ε we have Combining the estimates (33), (34), and (35) into (32) and summing over K we have Remark 2.8. In terms of the analysis from homogenization, the assumptions that perforations do not intersect the global boundary or the boundary of the coarse-grid may also be relaxed. This can be achieved through the method of periodic unfolding cf. [6] and references therein. In this methodology, assumptions on perforations intersecting boundaries may be relaxed. In addition, corrector estimates such as those derived in Appendix B may be extended to this setting by combining periodic unfolding methods to standard corrector estimate techniques [12]. However, this comes at the cost of much technical overhead and we focus on the theory without intersections and verify that it may be extended via our numerical results.

Numerical examples.
We will now demonstrate our error estimate from Theorem 2.6 on a few numerical examples. We begin by setting up how we will compute the numerical solutions by constructing multiscale basis functions. We will set up numerical experiments that vary both scales ε and coarse-grid value h. This will be done for coarse-grids not intersecting perforations and coarse-grids intersecting the perforations to show that the theory may be extended beyond the theory of non-intersection.
3.1. Problem setup. In our numerics we have the flexibility to add in an oscillatory coefficient, k ε , in addition to oscillations created by the multiscale geometry. We restate the main problem (1) here for our specific case. Recall, we wish to solve the following boundary value problem in the perforated domain Ω ε = Ω\S ε , where Ω = [0, 1] × [0, 1] ⊂ R 2 , S ε is the collection of the perforations, and Γ ε is the interface as shown in Figure 3. We wish to solve Again, we will not utilize the flexibility in our code to add in the oscillatory coefficient in our results, we just state it here to show that it is an easy extension.

Figure 3. Perforated solution domain with coarse mesh of finite elements
To begin construction of the multiscale basis functions, let T H be a partition of the perforated domain Ω into elements K, which are shown in Figure 4. Then, we denote the perforated element K = K ∩ Ω ε and we have that Ω ε = K∈T H K. For a single element, we denote with B the domain of the perforation inside K and ∂B is its boundary. With "•", as shown in Figure 4, we indicate the vertices of K in which the coarse-grid nodal values are calculated, and with N we denote the number of vertices in T H . Then, the multiscale basis functions φ M i , i = 1, 2, . . . , N , are solutions to the following local problems which we solve for each coarse finite element K ∈ T H and {φ L i (x)} N i=1 is the standard piecewise bilinear basis. We show an example of a perforated multiscale basis function with 9 holes in Figure 5. 3.2. Numerical setup and methods. Again, we do not consider the case when k ε (x) is a highly oscillating function of x because we are interested only in oscillations coming from the perforations. Therefore, we take k ε to be a constant. In the MsFEM framework, we use finite elements which are perforated squares and the coarse-grid size is denoted by h. The small parameter ε characterizes the size of the microstructures. For a consistent numerical analysis of the convergence of the method, we need to decrease uniformly h and ε. Hence, we consider solution domains Ω ε with periodically arranged identical perforations. However, the method as a computational tool is not restricted to periodic media, only the analysis is restricted to such simple cases.   For the reference solution, we use the total fine-scale grid used in the computation of the local problems with size n f . Thus, the fine-grid is also of order N micro , however, we lose the parallel structure and it must be solved in entirety. This is done using standard bilinear finite elements on the fine-grid rectangles. We run numerical simulations with different number of holes per coarse element. We will vary both h and ε and record the results in tables. In Figure 6 we show example of two geometries and the placement of the coarse-grid elements. In Figure  6(a) we show macro-elements with 4 holes entirely included in them. The second geometry has coarse elements with the edges and vertices inside the perforations, see Figure 6(b). The first example matches our theory, while the second demonstrates the applicability beyond the assumption of the perforations not intersecting the coarse-grids.
We solve the local problems (38) using the standard Finite Element Method with linear Lagrange triangular elements. For the numerical integration we use a Gaussian quadrature rule and we use the preconditioned Stabilized Biconjugate Gradient Method as a linear solver. We show numerical results for the L 2 norm and the L ∞ norm. We recall that for the standard MsFEM , the following L 2 error estimate holds, for the case of oscillatory coefficients in [16], and the proof of Theorem 2.6 in the case for perforated Neumann problems where this estimate can be achieved by a use of the standard Aubin-Nitsche lemma.
The following improved error estimate is shown numerically again in [16], 3.3. Numerical results. We now will present our numerical results. Throughout we take the coefficient to be k ε (x) = 5 and the right-hand side to be f (x) = 16. We begin by performing tests on the first geometry in Figure 6(a) where the coarse-grid does not intersect the geometry. Then, we will perform similar test on the second geometry where the coarse-grid nodes and edges intersect the perforations in Figure  6(b). We perform two tests on each geometry. Fixing h then decreasing ε and then fixing ε and decreasing h.

3.3.1.
Non-intersecting coarse-grid. We now begin with tests on Figure 6(a). Here we keep h = 1 16 fixed and decrease ε. In this test case we have 16 2 = 256 coarse-grid blocks with varying number of holes per block. We start with 1 perforation per block and we finish with 8 2 = 64 holes per macro finite element. The results from Table 1 show a convergence rate, which is even better than the theoretical one given by (39). We also observe that when ε becomes very small, which is equivalent to the ε h 1 2 term getting very small, the convergence rate starts to decline. This is most likely due to the fact that the h 2 term becomes dominant for very small values of ε. In Figure 7 and 8, we show numerical results for different values of ε. Now keeping ε = 1 128 fixed and decrease h. In this experiment we have 128 2 = 16384 perforations and a varying number of coarse-grid blocks. The convergence rates are given in Table 2. When h is relatively large, the error first decreases, and then, when h becomes small enough, the error starts to increase as predicted by the theoretical estimate with a convergence rate close to the improved one (40). Note that the final error in in Table 2, when h = 0.0078125, we are in the h = ε regime.   Here we may experience resonance errors from the order ( ε h ) 1 2 error terms. The microscale and the MsFEM solutions are shown in Figures 9 and 10. 3.3.2. Intersecting coarse-grid. Now we will perform the same tests, but this time on a geometry that has the perforations intersecting the coarse-grids, moreover, the nodes of the coarse-grid are directly centered on the perforations. Let us assume that we have some fixed coarse-grid size h and we want to vary the number of holes    Figures 11(a), 11(b), and 11(c), respectively.
We now begin with tests on 6(b). Here we fix h = 1 16 and decrease ε. In this test case we have 16 2 = 256 coarse-grid blocks with different number of holes per block.
Here h is sufficiently small and, as we can see from Table 3, we obtain a very good convergence rate, which is even better than the theoretical result (39) and coincides    decreases and then, when h becomes sufficiently small we obtain a convergence rate, which tends to the improved one (40). The microscale and the MsFEM solutions are given in Figures 14 and 15.

Conclusion.
In this work we developed and analyzed a multiscale finite element for domains with porous microstructures. The standard error estimates for the case of oscillatory coefficients were shown to also hold in the case of perforated media. Due to the complexity of the Poincaré constants in perforated domains a discussion on their dependence of the microstructure is warranted and was presented. It should be noted that for our case of isolated particles the geometry has no effect on this constant. Further, to develop and support the theory, we provided corrector estimates and boundary correction estimates. Finally, we implemented the algorithm on two different geometries and recorded the results in tables varying geometry size and coarse-grid sizes. The results were in good agreement with the theory. Future works include the incorporation of these methods into more complicated nonlinear and multiphysics problems such as lithium ion transport to speed up computation. In addition, the investigation of classical improvements methods such as oversampling would also be of interest.
Appendix A. Estimating boundary correctors θ ε , θ I,ε . We will need the following inequalities. The following Gagliardo-Nirenberg interpolation inequality [28] w H 1/2 (∂K) ≤ C w 1/2 and the trace inequality for w ∈ H 1 (K) Figure 11. Coarse-grid blocks with different number of perforations per block. we refer the reader to [2] for proof. From the global boundary correction equation, −∇θ ε · n = 0 on Γ ε (43b) we see that we may easily obtain the bounds using the interpolation inequality (41) on the global boundary ∂Ω ε \Γ ε = ∂Ω, (since the perforations do not intersect the global boundary) we have using the classical estimate from [22] θ ε H 1 (Ωε) ≤ C c 1 H 1/2 (∂Ω) ≤ C c 1 1/2  By smoothness, we have Q L ∞ ≤ C and ∇Q L ∞ ≤ C/ε , and so using a trace inequality we have the estimate Thus, since we suppose global regularity of the homogenized quantity c 0 ∈ C 2,1 (Ω), the solution to (7), and therefore is in H 3 (Ω). We have the bound on the quantity of interest εθ ε H 1 (Ωε) ≤ Cε 1/2 . Similarly on the element K, from (24) we see that using (41) and estimates of L ∞ estimates of Q, ∇Q we obtain Finally, using (42), we obtain θ I,ε H 1 ( K) ≤ h −1/2 ε −1/2 C|∇c 0 | H 1 (K) + C(h −1/2 |∇c 0 |) Here we have used the global regularity c 0 ∈ C 2,1 (Ω). Thus, we have after summing over K, εθ I,ε H 1 (Ωε) ≤ C ε h 1/2 . Remark B.2. In addition, to avoid technical details we suppose the perforations do not intersect the global boundary. These assumptions may be relaxed with the use of periodic unfolding techniques to handle perforations intersecting the boundaries [6,12].
Here, we used ω ε = 0 on ∂Ω, (6b), and (43b). On the right hand side, using (7), and for clarity moving to indicial notation we have Taking a closer look at the first term we see that using n j=1 ∂ ∂x j δ ij + ∂Q i ∂x j = 0, from (6) and using the symmetry of D * , we obtain This last term is precisely equation 2.18 of [27], where the oscillatory coefficients are constant and the oscillations arise from the domain alone. In the same way as [27],