Coursehero >> Texas >> Texas A&M >> ACCT 328

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

of Journal Computational Physics 220 (2006) 155 174 www.elsevier.com/locate/jcp Accurate multiscale nite element methods for two-phase ow simulations Y. Efendiev a b a,* , V. Ginting b, T. Hou c, R. Ewing d d Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, United States Institute for Scienti c Computation and Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, United States c Applied Mathematics, Caltech, Pasadena, CA 91125, United States Institute for Scienti c Computation and Department of Mathematics, Texas A&M University, College Station, TX 77843-3368, United States Received 10 February 2005; received in revised form 28 January 2006; accepted 5 May 2006 Available online 10 July 2006 Abstract In this paper we propose a modi ed multiscale nite element method for two-phase ow simulations in heterogeneous porous media. The main idea of the method is to use the global ne-scale solution at initial time to determine the boundary conditions of the basis functions. This method provides a signi cant improvement in two-phase ow simulations in porous media where the long-range e ects are important. This is typical for some recent benchmark tests, such as the SPE comparative solution project [M. Christie, M. Blunt, Tenth spe comparative solution project: a comparison of upscaling techniques, SPE Reser. Eval. Eng. 4 (2001) 308 317], where porous media have a channelized structure. The use of global information allows us to capture the long-range e ects more accurately compared to the multiscale nite element methods that use only local information to construct the basis functions. We present some analysis of the proposed method to illustrate that the method can indeed capture the long-range e ect in channelized media. 2006 Elsevier Inc. All rights reserved. Keywords: Multiscale; Finite element; Finite volume; Global; Two-phase; Upscaling 1. Introduction Subsurface ows, as occur in the production of hydrocarbons as well as in environmental remediation projects, are a ected by heterogeneities in a wide range of length scales. It is, therefore, very di cult to resolve numerically all of the scales that impact transport through such systems. Typically, upscaled or multiscale models are employed for such systems. The main idea of upscaling techniques is to form coarse-scale equations with a prescribed analytical form that may di er from the underlying ne-scale equations. In multiscale * Corresponding author. Tel.: +1 979 845 1972. E-mail address: efendiev@math.tamu.edu (Y. Efendiev). 0021-9991/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.jcp.2006.05.015 156 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 methods, the ne-scale information is carried throughout the simulation and the coarse-scale equations are generally not expressed analytically, but rather formed and solved numerically. Our purpose in this paper is to propose a modi ed multiscale nite element method (MsFEM) for the computations of two-phase ows. MsFEM is rst introduced in [18]. Its main idea is to incorporate the small-scale information into nite element basis functions and capture their e ect on the large scale via nite element computations. Recently, a number of multiscale numerical methods, such as residual free bubbles [6,23], variational multiscale method [18], multiscale nite element method (MsFEM) [17], two-scale nite element methods [22], two-scale conservative subgrid approaches [2,21], and heterogeneous multiscale method (HMM) [14] have been proposed. We remark that special basis functions in nite element methods have been used earlier in [4]. The generalized nite element method has also been introduced in [3] using special basis function. Multiscale nite element methodology has been modi ed and successfully applied to two-phase ow simulations in [19,20] and later in [8,1]. Arbogast [2] used variational multiscale strategy and constructed a multiscale method for two-phase ow simulations. In this paper, we propose a multiscale nite element approach in which the basis functions are constructed using the solution of the global ne-scale problem at initial time (only). The heterogeneities of the porous media are typically well represented in the global ne-scale solutions. In particular, the connectivity of the media is properly embedded into the global ne-scale solution. Thus, for the porous media with channelized features (where the high/low permeability region has long-range connectivity), this type of approach is expected to work better. Indeed, our computations show that our modi ed approach performs better, for porous media with channelized structure, than the approaches in which the basis functions are constructed using only local information. We present some analysis to justify our numerical observations. For the analysis, we use a pressure-streamline coordinate system at initial time for a simpli ed channelized media. In this coordinate system, one can perform asymptotic expansion and show that the variations of leading order pressure across streamlines are negligible, and the pressure depends smoothly on the pressure at the initial time. Furthermore, we show that global basis functions can represent the leading order pressure accurately. In our numerical simulations, we have used cross-sections of recent benchmark permeability elds, such as the SPE comparative solution project [10], in which the porous media have a channelized structure and a large aspect ratio. We would like to remark that our proposed approach is di erent from the oversampling method for multiscale nite element methods [18]. In particular, we use only the global solution of the ne-scale problem at an initial time to extract boundary conditions for the basis functions. On the other hand, the oversampling technique uses the solutions of the larger problems to construct the basis functions directly. Moreover, the proposed multiscale nite element solution is accurate at the initial time. Finally, we would like to note that the global solutions in upscaling procedures have been previously used in [7], which motivated our work. The authors in [7] show that the upscaled models that use the local information tend to perform worse for channelized porous media. Global information within mixed multiscale nite element methods was rst used in [1]. In this paper, we perform numerical tests where the global ow direction has changed. We have also tested various ranges of mobility and observed very good agreement when modi ed basis functions are used. Finally, we have used modi ed basis functions for certain linear and nonlinear parabolic equations. We have observed an order of magnitude of improvement in the error for permeability elds from the SPE comparative solution project. This paper is organized in the following way. In the next section we present the details of modi ed multiscale nite element methods. The numerical results are presented in section three. In Appendix A, we present some theoretical results related to the modi ed multiscale nite element method. 2. Modi ed multiscale nite element methods We consider two-phase ow in a reservoir X under the assumption that the displacement is dominated by viscous e ects; i.e. we neglect the e ects of gravity, compressibility, and capillary pressure. Porosity will be considered to be constant. The two phases will be referred to as water and oil, designated by subscripts w and o, respectively. We write Darcy s law, with all quantities dimensionless, for each phase as follows: vj k rj S k rp; lj 2:1 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 157 where vj is the phase velocity, k is the permeability tensor, krj is the relative permeability to phase j (j = o, w), S is the water saturation (volume fraction), p is pressure and lj is the viscosity of phase j (j = o, w). In this work, a single set of relative permeability curves is used and k is assumed to be a diagonal tensor. Combining Darcy s law with a statement of conservation of mass allows us to express the governing equations in terms of the socalled pressure and saturation equations: r k S k rp q; oS v rf S 0; ot 2:2 2:3 where k is the total mobility, f is the fractional ow of water, q is a source term and v is the total velocity, which are respectively given by: k rw S k ro S k rw S =lw ; f S ; lw lo k rw S =lw k ro S =lo v vw vo k S k rp: k S 2:4 2:5 The above descriptions are referred to as the ne model of the two-phase ow problem. Typical boundary conditions for (2.2) considered in this paper are xed pressure at some portions of the boundary and no- ow on the rest of the boundary. For the saturation Eq. (2.3), we impose S = 1 on the in ow boundaries. For simplicity, in further analysis we will assume q = 0. The upscaling of two-phase ow systems is discussed by many authors [9,5,13]. In most upscaling procedures, the coarse-scale pressure equation is of the same form as the ne-scale Eq. (2.2), but with an equivalent grid block permeability tensor k* replacing k. For a given coarse-scale grid block, the tensor k* is generally computed through the solution of the pressure equation over the local ne-scale region corresponding to the particular coarse block [12]. Coarse-grid k* computed in this manner has been shown to provide accurate solutions to the coarse-grid pressure equation. As we mentioned in Section 1, for channelized porous media, the global information can be used in calculation of e ective coarse-grid permeability [7], but these upscaling approaches are not exact at the initial time. The objective of this work is to propose an accurate multiscale method. We will use the multiscale nite element framework, though a nite volume element method is chosen as a global solver. Finite volume method is chosen because, by its construction, it satis es the numerical local conservation which is important in groundwater and reservoir simulations. Let Kh denote the collection of coarse elements/rectangles K. Consider a coarse element K, and let nK be its center. The element K is divided into four rectangles of equal area by connecting nK to the midpoints of the element s edges. We denote these quadrilaterals by Kn, where n 2 Zh(K), are the vertices of K. Also, we denote Zh = KZh(K) and Z 0 & Z h the vertices which do not lie on the Dirichlet h boundary of X. The control volume Vn is de ned as the union of the quadrilaterals Kn sharing the vertex n. The key idea of the method is the construction of basis functions on the coarse grids, such that these basis functions capture the small-scale information on each of these coarse grids. The method that we use follows its nite element counterpart presented in [18]. The basis functions are constructed from the solution of the leading order homogeneous elliptic equation on each coarse element with some speci ed boundary conditions. Thus, if we consider a coarse element K that has d vertices, the local basis functions /i, i = 1, . . . ,d are set to satisfy the following elliptic problem: r k r/i 0 /i gi on oK; in K; 2:6 for some function gi de ned on the boundary of the coarse element (or representative volume element, RVE) K. Hou et al. [18] have demonstrated that a careful choice of boundary conditions would improve the accuracy of the method. In previous ndings, the function gi for each i is chosen to vary linearly along oK or to be the solution of the local one-dimensional problems [19] or the solution of the problem in a slightly larger domain is chosen to de ne the boundary conditions. The boundary conditions for the basis functions that are used in this paper will be discussed later. We will require /i(xj) = dij. Finally, a nodal basis function associated with the vertex xi in the domain X is constructed from the combination of the local basis functions that share this xi 158 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 and zero elsewhere. We would like to note that one can use an approximate solution of (2.6) when it is possible. For example, in the case of periodic or scale separation cases, the basis functions can be approximated using homogenization expansion (see [15]). This type of simpli cation is not applicable for problems considered in this paper. Next, we denote by Vh the space of our approximate pressure solution, which is spanned by the basis functions f/j gxj 2Z 0 . Then we formulate the nite dimensional problem corresponding to nite volume element h formulation of (2.2). A statement of mass conservation on a coarse control volume Vx is formed from (2.2), where now the approximate solution is written as a linear combination of the basis functions. Assembly of this conservation statement for all control volumes would give the corresponding linear system of equations that can be solved accordingly. The resulting linear system has incorporated the ne-scale information through the involvement of the nodal basis functions on the approximate solution. To be speci c, the problem P now is to seek ph 2 Vh with ph xj 2Z 0 pj /j such that h Z h k S k rp n dl 0; 2:7 oV n for every control volume Vn X. Here ~ de nes the normal vector on the boundary of the control volume, oVn n and S is the ne-scale saturation eld at this point. We note that concerning the basis functions, a vertex-centered nite volume di erence is used to solve (2.6), and using the harmonic average to approximate the permeability k at the edges of ne control volumes. The main idea of the modi ed multiscale nite volume element method (MsFVEM) is to use the solution of the ne-scale problem at time zero to determine the boundary conditions for the basis functions. The basis functions are constructed using these boundary conditions. To describe the method, we denote the solution of (2.2) at time zero by pinit(x). For simplicity, we will assume S = 0 at time zero. In de ning pinit(x), we use the actual boundary conditions of the global problem. pinit(x) depends on global boundary conditions, and, generally, is updated each time when global boundary conditions change. For some special cases, one does not necessarily need to update pinit when boundary conditions change. We will discuss it later. The boundary conditions in (2.6) for modi ed basis functions are de ned in the following way. For each rectangular element K with vertices xi (i = 1, 2, 3, 4) denote by /i(x) a restriction of the nodal basis on K, such that /i(xj) = dij. At the edges where /i(x) = 0 at both vertices, we take boundary condition for /i(x) to be zero. Consequently, the basis functions are localized. We only need to determine the boundary condition at two edges which have the common vertex xi (/i(xi) = 1). Denote these two edges by [xi 1, xi] and [xi, xi+1] (see Fig. 2.1). We only need to describe the boundary condition, gi, for the basis function /i, along the edges [xi, xi+1] and [xi, xi 1]. If pinit(xi) 6 pinit(xi+1), then gi x j xi ;xi 1 pinit x pinit xi 1 ; pinit xi pinit xi 1 1 2pinit xi gi x j xi ;xi 1 pinit x pinit xi 1 : pinit xi pinit xi 1 If pinit(xi) = pinit(xi+1) 6 0 then gi j xi ;xi 1 /0 x i pinit x pinit xi 1 ; where /0 x is a linear function on [xi, xi+1] such that /0 xi 1 and /0 xi 1 0. Similarly, i i i gi 1 j xi ;xi 1 /0 x i 1 1 2pinit x i 1 pinit x pinit xi 1 ; 2:8 where /0 x is a linear function on [xi, xi+1] such that /0 xi 1 1 and /0 xi 0. If pinit(xi) = pinit(xi+1) 0, 6 i 1 i 1 i 1 init init then one can also use simply linear boundary conditions. If p (xi) = p (xi+1) = 0 then linear boundary conditions are used. In the applications considered in this paper, the initial pressure is always positive. Finally, the basis function /i is constructed by solving (2.6). The choice of the boundary conditions for the basis functions is motivated by the analysis. In particular, we would like to recover the exact ne-scale solution along each edge if the nodal values of the pressure are equal to the values of exact ne-scale pressure. This is the underlying idea for the choice of boundary conditions. Using this property and Cea s lemma one can show that the pressure obtained from the numerical solution is equal to the underlying ne-scale pressure. Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 159 i x )=1 i i x i+1 )=0 x i+1 xi i x x i x )=0 x i Fig. 2.1. Schematic description of nodal points. First, we would like to note that these basis functions are local. We only use the global solution, pinit, to construct the boundary conditions of the local multiscale bases. The local multiscale bases cannot be constructed directly from the global solution, pinit. It can be easily shown that these basis functions are linearly independent, and thus form a basis. Moreover, the sum of these basis functions is equal to 1 in each coarse init element, except in the elements where pP (xi) = pinit(xi+1) 6 0. Indeed, it can be directly veri ed that P4 4 the i 1 /i x 1 on P boundary. Because i 1 /i x satis es the linear elliptic equation within an element 4 K, it follows that i 1 /i x 1 in each coarse element K. One can easily modify one of the basis functions in elements with pinit(xi) = pinit(xi+1) 6 0 to guarantee their sum is equal to one. For example, changing gi+1 (see (2.8)) to gi 1 j xi ;xi 1 /0 x 2pinit1 i 1 pinit x pinit xi 1 will guarantee that their sum is equal to i 1 x one. However, as we will show next, with this modi cation the multiscale nite element solution is not exact at time zero, which is important in our applications. Next, we show that if these basis functions are used for linear elliptic equations (with k(S) = 1), then the resulting multiscale nite element solution is exact. We will show this for the multiscale nite element method. In the multiscale nite element method, the coarse-scale formulation (2.7) is given by XZ pi k S kr/i r/j dx 0; i X where pi are nodal pressure values on the coarse-grid. From the stability of multiscale nite element methods (see [16]), we have kp ph k 6 C inf kp qh k; qh qi /i . Choosing the nodal values of qi equal to the value of the ne-scale solution, one can easily where qh show that qh is equal to the ne-scale solution on the boundary of coarse blocks. This can be veri ed by direct computation. If pinit(xi) 6 pinit(xi+1), then on [xi, xi+1] we have pinit xi gi x j xi ;xi 1 pinit xi 1 gi 1 x j xi ;xi 1 pinit xi pinit x pinit xi 1 pinit x pinit xi pinit xi 1 init pinit x : init x p init x p p xi 1 pinit xi i i 1 2:9 P If pinit(xi) = pinit(xi+1) 6 0, then on [xi, xi+1] we have 160 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 pinit xi gi x j xi ;xi 1 pinit xi 1 gi 1 x j xi ;xi 1 1 0 init init init p x p xi 1 pinit xi 1 /0 x p xi /i x init i 1 2p xi 1 pinit x pinit xi 1 pinit x : 2:10 init 2p xi 1 One can show similar equalities for other edges. Furthermore, because qh satis es the underlying ne-scale equation in any coarse block and is equal to the underlying ne-scale solution on the boundary, thus it is equal to the ne-scale solution. For nite volume methods, this statement can also be proved assuming the uniqueness of the discrete solution. We omit this proof here. Our numerical results will demonstrate this. As we mentioned earlier, if pinit(xi) = pinit(xi+1) 6 0, then the basis functions do not sum up to one. To achieve the latter, one can modify the boundary conditions. But in this case, the multiscale solution at the initial time is not exact. For our computations, it is important to recover the exact ne-scale solution at the initial time. We note that even though the sum of the basis functions may not be one, in some coarse blocks where pinit(xi) = pinit(xi+1) 6 0, the basis functions still span a function that is approximately one. With direct comP4 putations, one can show that i 1 /i will be equal to pinit(x)/pinit(xi) on the edge [xi, xi+1]. Thus, it will be equal to 1 at the vertices, though slightly di erent from one along the edge. We would like to note that in the presence of source terms on the right hand side of (2.2) in order to recover the exact ne-scale solution, the basis functions need to be modi ed by incorporating the source term into the right hand side of (2.6). In particular, the corresponding right hand side in (2.6) is an original source term divided by the value of the pressure at the node i, q/pinit(xi). We would clarify the di erence between the proposed approach and the oversampling method for multiscale nite element methods [18]. Note, we use only the global solution of the ne-scale problem at the initial time and our multiscale nite element solution. On the other hand the main idea of the oversampling method is to use the solutions of the larger problems with some a priori boundary conditions. Typically, four independent solutions are constructed with some known boundary conditions. Then using these solutions, the multiscale basis functions are constructed. In the proposed approach, we use the global solution to obtain only boundary conditions for the multiscale basis functions. For two-phase ow simulations, we will use IMPES formulation (implicit pressure and explicit saturation) for the computations. Each time the pressure equation is solved and the velocity is computed. Then the velocity is used to update the saturation. For a linear problem, our approach has redundancy because it uses the ne-scale solution. Whereas for two-phase ow simulations, the pressure equation is solved many times. With our modi ed multiscale nite element method, the pressure equation will be solved using the pre-computed multiscale basis functions at the initial time. We would like to note that the permeability eld is the only function that induces the small-scale features of the ow. This small-scale information is incorporated into the multiscale basis functions. Moreover, the saturation dynamics is governed by the permeability eld and the saturation eld is generally smooth except near sharp fronts with locations determined from the heterogeneous permeability eld. Next, we would like to add a comment how to achieve the low computational complexity with multiscale nite element methods when applied to two-phase ow problems. For every saturation eld, MsFVEM produces a corresponding velocity eld. The multiscale basis functions should therefore be re-computed each time the saturation pro le changes. However, it can be shown that if the saturation is smooth within the coarse block, then the basis functions that take into account the saturation variation within the coarse block are approximately the same as the basis functions that neglect the saturation variation in the coarse block. The error made with this approximation is of order coarse mesh size. Because of this, typically in multiscale simulations (e.g., [20]), one updates the basis functions in time near a sharp front. We have observed that there is only a slight improvement if the basis functions are updated near sharp fronts. In the calculations below, the basis functions are not updated. We only update the basis functions if the global boundary conditions are changed. Below, we present a representative numerical result that compares the simulations when the basis functions are updated everywhere with the results when no update of the basis functions is performed. Using homogenization techniques for periodic media, one can show that the global multiscale nite element method does not contain the resonance errors, typically observed in multiscale nite element methods that use Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 161 local basis functions. This result will be presented elsewhere. However, this analysis does not reveal the capability of the modi ed multiscale nite element method in capturing long-range features of the ow. In the Appendix A, we present some analysis using the pressure-streamline framework that demonstrates that the modi ed multiscale nite element method is more e cient for porous media ows with long-range interactions than the standard multiscale nite element method. In a channelized medium, the dominant ow is within the channels. Our analysis assumes a single channel. Here, we brie y mention the main ndings presented in the Appendix A. Denote the initial stream function (see (A.1) in Appendix A) and pressure by g = w(x, t = 0) and f = p(x, t = 0) (f is also denoted by pinit previously). Then the equation for the pressure can be written as o op o op 2 jkj k S k S 0: 2:11 og og of of For simplicity, S = 0 at time zero is assumed. We consider a typical boundary condition that gives high ow within the channel, such that the high ow channel will be mapped into a large slab in (g, f) coordinate system (see Fig. A.1). If the heterogeneities within the channel in g direction is not strong (e.g. narrow channel in Cartesian coordinates), the saturation within the channel will depend on f. In this case, the leading order pressure will depend only on f, and it can be shown that p g; f; t p0 f; t high order terms; where p0(f, t) is the dominant pressure. The explanation of higher order terms is presented in the Appendix A. This asymptotic expansion shows that the time-varying pressure strongly and smoothly depends on the initial pressure (i.e. the leading order term in the asymptotic expansion is a function of initial pressure and time only). Because the global basis functions can recover the initial pressure exactly, the modi ed basis can capture the global pressure more accurately. In the Appendix A, we discuss more extensively the advantages of the modi ed multiscale basis functions. We would like to note that our goal is to construct a set of basis functions for the ow equation at initial time that can be used to solve the ow equation on the coarse grid at later times. This is a very important for multi-phase ow simulations, because the ow equations are solved many times and solving the ow equations is CPU demanding. By constructing a set of basis functions at initial time, the ow equation is solved on the coarse grid at later times. Note that for k(x) with scale-separation, this question can be answered within homogenization theory. Moreover, one can construct the homogenized coe cients, k*. Here, we show that for the elds with strong non-local e ects, one can construct basis functions and project the solution into the coarse dimensional space. We would like to note that the analysis presented in the Appendix A is for a single-channel ow, and can be extended to some more complicated ow scenarios. This is only a simpli ed model that allows us to demonstrate the importance of global information in constructing the basis functions. We would like to note that global basis functions capture small-scale information similar to the standard multiscale basis functions, i.e. in the case of scale separation, the convergence of the modi ed multiscale nite element methods is similar to that of the standard multiscale nite element method. This can be established using homogenization techniques. It is important to remark that, one does not know, in general, which channels will be active uid carriers and that the latter depends on global boundary conditions. The modi ed multiscale basis functions embed these global features as well as global boundary conditions into the basis functions. 3. Numerical results In this section, we present representative simulation results for ux functions f(S) with viscosity ratio lo/lw = 5. We have tested higher viscosity ratios and observed very similar results. In all cases the systems are considered to be one of the layers of the benchmark test, the SPE comparative project [10] (upper Ness layers). These permeability elds are highly heterogeneous, channelized, and di cult to upscale. In Fig. 3.1 we depict the log-permeability for one of the layers. Simulation results are presented for the total ow rate and the oil cut as a function of pore volume injected (PVI). Note that the oil cut is also referred to as the fractional ow of oil. The oil cut (or fractional ow) is de ned as the fraction of oil in the produced uid and is given by qo/qt, where qt = qo + qw, R with qo and qw being the ow rates of oil and water at the production edge R of the model. In particular, qw oXout f S v n dx, qt oXout v n dx, and qo = qt qw, where oXout is the 162 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 8 6 4 2 0 Fig. 3.1. Log-permeability for one of the layers of upper Ness. outer ow boundary. We will use the notation Q for total ow qt and F for fractional ow qo/qt in numerical Rt results. Pore volume injected, de ned as PVI V1p 0 qt s ds, with Vp being the total pore volume of the system, provides the dimensionless time for the displacement. When using multiscale nite element methods for two-phase ow, one can update the basis functions near the sharp fronts. Indeed, sharp fronts modify the local heterogeneities and this can be taken into account by re-solving the local equations, (2.6), for basis functions. If the saturation is smooth in the coarse block, it can be approximated by its average in (2.6), and consequently, the basis functions are not needed to be updated. It can be shown that this approximation yields rst-order errors (in terms of coarse mesh size). In our simulations, we have found only a slight improvement if the basis functions are updated, thus the numerical results for the modi ed MsFVEM presented in this paper do not include the basis function update near the sharp fronts. In all numerical examples, related to the SPE comparative solution project, the ne-scale eld is 220 60, while the coarse-scale eld is 22 6. We have observed similar results for other coarse grids. We consider two types of boundary conditions. For the rst type of boundary conditions, we specify p = 1, S = 1 along the x = 0 edge and p = 0 along the x = 1 edge. On the rest of the boundaries, we assume no ow boundary condition. We call this type of the boundary condition as side-to-side. The other type of boundary conditions is obtained by specifying p = 1, S = 1 along the x = 0 edge for 0.5 6 z 6 1 and p = 0 along the x = 1 edge for 0 6 z 6 0.5. On the rest of the boundaries, we assume no ow boundary condition. We will be also considering changing boundary conditions, where the boundary conditions are changed from one type to another at certain time. The objective of our rst set of results is demonstrate that the proposed procedure is exact for single phase simulations. In Fig. 3.2, the crossplot between the total ow rate (qt) for ne-scale solutions and the corresponding multiscale solutions for 50 layers of the upper Ness is plotted. In the left gure the crossplot is depicted for modi ed MsFVEM and in the right plot it is depicted for the standard MsFVEM. Every point in this gure corresponds to one of the layers of the upper Ness (total 50 layers) of the SPE comparative solution project [10]. The results corresponding to the modi ed MsFVEM are exact, while there is a deviation in the results of the standard MsFVEM. Next, we present some representative ow results. In Fig. 3.3, the fractional ow (F = qo/qt, left gure) and the total ow (Q = qt, right gure) curves are plotted for the layer number 43. One can see clearly that the modi ed MsFVEM method gives nearly exact results for these integrated responses. The standard MsFVEM to tends overpredict the total ow rate at time zero. This initial error persists at later times, and gives about 15% error at later times for both the total production and fractional ow rates. This phenomena is often observed in upscaling of two-phase ows. In Fig. 3.4, the saturation maps are plotted at PVI = 0.5. This result is representative for saturation pro les at earlier and later times. We can see from these gures that standard MsFVEM (middle gure) tends to miss some of the ne features of the ow. For example, at the lower left corner we observe an overestimation of a small saturation pocket, which does not exist in the ne-scale saturation map (left gure). Moreover, the closer look at the result obtained using the standard MsFVEM shows Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 Modified MsFVEM 500 450 163 Standard MsFVEM 600 Error = 16 % 500 400 350 300 Q m s fv e m 250 200 400 s fv e m Qm 100 200 Q fine 300 200 150 100 50 0 0 100 300 400 500 0 0 100 200 Q 300 400 500 fine Fig. 3.2. Total ow comparison for 50 layers of upper Ness. 1 0.9 0.8 0.7 250 fine modified MsFVEM standard MsFVEM 350 fine modified MsFVEM standard MsFVEM 300 0.6 0.5 0.4 150 Q 200 F 0.3 0.2 0.1 50 100 0 0 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.3. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for side-to-side ow. that there is an overprediction near the lower boundary of the layer. On the other hand, the result obtained using the modi ed multiscale nite element method looks exactly the same as the ne-scale solution. The relative L2 error for the modi ed MsFVEM is less than 5%, while the relative L2 error for standard MsFVEM is about 27%. In the next set of results, we repeat these calculations for the corner-to-corner ow. In Fig. 3.5, the fractional ow as well as the production curves are plotted. In Fig. 3.6, the saturation plots are depicted. These results are very similar to the ones obtained with the side-to-side boundary condition. For the next set of results, we consider another layer of the upper Ness (layer 59). In Fig. 3.7, both fractional ow (left gure) and total ow (right gure) are plotted. We observe that the modi ed MsFVEM gives almost the exact results for these quantities, while the standard MsFVEM overpredicts the total ow rate, and there are deviations in the fractional ow curve around PVI % 0.6. Note that unlike the previous case, fractional ow for standard MsFVEM is nearly exact at later times (PVI % 2). In Fig. 3.8, the saturation maps are plotted at PVI = 0.5. The left gure represents the ne-scale, the middle gure represents the results obtained using standard MsFVEM, and the right gure represents the results obtained using the modi ed 164 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 saturation plot at PVI=0.5 using standard MsFVEM 1 1 saturation plot at PVI=0.5 using modified MsFVEM 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 Fig. 3.4. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Side-to-side boundary condition is used. 1 0.9 0.8 280 fine modified MsFVEM standard MsFVEM fine modified MsFVEM standard MsFVEM 0.7 0.6 0.5 0.4 150 Q 200 F 0.3 0.2 0.1 50 100 0 0 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.5. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for corner-to-corner ow. saturation plot at PVI=0.5 using standard MsFVEM 1 1 saturation plot at PVI=0.5 using modified MsFVEM 1 0.9 0.9 0.9 0.8 0.8 0.8 0.7 0.7 0.7 0.6 0.6 0.6 0.5 0.5 0.5 0.4 0.4 0.4 0.3 0.3 0.3 0.2 0.2 0.2 0.1 0.1 0.1 0 0 0 Fig. 3.6. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Corner-to-corner boundary condition is used. Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 1 0.9 0.8 0.7 0.6 600 1000 165 fine modified MsFVEM standard MsFVEM 800 fine modified MsFVEM standard MsFVEM 0.5 0.4 0.3 0.2 0.1 0 0 Q 400 200 0 F 0.5 1 PVI 1.5 2 0.5 1 PVI 1.5 2 Fig. 3.7. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for corner-to-corner ow. saturation plot at PVI=0.5 using standard MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 saturation plot at PVI=0.5 using modified MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Fig. 3.8. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Corner-to-corner boundary condition is used. MsFVEM. We observe from this gure that the saturation map obtained using standard MsFVEM has some errors. These errors are more evident near the lower left corner. The results of the saturation map obtained using the modi ed MsFVEM is nearly the same as the ne-scale saturation eld. It is evident from these gures that the modi ed MsFVEM performs better than the standard MsFVEM. The next case considered involves the same permeability elds, but with changing boundary conditions. The ow is initially from left to right, as speci ed in the previous example. However, at a time of 0.6 PVI, the global boundary condition is changed such that the ow is driven from the lower left corner of the model to the upper right corner. This is achieved by specifying p = 1, S = 1 along the x = 0 edge for 0 6 z 6 0.5 and p = 0 along the x = 1 edge for 0.5 6 z 6 1 for t > 0.6 PVI. At the time when the boundary condition is changed, we update the global basis functions by using the ne-grid pressure at 0.6 PVI, and further calculations are performed using the updated basis functions. Simulation results for layer 43 are shown in Fig. 3.9. Note that at PVI = 0.6, one can observe a discontinuity ( kink ) in the fractional ow curve. This is caused by the change in boundary conditions. In Fig. 3.10, we plot saturation pro les after boundary conditions are changed, at 0.7 PVI. The multiscale simulations using the global basis functions again track the ne-grid solution much more closely than standard multiscale method with local basis functions. In Fig. 3.11, the boundary conditions at 0.6 PVI are changed to p = 1, S = 1 along the x = 0 edge for 0.5 6 z 6 1, p = 0 along the x = 1 edge for 0 6 z 6 0.5, and no ow on the rest of the boundary, i.e. the ow direction has changed from top corner 166 1 0.9 0.8 0.7 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 fine modified MsFVEM standard MsFVEM 350 fine modified MsFVEM standard MsFVEM 300 250 0.6 Q 0.5 0.4 200 F 150 0.3 0.2 0.1 50 100 0 0 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.9. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for changing boundary conditions. 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Fig. 3.10. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Changing boundary condition is used. to bottom corner. We found this to be one of the extreme cases, where the total ow rate drops signi cantly. In this case, the standard multiscale nite element method has large errors near 0.6 PVI, while we see from Fig. 3.11 that modi ed multiscale method performs very well. Again, we see an improvement when global basis functions are used in the simulations. In Figs. 3.12 and 3.13, the same test is performed for the layer 50. Again, we see an improvement rendered by the modi ed multiscale nite element method. In our next set of numerical results, we compare modi ed MsFVEM and standard MsFVEM where the basis functions are updated everywhere at each time step. The objective of these numerical results is to show that the update of basis functions does not give signi cant improvement. We only present fractional ow and production curves. For saturation plots, we have also observed almost no improvement when basis functions are updated. In Figs. 3.14 and 3.15, the fractional ow and total production curves for layers 67 and 68 are plotted. We observed from this gure that the improvement achieved by updating the basis functions is insigni cant. We have tested many of the layers and observed almost no improvement when the basis functions are updated everywhere. Our nal numerical results are for the permeability elds that are generated using two-point statistics. In general, it is easier to handle these types of heterogeneities using upscaling or multiscale methods. To generate this permeability eld, we have used GSLIB algorithm [11]. The permeability is log-normally distributed with Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 1 0.9 0.8 0.7 150 167 fine modified MsFVEM 200 fine modified MsFVEM 0.6 0.5 0.4 0.3 0.2 0.1 0 0 0 0 50 F Q 100 0.5 1 PVI 1.5 2 0.5 1 PVI 1.5 2 Fig. 3.11. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for changing boundary conditions. 1 0.9 0.8 450 fine modified MsFVEM standard MsFVEM fine modified MsFVEM standard MsFVEM 550 500 0.7 400 0.6 350 0.5 0.4 F Q 300 250 200 150 100 0.3 0.2 0.1 0 0 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.12. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for changing boundary conditions. saturation plot at PVI=0.5 using standard MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 saturation plot at PVI=0.5 using modified MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Fig. 3.13. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Changing boundary condition is used. 168 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 fine modified MsFVM standard MsFVM with update 1200 1100 1000 900 800 700 600 500 400 300 200 Q fine modified MsFVM standard MsFVM with update F 0 0 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.14. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM when the basis functions are updated everywhere and modi ed MsFVEM for layer 67. 1 0.9 0.8 45 0.7 40 0.6 35 0.5 0.4 0.3 0.2 0.1 0 0 Q F 30 25 20 15 10 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 fine modified MsFVM standard MsFVM with update 55 50 fine modified MsFVM standard MsFVM with update Fig. 3.15. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM when the basis functions are updated everywhere and modi ed MsFVEM for layer 68. prescribed variance r2 = 1.5 (r2 here refers to the variance of log k) and some correlation structure. The correlation structure is speci ed in terms of dimensionless correlation lengths in the x and z-directions, lx = 0.4 and lz = 0.04, nondimensionalized by the system length. Spherical variogram is used [11]. In this numerical example, the ne-scale eld is 120 120, while the coarse-scale eld is 12 12. In Fig. 3.16, we plot the fractional ow (left gure) as well as the total production (right gure). One can see the improvement obtained using the modi ed MsFVEM, though standard MsFVEM also performs very well. For saturation plots, we observed smaller L2 relative errors. In particular, the standard MsFVEM gives nearly 7% errors, while the modi ed MsFVEM gives less than 3% errors. Our analysis presented in the Appendix A explains why the standard multiscale nite element method work better for permeability elds generated using two-point geostatistics, where layering is parallel to Cartesian grid (see Fig. 3.17). Finally, we would like to note that we have applied the modi ed multiscale nite element methods to certain linear and nonlinear parabolic equations, where channelized permeability elds are used. We have Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0 F fine modified MsFVEM standard MsFVEM 169 fine modified MsFVEM standard MsFVEM 11 10 9 8 7 6 5 4 3 Q 0.5 1 PVI 1.5 2 0 0.5 1 PVI 1.5 2 Fig. 3.16. Fractional ow (left gure) and total production (right gure) comparison for standard MsFVEM and modi ed MsFVEM for corner-to-corner ow. saturation plot at PVI=0.5 using standard MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 saturation plot at PVI=0.5 using modified MsFVEM 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 Fig. 3.17. Saturation maps at PVI = 0.5 for ne-scale solution (left gure), standard MsFVEM (middle gure), and modi ed MsFVEM (right gure). Corner-to-corner boundary condition is used. observed an order of magnitude improvement when using the modi ed multiscale nite element methods. For example, for linear parabolic equations, the relative L2 error for standard multiscale nite element is about 3%, while for the modi ed multiscale nite element method we have observed a 0.4% error. Similarly, for nonlinear ows, such as Richards equation, we have observed an order of magnitude improvement in the relative errors. We would like also to note that the modi ed MsFVEM can be extended to 3-D. 4. Concluding remarks In this paper we propose a modi ed multiscale nite element method for two-phase ow simulations using the global ne-scale solution of single-phase equations. The main goal of this paper is to better capture the long-range features that occur in two-phase ow simulations. For this purpose, we choose permeability elds from the SPE comparative solution project [10], which have the channelized structure. We demonstrate numerically that the proposed method is capable of capturing the long-range ow features accurately for these elds. On the other hand, for more regular elds generated using two-point statistics, the standard MsFVEM works as well as the modi ed MsFVEM, though there is some slight improvement. We present some analysis that allows us to explain why the modi ed multiscale nite element method captures the global information better. We would like to note that the modi ed basis functions depend on global boundary conditions and need to be modi ed if the global boundary conditions are changed. Consequently, they are more applicable 170 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 for problems where global boundary conditions do not change frequently. Finally, we have observed an order of magnitude improvement if the global basis functions are used for some linear and nonlinear parabolic equations in heterogeneous media, where the heterogeneities have channelized structure. Acknowledgments Efendiev acknowledges Louis Durlofsky, Patrick Jenny, and Xiao-Hui Wu for many useful discussions and suggestions. The research of Y.E. is partially supported by NSF Grants DMS-0327713 and EIA-0218229 and DOE Grant DE-FG02-05ER25669. The research of T.Y.H. is partially supported by the NSF ITR Grant ACI-0204932. Appendix A. Capturing non-local e ects with global basis functions In this appendix, we present some analysis for the modi ed multiscale nite element methods for channelized porous media. For the analysis, we will use streamline-pressure coordinates. We will show that the modi ed multiscale nite element method captures the non-local e ects induced by high ow channels. To show this, we present some asymptotic results. These results basically show that time-varying pressure is strongly in uenced by the initial pressure eld. Then we show that the modi ed multiscale methods can capture non-local e ects more e ciently, because long-range information along these channels is accurately incorporated into the modi ed basis functions. We will restrict our analysis to a two-dimensional case and assume the heterogeneous porous media are isotropic, k(x) = k(x)I. The stream function is de ned as $ w = v = (v1, v2). Note that the stream function w is a scalar eld in 2-D de ned by ow=ox1 v2 ; ow=ox2 v1 : A:1 Using incompressibility, one can easily show that 1 rw 0: r k S k We will assume that the boundary conditions are prescribed by no ow on the lateral sides and p = 1 on the left vertical edge and p = 0 on the right vertical edge, and initially S = 0 inside the domain and S = 1 at the left vertical edge. The boundary conditions on the stream function depend on the pressure and velocity eld. For our speci ed boundary conditions, the stream function will be constant on lateral edges. Because stream function is de ned up to a constant, we can de ne it by zero at the bottom edge. The value at the top edge is determined by the total ow rate. It can be easily shown that op op op op rw rp k k 0: ox2 ox1 ox2 ox1 Consequently, (w,p) de ne an orthogonal curvilinear system of coordinates. Because of the orthogonality of the coordinate system, the associated Euclidean metric tensor, g is diagonal and g11 jrpj ; g22 jrwj : A:2 It can be easily shown that in the streamline-pressure coordinate system the elliptic equation becomes: ! ! p o jrwj2 oz o2 z r krz det g 2; ow jrpj2 ow op !! A:3 2 p o2 z 1 o jrpj oz rz det g : r k ow2 op jrwj2 op Furthermore, it can be easily veri ed that k jrwj : jrpj A:4 2 2 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 171 To understand how the multiscale nite element method captures the non-local e ects, we will assume that there is a high permeability channel in the porous media. Moreover, we assume that the ow along the channel is a dominant ow (i.e. we neglect the cross ow due to the mobility). This assumption holds if the permeability k(S)k is much larger in the channel compared to the permeability outside the channel. We would like to note that although this assumption can be considered a good approximation for the problems with adverse mobility ratio and high ow channels of permeability, here we use it to show that our proposed approach can capture the non-local e ects induced by the high permeability channels e ciently. Denote the initial stream function and pressure by g = w(x, t = 0) and f = p(x, t = 0). For simplicity, we will assume S = 0 at time zero. Then the equation for pressure and stream function can be written down in this curvilinear orthogonal coordinate system using (A.2), (A.4), and (A.3) o op o op jkj2 k S k S 0; og og of of ! A:5 o 1 ow o 1 ow 0: og k S og of k S jkj2 of Applying the same change of variables to saturation equation, we get oS of S of S v rg v rf 0: ot og of 2 2 A:6 2 2 0 Furthermore, it can be calculated that v rg k S kjrgj op jv0 j vg and v rf k S kjrfj op jvkj vf , k og of where jv0j = j$wj is the absolute value of the initial velocity. As for the boundary conditions, we have p(f = 0, t) = 0, p(f = 1, t) = 1, and zero Neumann boundary conditions on g = 0 and g = g0, where g0 represent the total ux at the initial time. Next, we consider mapping of the permeability eld into (g, f) coordinate system (see Fig. A.1). We assume that the ow within the channel is su ciently high such that high ow channel is mapped into a large slab (see Fig. A.1) in (g, f) coordinate system. We denote by 1 d(t) the fraction of the total ow within the channel, and assume that d(t) remains small (for further calculations we simply denote it by d). In the (g, f) coordinate system, the ow within the channel will occupy 1 d portion of the whole domain (see Fig. A.1). Note that small d implies some relation between the value of the permeability within the channel and the thickness of the channel. This relation can be easily obtained for simple channels. For our calculations, we only assume that the image of the channel in (g, f) coordinate system occupies a large slab as discussed above. Furthermore, we neglect the heterogeneities within the channel (e.g. channel s thickness in x y coordinates is small). Consequently, the saturation within the channel can be assumed to depend only on f at any time. In this case, the coe cients of the elliptic equation for pressure can be written as y ( ) x ( ) Fig. A.1. Schematic description of high ow channel and its representation in streamline-pressure coordinate system. 172 2 2 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 jkj k S jk 0 j k0 f; t 1Q1 d jk 1 j k1 g; f; t 1Qd ; k S k0 f; t 1Q1 d k1 g; f; t 1Qd ; where Q1 d denotes the region representing the channel and k0 the permeability within this channel, and Qd is the outside region and k1 the permeability outside this channel. Moreover, k0 can be assumed to vary only along the streamline, i.e. k0 = k0(f) and is much larger than k1. In this linear setting (when the pressure can be treated separately) one can perform formal expansion for pressure and show that p g; f; t p0 f; t dp1 g; f; t ; A:7 where p0 is the solution of o op0 k0 f; t 0: A:8 of of 2 o o o Indeed, it can be shown that og jkj k S o p p0 of k S o p p0 of k1 g; f; t k0 f; t 1Qd op0 . of og of From here using standard estimates, one can show that p p0 is small provided d is small. Using (A.7) one can show S g; f; t S 0 f; t dS 1 g; f; t One can also show that w g; f; t w0 g; t dw1 g; f; t Rigorous justi cations of these asymptotic expansions for coupled pressure and saturation equation is beyond the scope of our main goal and is currently under investigation. Note that in the asymptotic expansion presented above, we have considered linear pressure equation assuming that saturation mainly depends on f. The expansions (A.9) and (A.7) can be explained physically. Because the channel has high permeability, the dominant ow will be within the channel. The saturation will change in the channel rapidly (at much faster time scales compared to the saturation change outside the channel), and thus will control the pressure change in the channel. As a result, pressure will mainly vary along f. Note that p0 and S0 will vary on much faster time scales compared to other terms of the expansions. Next, we show that the basis functions of the modi ed multiscale nite element method can capture p0(f, t) e ciently. The asymptotic result shows that the dynamic pressure depends strongly and smoothly on the initial pressure eld. Because the initial pressure eld can be accurately approximated by the modi ed basis functions, we can show that the modi ed multiscale basis functions can accurately approximate the pressure eld at later times. For this purpose, we need to show that the span of the modi ed basis functions contains an appropriate one dimensional basis. Without a loss of generality, we consider a coarse element that is on the fast ow trajectory (channel) (see Fig. A.2). Denote by /i(x) the modi ed basis functions. The pressures A:9 2 1 2 Fig. A.2. Schematic description of a streamline with a coarse block. Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 173 at the points T1 and T2 are di erent because iso-pressure lines are perpendicular to the streamline. Denote the restriction of the pressure by pinit x . Consider loc b1 x pinit x pinit T 2 loc loc ; pinit T 1 pinit T 2 loc loc b2 x pinit x pinit T 1 loc loc : pinit T 2 pinit T 1 loc loc Because 1 and pinit x is in span(/1, . . . , /4), b1(x) and b2(x) are also in the span of /i, i = 1, . . . , 4. Clearly, loc b1(x) and b2(x) are linear with respect to f, and bi(Tj) = dij, because pinit is f. Moreover, b1 and b2 are linearly loc independent and b1 + b2 = 1. Consequently, the linear approximation of p0(f, t), which is the solution of (A.8), in the span of b1 and b2. There is another way to show that the span of global basis in a coarse element captures p0(f, t). Because the sum of the multiscale basis functions is one, the basis functions span 1. Moreover, the basis functions also span f because their linear combinations using the ne-scale nodal pressure values gives f. Thus, 1 and f are in the span of the global basis functions. Because p0(f, t) is the solution of (A.8), it can be approximated with linear functions with respect f. Here, we assume S0 is a smooth function, and consequently linear approximation of p0(f, t) gives rst-order accuracy (cf. (A.8)). In general, S0 can have sharp fronts in some coarse blocks, and linear approximations will not be accurate in these coarse blocks. As we mentioned earlier, the basis functions can be updated near the front. The update of basis functions will allow us to achieve better accuracy in approximating p0(f, t). As it was mentioned earlier, we have found only a very slight improvement when the basis functions are updated near the front. Here, we separate the issue of basis update near sharp interfaces from the issue of capturing of global e ects with multiscale basis functions, and currently, we are investigating how the update of the basis functions may a ect the accuracy of the multiscale nite element method. We see from the numerical results (e.g. Fig. 3.4) that when using the standard multiscale nite element methods, the saturation pro les are noisy and do not follow the streamlines very well. The explanation for this is that the span of the standard multiscale nite element basis functions does not contain the functions that only depend on f. Consequently, they can not accurately represent p0(f, t). For example, we see from Fig. 3.4 that the standard multiscale nite element method introduces an arti cial channel at the bottom left corner. This is because the locality of the basis functions misses the global connectivity of the media. Moreover, we observe that the boundaries of the main channel is not accurately represented by the standard multiscale basis functions, because they cannot simply span f. For the permeability elds generated by using two-point geostatistics with long correlation length in horizontal direction, the high permeability channels are parallel to the Cartesian coarse-grid. Since the multiscale nite element bases are linear functions along the edges of coarse-grid blocks, it can be easily shown that the span of the multiscale nite element basis functions contains the appropriate linear functions along high ow channels. Indeed, assume T2 in Fig. A.2 is on the opposite edge and the segment (high ow channel) [T1, T2] is parallel to horizontal axis. Then, we consider the span of the multiscale basis functions with 1 at two vertices of the edge containing T1 and 0 at the other two vertices. The basis functions in this span accurately represent the ow along this horizontal channel. Consequently, the standard multiscale basis functions are capable of capturing channels parallel to Cartesian grids. This is one of the reasons why the multiscale nite element method works better if high permeability channels are horizontal (or vertical). We would like to note that the above analysis is for a single-channel ow and can be extended to some more complicated ow scenarios. It is important to remark that, one does not know, in general, which channels will be active uid carriers and the latter depends on global boundary conditions. The modi ed multiscale basis functions embed these global features as well as global boundary conditions into the basis functions. References [1] J. Aarnes, On the use of a mixed multiscale nite element method for greater exibility and increased speed or improved accuracy in reservoir simulation, SIAM MMS 2 (2004) 421 439. [2] T. Arbogast, Implementation of a locally conservative numerical subgrid upscaling scheme for two-phase Darcy ow, Comput. Geosci. 6 (2002) 453 481, Locally conservative numerical methods for ow in porous media. [3] I. Babuska, G. Caloz, E. Osborn, Special nite element methods for a class of second order elliptic problems with rough coe cients, SIAM J. Numer. Anal. 31 (1994) 945 981. [4] I. Babuska, E. Osborn, Generalized nite element methods: Their performance and their relation to mixed methods, SIAM J. Numer. Anal. 20 (1983) 510 536. 174 Y. Efendiev et al. / Journal of Computational Physics 220 (2006) 155 174 [5] J.W. Barker, S. Thibeau, A critical review of the use of pseudo-relative permeabilities for upscaling, SPE Res. Eng. 12 (1997) 138 143. [6] F. Brezzi, Interacting With The subgrid World, In Numerical Analysis 1999, (Dundee) Chapman & Hall/CRC, Boca Raton, FL, 2000, pp. 69 82. [7] Y. Chen, L.J. Durlofsky, M. Gerritsen, X.H. Wen, A coupled local-global upscaling approach for simulating ow in highly heterogeneous formations, Adv. Water Resour. 26 (2003) 1041 1060. [8] Z. Chen, T.Y. Hou, A mixed multiscale nite element method for elliptic problems with oscillating coe cients, Math. Comp. 72 (2002) 541 576 (electronic). [9] M. Christie, Upscaling for reservoir simulation, J. Pet. Tech. (1996) 1004 1010. [10] M. Christie, M. Blunt, Tenth spe comparative solution project: a comparison of upscaling techniques, SPE Reser. Eval. Eng. 4 (2001) 308 317. [11] C.V. Deutsch, A.G. Journel, GSLIB: Geostatistical Software Library and User s Guide, second ed., Oxford University Press, New York, 1998. [12] L.J. Durlofsky, Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media, Water Resour. Res. 27 (1991) 699 708. [13] L.J. Durlofsky, Coarse scale models of two phase ow in heterogeneous reservoirs: Volume averaged equations and their relationship to the existing upscaling techniques, Comput. Geosci. 2 (1998) 73 92. [14] W. E, B. Engquist, The heterogeneous multi-scale methods, Commun. Math. Sci. 1 (1) (2003) 87 133. [15] Y. Efendiev, T. Hou, V. Ginting, Multiscale nite element methods for nonlinear problems and their applications, Commun. Math. Sci. 2 (2004) 553 589. [16] Y.R. Efendiev, T.Y. Hou, X.H. Wu, Convergence of a nonconforming multiscale nite element method, SIAM J. Num. Anal. 37 (2000) 888 910. [17] T.Y. Hou, X.H. Wu, A multiscale nite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1997) 169 189. [18] T. Hughes, G. Feijoo, L. Mazzei, J. Quincy, The variational multiscale method a paradigm for computational mechanics, Comput. Methods Appl. Mech. Eng. 166 (1998) 3 24. [19] P. Jenny, S.H. Lee, H. Tchelepi, Multi-scale nite volume method for elliptic problems in subsurface ow simulation, J. Comput. Phys. 187 (2003) 47 67. [20] P. Jenny, S.H. Lee, H. Tchelepi, Adaptive multi-scale nite volume method for multi-phase ow and transport in porous media, Multiscale Model. Simul. 3 (2005) 30 64. [21] I.G. Kevrekidis, C.W. Gear, J.M. Hyman, P.G. Kevrekidis, O. Runborg, C. Theodoropoulos, Equation-free, coarse-grained multiscale computation: enabling microscopic simulators to perform system-level analysis, Commun. Math. Sci. 1 (4) (2003) 715 762. [22] A.-M. Matache, C. Schwab, Homogenization via p-FEM for problems with microstructure, in: Proceedings of the Fourth International Conference on Spectral and High Order Methods (ICOSAHOM 1998) (Herzliya), vol. 33, 2000, pp. 43 59. [23] G. Sangalli, Capturing small scales in elliptic problems using a residual-free bubbles nite element method, Multiscale Model. Simul. 1 (2003) 485 503 (electronic).

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: