Mass balance eqs. for multi-component multi-phase flow For multi-phase flow we require mass balance of each component: ∂ � � dt ( φm α,i ρ i s i ) + ∇ · ( m α,i ρ i v i ) = q α , i i Thus, we accumulate the mass of component α in each phase. Pressure equation The system of mass balance equations may be written on the form ∂ dt ( φA [ s i ]) + ∇ · ( A [ v i ]) = [ q α ] , i = a, l, v, α = w, o, g. Multiplying with A − 1 and summing the eqs. gives the pressure eq: ∂φ � A − 1 ∂A � � ∇· v i +[ v i ] · ( A − 1 ∇ A ) 1 = 1 A − 1 [ q α ] . dt + φ [ s i ] · 1 + dt i Applied Mathematics 21/89
Production processes Primary production: puncturing the balloon When the first well is drilled and opened for production, trapped hydrocarbons start flowing (toward the well) due to over-pressure. Secondary production: maintaining reservoir flow As pressure drops, less oil and gas is flowing. To maintain pressure and push out more profitable hydrocarbons one starts injecting water or gas (or CO 2 in liquid form) into the reservoir. Enhanced oil recovery (EOR): altering reservoir flow EOR is used to alter flow patterns and fluid properties. E.g.,: Inject foam to more efficiently push out oil. Inject polymers to change the flow properties of water. Inject solvents to, e.g., develop miscibility with gas. Applied Mathematics 22/89
Discretizing the pressure equation (for incompressible flow) Three-phase three-component flow (black-oil model) v i = − Kk ri ( ∇ p i + ρ i g ∇ z ) , µ i �� � ∂φ ∂p l � � dt + ∇ · c i v i · ∇ p l = q, + φ c i s i v i + dp l i i i where the c j s are phase compressibilities, and q = 1 A − 1 [ q α ]. Incompressible flow with constant capillary pressure v = − λ ( ∇ p l + ωg ∇ z ) , ∇ · v = q. where v = � i v i , λ i = Kk ri /µ i , λ = � i λ i , and ω = � i λ i ρ i /λ . Applied Mathematics 23/89
Mass conservative methods for the pressure equation Mass conservation on a grid G = { Ω i } A method is (mass) conservative if � � � ∇ · v dx = v · n ds = q dx. Ω i ∂ Ω i Ω i for each grid cell Ω i in Ω (the reservoir). Mass conservative methods Finite-volume methods Mixed finite-element methods Mimetic finite-difference methods Applied Mathematics 24/89
Finite-volume methods Finite volume methods have the following characteristics: Each equation imposes mass conservation for one grid cell: � � v · n ds = q dx. ∂ Ω i Ω i Pressure is cell-wise constant. The outflux � v · n ds ∂ Ω i is estimated using Darcy’s law with the pressure gradient defined in terms of neighboring pressure values. Applied Mathematics 25/89
The two-point flux-approximation (TPFA) scheme TPFA: Flux across γ ij = ∂ Ω i ∩ ∂ Ω j is expressed as a function of p i and p j , the cell pressures in Ω i and Ω j , respectively. v = t (p − p ) ij ij ι j Ω j Ω ι Applied Mathematics 26/89
The two-point flux-approximation (TPFA) scheme Discretized equations Write the pressure equation (for incompressible flow) as follows: −∇ · λ ∇ p = f, where f = q + g∂ z ( λω ). Then � � − λ ∇ p · n ds f i = f dx = Ω i ∂ Ω i � � − λ ∇ p · n ij ds = ∂ Ω i ∩ ∂ Ω j j � � = v ij ≈ t ij ( p i − p j ) . j j Here n is the unit normal on ∂ Ω i pointing into Ω j and n ij is the unit normal on ∂ Ω i ∩ ∂ Ω j pointing into Ω j . Applied Mathematics 27/89
The two-point flux-approximation (TPFA) scheme Derivation of transmissibilities On γ ij = ∂ Ω i ∩ ∂ Ω j the TPFA scheme estimates � p j − p i � ∇ p ≈ ∂p ij = n ij , dist( c i , c j ) and assumes that on each interface λ = λ ij where λ ij is a distance-weighted harmonic average of λ i,ij = n ij · λ i n ij and λ j,ij = n ij · λ j n ij . E.g., for an interface between two cells in the x –coordinate direction in a Cartesian grid with cell dimensions (∆ x i , ∆ y i , ∆ z i ): � − 1 � ∆ x i + ∆ x j λ ij = (∆ x i + ∆ x j ) . λ i,ij λ j,ij Applied Mathematics 28/89
The two-point flux-approximation (TPFA) scheme Derivation of transmissibilities, cont. Thus, � v ij = λ ∇ p · n ij ds ≈ −| γ ij | λ ij ∂p ij = t ij ( p i − p j ) , ∂ Ω i ∩ ∂ Ω j where � − 1 � ∆ x i + ∆ x j t ij = 2 | γ ij | . λ i,ij λ j,ij The TPFA scheme results in a linear system of the form � � j t ij if k = i, Ap = f where a ik = − t ik if k � = i. Applied Mathematics 29/89
The two-point flux-approximation (TPFA) scheme Main limitation: Example: Assume that K is a constant tensor and let γ ij be an interface between two adjacent grid cells Ω i and Ω j in the x –coordinate direction in a Cartesian grid. Then � � � v ij = − k xx p x + k xy p y + k xz p z ds. γ ij Since p i and p j can not be used to estimate p y and p z , we see that the TPFA scheme neglects the contribution from k xy p y and k xz p z . Applied Mathematics 30/89
The two-point flux-approximation (TPFA) scheme K -orthogonality The TPFA scheme only convergent for K -orthogonal grids, i.e., grids where each cell is a parallelepiped and n ij · Kn ik = 0 , ∀ Ω i ⊂ Ω , n ij � = ± n ik . Examples: n 2 n 1 γ K K K i ij j n T K n 2 =0 1 ∆ ∆ i j Applied Mathematics 31/89
Multi-point flux-approximation (MPFA) schemes To obtain consistent interface fluxes for grids that are not K -orthogonal, one must also estimate partial derivatives in coordinate directions parallel to the interfaces. Solid line Ω 3 = primal grid. Ω 1 Ω 2 Dashed line Ω 6 Ω 4 Ω 5 = dual grid. To compute the flux across γ 25 one uses cell pressures in cells 1–6. Applied Mathematics 32/89
Multi-point flux-approximation (MPFA) schemes To derive an MPFA scheme it is common to introduce a dual grid. Ω 3 Ω 1 Ω 2 Ω 6 Ω 4 Ω 5 Example: Control-volume finite element method Let V be the bilinear (trilinear) finite element space with respect to the dual grid, i.e., V is spanned by bilinear (trilinear) nodal basis functions where the nodes are cell centers in the primal grid. � � Find p ∈ V such that − λ ∇ p · n ds = f dx. Ω i Ω i Applied Mathematics 33/89
Multi-point flux-approximation (MPFA) schemes The O-method Ω 1 Ω 4 O-method interaction region: x 1 x 23 x 4 Region confined by edges or faces II Ω 1 I Ω 4 connecting the cell centers. x 34 x 12 IV The shaded region is the interaction Ω 2 III Ω 3 x 2 region associated with Ω 1 –Ω 4 . Ω 2 x 3 x 14 Ω 3 Function space V for the O-method (quadrilateral grid): Each v ∈ V is linear on each quadrant of IR, and continuous at junction between the dual and primal grid. The flux normal to the interfaces of the primal grid − λ ∇ v · n ij is continuous at junction between the dual and primal grid. Applied Mathematics 34/89
Mixed finite-element method Mixed formulation of pressure equation Pressure equation (incompressible flow): v = − λ ( ∇ p + ωg ∇ z ) , ∇ · v = q. Mixed formulation (no-flow boundary conditions): Find ( p, v ) ∈ L 2 (Ω) × H div 0 (Ω) such that � � � v · λ − 1 u dx − p ∇ · u dx = − ωg ∇ z · u dx, Ω Ω Ω � � l ∇ · v dx = ql dx, Ω Ω for all u ∈ H div 0 (Ω) and l ∈ L 2 (Ω). Applied Mathematics 35/89
Mixed finite-element method Raviart-Thomas mixed finite-element method In the Raviart–Thomas mixed FEM of lowest order (for triangular, tetrahedral, or regular parallelepiped grids), L 2 (Ω) is replaced by � 1 if x ∈ Ω i , � U = { p = p i χ i } , χ i ( x ) = 0 otherwise . i and H div 0 (Ω) is replaced by V = { v ∈ H div 0 (Ω) such that v | Ω i have linear components ∀ Ω i ∈ Ω , ( v · n ij ) | γ ij is constant ∀ γ ij ∈ Ω , and v · n ij is continuous across γ ij } . Applied Mathematics 36/89
Mixed finite-element method Raviart-Thomas mixed finite-element method for triangular grids Restricted to a single triangle with corner-points P 1 , P 2 , P 3 in a triangular grid, the Raviart-Thomas mixed FEM space V is spanned by basis functions ψ 1 , ψ 2 , and ψ 3 . P � x 3 � ψ 1 = α 1 y − P 1 , � x � y − P 2 ψ 2 = α 2 , P 2 � x � y − P 3 ψ 3 = α 3 . P 1 ψ 1 | E 23 · n E 23 = α 1 ( P 2 + s ( P 3 − P 2 ) − P 1 ) · n E 23 = α 1 ( P 2 − P 1 ) · n E 23 . Applied Mathematics 37/89
Mixed finite-element method Raviart-Thomas mixed finite-element method Find p = � p i χ i and v = � i ψ i such that � � � v · λ − 1 ψ j dx − p ∇ · ψ j dx = − ωg ∇ z · ψ j dx, ∀ j, Ω Ω Ω � � ∇ · v dx = q dx, ∀ j. Ω j Ω j Evaluation of integrals: For the Raviart-Thomas mixed FEM, ∇ · ψ j is cell-wise constant. Numerical integration is therefore employed only to evaluate � � v · λ − 1 ψ j dx and ωg ∇ z · ψ j dx. Ω Ω The integration is usually done using so-called Piola mappings to reference elements and numerical quadrature rules. Applied Mathematics 38/89
Are any of these methods satisfactory? TPFA: Used by most commercial reservoir simulators, but are inaccurate because they are not designed for the type of grid models that are built today using modern geomodeling tools. MPFA: More accurate than TPFA, but less robust and hard to implement for general grids, especially if the grid is non-conforming with non-matching faces. MFEM: Accurate and quite robust, but cumbersome to implement – because different cells in geological models are generally not diffeomorphic, one needs a reference element and a corresponding Piola transform for each topological case. Applied Mathematics 39/89
Grid complexity issues to consider Applied Mathematics 40/89
Mimetic finite difference methods (FDMs) Mimetic FDMs: Finite-difference-like mixed FEMs that avoid quadrature rules, reference elements, and Piola mappings: easy to implement, also for grids with irregular cell geometries and non-matching faces. Accuracy as low order MFEMs. Evaluates the integrals � � v · λ − 1 ψ j dx and ωg ∇ z · ψ j dx. Ω Ω with an approximate bilinear form. Applied Mathematics 41/89
Upscaling: Building a coarse grid model Why? Geological parameters are given on a grid that is too fine for direct simulation with current simulators — simple averaging generally gives inaccurate results. Applied Mathematics 42/89
Upscaling: Building a coarse grid model Aim of upscaling procedures for the pressure equation: Derive effective grid-block permeability tensors that produce the same flow response as the underlying geological grid-model. Let p be the solution that we obtain by solving −∇ · λ ∇ p = f, in Ω . For each coarse grid-block V we seek a tensor λ ∗ V such that � � λ ∇ p dx = λ ∗ ∇ p dx. V V V I.e., the net flow-rate ¯ v through V is related to the average v = − λ ∗ ∇ p . pressure gradient ∇ p in V through Darcy’s law ¯ Applied Mathematics 43/89
Upscaling: Building a coarse grid model The one-dimensional case Consider the one dimensional pressure equation: − ∂ x ( λp ′ ( x )) = 0 in (0 , 1) , p (0) = p 0 , p (1) = p 1 . Here the Darcy velocity v is constant and p ′ ( x ) ∝ λ − 1 . Hence, � 1 �� 1 � − 1 p ′ ( x ) = p 1 − p 0 dx p ′ ( x ) dx = p 1 − p 0 = ⇒ . λ λ 0 0 Thus, for ( a, b ) ⊂ (0 , 1) we have � b � b � − 1 � � b � 1 dx 1 � λp ′ ( x ) dx = p ′ ( x ) dx . b − a λ b − a a a a I.e., λ ∗ V is identical to the harmonic mean for all V ⊂ (0 , 1). Applied Mathematics 44/89
Upscaling: Building a coarse grid model Two special multi-dimensional cases Flow from top to bottom Flow from left to right p = 1 v . n = 0 v . n = 0 v . n = 0 p = 1 p = 0 v . n = 0 p = 0 λ ∗ = arithmetic average λ ∗ = harmonic average Conclusion: correct upscaling depends on the flow. Applied Mathematics 45/89
Upscaling: Building a coarse grid model Harmonic-arithmetic averaging Harmonic-arithmetic averaging To model flow in more than one direction, define a diagonal permeability tensor with the following diagonal components: k xx = µ z a ( µ y a ( µ x k yy = µ z a ( µ x a ( µ y k zz = µ x a ( µ y a ( µ z h )) , h )) , h )) . Here µ ξ a and µ ξ h is the arithmetic and harmonic means in the ξ -coordinate direction. Harmonic-arithmetic averaging gives correct upscaling for perfectly stratified media with flow parallel to, or perpendicular to the layers. Applied Mathematics 46/89
Numerical example Stratified model Unstratified model 5 4 3 2 1 Boundary conditions: BC1: p = 1 at ( x, y, 0), p = 0 at ( x, y, 1), else no-flow. BC2: p = 1 at (0 , 0 , 0), p = 0 at (1 , 1 , 1), else no-flow. Relative flow Stratified model Unstratified model rate error BC1 BC2 BC1 BC2 5.52e − 02 1.10e − 02 9.94e − 04 Harmonic 1 Arithmetic. 4.33e+03 2.39e+02 2.33e+04 2.13e+03 8.14e − 02 1.55e − 01 Harm. – Arit. 1 1.14 Applied Mathematics 47/89
Flow-based upscaling Flow based upscaling: For each grid block V , solve the homogeneous equation −∇ · λ ∇ p = 0 in V, with three sets of boundary conditions, one for each coordinate direction. Compute an upscaled tensor λ ∗ with components λ xξ = − Q ξ L ξ / ∆ P x , λ yξ = − Q ξ L ξ / ∆ P y , λ zξ = − Q ξ L ξ / ∆ P z . Here Q ξ , L ξ and ∆ P ξ are the net flow, the length between opposite sides, and the pressure drop in the ξ -direction. Fundamental problem: What kind of boundary conditions should be imposed? Applied Mathematics 48/89
Flow-based upscaling Fixed and periodic boundary conditions Fixed boundary conditions (yields a diagonal tensor) v.n=0 p=0 V v.n=0 V v.n=0 i i p=1 p=1 v.n=0 p=0 Periodic boundary conditions, x -direction p (1 , y, z ) = p (0 , y, z ) − ∆ p, v (1 , y, z ) = v (0 , y, z ) , p ( x, 1 , z ) = p ( x, 0 , z ) , v ( x, 1 , z ) = v ( x, 0 , z ) , p ( x, y, 1) = p ( x, y, 0) , v ( x, y, 1) = v ( x, y, 0) . Periodic BC yields a symmetric and positive definite tensor. Applied Mathematics 49/89
Numerical example Boundary conditions: BC1: p = 1 at ( x, y, 0), p = 0 at ( x, y, 1), else no-flow (flow from left to right). BC2: p = 1 at (0 , 0 , 0), p = 0 at (1 , 1 , 1), else no-flow (flow from corner to opposite corner). Relative flow Stratified model Unstratified model rate error BC1 BC2 BC1 BC2 Harm.–Arit. 1 1.14 0.081 0.155 Fixed BC 1 1.14 1 1.893 Periodic BC 1 1.14 0.986 1.867 Applied Mathematics 50/89
Remarks on upscaling Upscaling by design of a coarse simulation model has been, and still is, the default (and only) way that industry employs to fit the simulation model to the capabilities of commercial reservoir simulators. But ... Upscaling depends on flow regime. The favorite and most robust upscaling methods are designed for grids with shoe-box shaped grid blocks. Upscaling is inadequate for modeling flow in carbonate reservoirs, channelized reservoirs, or reservoirs where fine scale features that are not resolved by the coarse grid have dominant impact on large scale flow patterns. Applied Mathematics 51/89
Examples of heterogeneous structures Logarithm of permeability in a shallow-marine Tarbert formation. Logarithm of permeability in a fluvial Upper Ness formation. Applied Mathematics 52/89
Exhale – inhale Multiscale methods coming up. Applied Mathematics 53/89
Multiscale methods for the pressure equation Multiscale methods: Numerical methods that attempt to model physical phenomena on coarse grids while honoring small-scale features that impact the coarse grid solution in an appropriate way. Heterogeneous Multiscale Methods Local global upscaling Multiscale discontinuous Galerkin Methods t s d e n d e t e o i m h z Two−scale locally conservative upscaling n e t i i l e l f e m a r e m Multiscale mixed finite element method e Multiscale finite element methods G Residual free bubbles Variational multiscale methods Multiscale finite volume method Applied Mathematics 54/89
Flow based upscaling versus multiscale method Standard upscaling: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: ⇓ Coarse grid blocks: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: ⇓ Coarse grid blocks: ⇓ Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: ⇓ Coarse grid blocks: ⇓ ⇑ Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: ⇓ ⇑ Coarse grid blocks: ⇓ ⇑ Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: ⇓ ⇑ Coarse grid blocks: ⇓ ⇑ Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: Multiscale method: ⇓ ⇑ Coarse grid blocks: ⇓ ⇑ Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: Multiscale method: ⇓ ⇓ ⇑ Coarse grid blocks: Coarse grid blocks: ⇓ ⇑ ⇓ Flow problems: Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: Multiscale method: ⇓ ⇓ ⇑ Coarse grid blocks: Coarse grid blocks: ⇓ ⇑ ⇓ ⇑ Flow problems: Flow problems: Applied Mathematics 55/89
Flow based upscaling versus multiscale method Standard upscaling: Multiscale method: ⇓ ⇑ ⇓ ⇑ Coarse grid blocks: Coarse grid blocks: ⇓ ⇑ ⇓ ⇑ Flow problems: Flow problems: Applied Mathematics 55/89
Multiscale finite-element method (MsFEM), 1D Model problem: Consider the following elliptic problem: ∂ x ( k ( x ) ∂ x p ) = f, in Ω = (0 , 1) , p (0) = p (1) = 0 , where f, k ∈ L 2 (Ω) and 0 < α < k ( x ) < β for all x ∈ Ω. Variational formulation: Find p ∈ H 1 0 (Ω) such that for all v ∈ H 1 a ( p, v ) = ( f, v ) 0 (Ω) , (1) where ( · , · ) is the L 2 inner-product and � k ( x ) p ′ ( x ) v ′ ( x ) dx. a ( p, v ) = Ω Applied Mathematics 56/89
Multiscale finite-element method (MsFEM), 1D Let N K = { 0 = x 0 < x 1 < . . . < x n − 1 < x n = 1 } be a set of nodal points and define K i = ( x i − 1 , x i ). For i = 1 , . . . , n − 1 define a basis function φ i ∈ H 1 0 (Ω) by a ( φ i , v ) = 0 for all v ∈ H 1 0 ( K i ∪ K i +1 ) , φ i ( x j ) = δ ij , where δ ij is the Kronecker delta. Basis functions: 0 1 = Linear FEM basis functions = Multiscale FEM basis functions Applied Mathematics 57/89
Multiscale finite-element method (MsFEM), 1D MsMFEM: Find the unique function p 0 in V ms = span { φ i } { u ∈ H 1 0 (Ω) : a ( u, v ) = 0 for all v ∈ H 1 = 0 ( ∪ i K i ) } satisfying for all v ∈ V ms . a ( p 0 , v ) = ( f, v ) Theorem Assume p solves the variational formulation. Then p = p 0 + � n i =1 p i where p i ∈ H 1 0 ( K i ) is defined by for all v ∈ H 1 a ( p i , v ) = ( f, v ) 0 ( K i ) . Applied Mathematics 58/89
Proof I of super convergence property Galerkin projection property: Assume p solves the variational formulation and v ∈ V ms . Then a ( p − p 0 , v ) = a ( p, v ) − a ( p 0 , v ) = ( f, v ) − ( f, v ) = 0 . Thus, p 0 is the orthogonal projection of p onto V ms . 0 (Ω) = V ms ⊕ H 1 Since H 1 0 ( ∪ i K i ) it follows that p 0 ( x i ) = p ( x i ) for all i, i.e., that p 0 is the interpolant of p in V ms . Applied Mathematics 59/89
Proof II of super convergence property Let p I be the interpolant of p in V ms . Then p − p I ∈ H 1 0 ( ∪ i K i ) and it follows from the mutual orthogonality of V ms and H 1 0 ( ∪ i K i ) with respect to a ( · , · ) that for all v ∈ V ms . a ( p − p I , v ) = 0 Hence, for all v ∈ V ms , a ( p I , v ) = a ( p, v ) = ( f, v ) = a ( p 0 , v ) and we see that a ( p I − p 0 , v ) = 0 for all v ∈ V ms . Thus, in particular, by choosing v = p I − p 0 we obtain a ( p I − p 0 , p I − p 0 ) = 0 , which implies p 0 = p I . Applied Mathematics 60/89
Multiscale finite element method, 1D Super convergence property: Solution of the variational problem is decomposed into the MsFEM solution and solutions of independent local subgrid problems. MsFEM in 1D = a Schur complement decomposition. Applied Mathematics 61/89
Multiscale finite element method, 1D Super convergence property: Solution of the variational problem is decomposed into the MsFEM solution and solutions of independent local subgrid problems. MsFEM in 1D = a Schur complement decomposition. Does the result extend to higher dimensions? Applied Mathematics 61/89
Multiscale finite element method, 1D Super convergence property: Solution of the variational problem is decomposed into the MsFEM solution and solutions of independent local subgrid problems. MsFEM in 1D = a Schur complement decomposition. Does the result extend to higher dimensions? No, but the basic construction applies and helps us understand how subgrid features of the solution can be embodied into a coarse grid approximation space. Applied Mathematics 61/89
Multiscale finite element method, 2D MsMFEM (for −∇ · K ( x ) ∇ p = f ): Find the unique function p 0 in V ms = span { φ i,j } { u ∈ H 1 0 (Ω) : a ( u, v ) = 0 for all v ∈ H 1 = 0 ( ∪ i K i ) } satisfying for all v ∈ V ms . a ( p 0 , v ) = ( f, v ) Here ( · , · ) is the L 2 inner-product and � a ( p, v ) = ∇ p · K ( x ) ∇ v dx. Ω Applied Mathematics 62/89
Multiscale finite element method, 2D Definition of basis functions MsFEM basis functions in 2D x x x i−1,j+1 i,j+1 i+1,j+1 p ∈ V ms implies that ∇· K ∇ φ i,j = 0 in all Ω m . Ω 1 Ω 2 x x x i−1,j i,j i+1,j Ω Ω 3 4 x x x i−1,j−1 i,j−1 i+1,j−1 Note: If p ∈ V ms then p 0 = p by projection property. Applied Mathematics 63/89
Multiscale finite element method, 2D Definition of basis functions MsFEM basis functions in 2D x x x i−1,j+1 i,j+1 i+1,j+1 p ∈ V ms implies that ∇· K ∇ φ i,j = 0 in all Ω m . Ω 1 Ω 2 φ i,j = 0 on edges not x x x i−1,j i,j i+1,j emanating from x i,j . Ω Ω 3 4 x x x i−1,j−1 i,j−1 i+1,j−1 Note: If p ∈ V ms then p 0 = p by projection property. Applied Mathematics 63/89
Multiscale finite element method, 2D Definition of basis functions MsFEM basis functions in 2D x x x i−1,j+1 i,j+1 i+1,j+1 p ∈ V ms implies that ∇· K ∇ φ i,j = 0 in all Ω m . Ω 1 Ω 2 φ i,j = 0 on edges not x x x i−1,j i,j i+1,j emanating from x i,j . φ i,j ( x m,n ) = δ i,m δ j,n . Ω Ω 3 4 x x x i−1,j−1 i,j−1 i+1,j−1 Note: If p ∈ V ms then p 0 = p by projection property. Applied Mathematics 63/89
Multiscale finite element method, 2D Definition of basis functions MsFEM basis functions in 2D x x x i−1,j+1 i,j+1 i+1,j+1 p ∈ V ms implies that ∇· K ∇ φ i,j = 0 in all Ω m . Ω 1 Ω 2 φ i,j = 0 on edges not x x x i−1,j i,j i+1,j emanating from x i,j . φ i,j ( x m,n ) = δ i,m δ j,n . Boundary conditions Ω Ω 3 4 on edges emanating from x i,j ? x x x i−1,j−1 i,j−1 i+1,j−1 Note: If p ∈ V ms then p 0 = p by projection property. Applied Mathematics 63/89
Multiscale finite volume method, 2D Multiscale finite volume method (MsFVM) is a control-volume finite-element version of the MsFEM. A CVFE method seeks solutions in a finite-element space (on a dual-grid), but rather than formulating the problem in a variational framework, it employs instead a finite-volume formulation (on a primal grid). MsFEM basis function on a dual grid connecting cell centers of the primal grid. Boundary conditions on “interior edges” derived by solving reduced dimen- sional pressure equation. Applied Mathematics 64/89
Multiscale finite volume method (for −∇ · K ∇ p = f ) Step 1: Compute basis functions φ i ( φ i denotes the basis function associated with the center of Ω i ). Step 2: Find p = � i p i φ i such that � � − K ∇ p · n ds = f dx. ∂ Ω j Ω j Step 3: Downscale to obtain a mass conservative velocity field v = − K ∇ u in all Ω j , ∇ · v = f in all Ω j , − K ∇ p v = on all ∂ Ω j . Applied Mathematics 65/89
Multiscale finite volume method Complexity and robustness issues Complexity and robustness issues for MsFVM Inaccurate for high-aspect ratios or cases with large anisotropy. Difficult to implement for complex grids, e.g., it is difficult to define a dual grid if primal grid honors geological features. No escaping downscaling for multi-phase flows; computational cost comparable to a domain decomposition sweep. Accuracy depends highly on proper boundary conditions. Applied Mathematics 66/89
Multiscale mixed finite element method (MsMFEM) Pressure eq.: v = − λ ( ∇ p + ωg ∇ z ) and ∇ · v = q . Mixed formulation (no-flow boundary conditions): Find p ∈ U ⊂ L 2 (Ω) and v ∈ V ⊂ H div 0 (Ω) such that � � � v · λ − 1 u dx − p ∇ · u dx = − ωg ∇ z · u dx, Ω Ω Ω � � l ∇ · v dx = ql dx, Ω Ω for all u ∈ V and l ∈ U . Multiscale mixed finite element method (MsMFEM): Mixed finite element method for which V is designed to embody the impact of fine scale structures. Applied Mathematics 67/89
Multiscale mixed finite element method Basis functions Associate a basis function χ m for pressure with each grid block K : � 1 if x ∈ K m , U = span { χ m : K m ∈ K} where χ m = 0 else, and a velocity basis function ψ ij with each interface ∂K i ∩ ∂K j : span { ψ ij = − k ∇ φ ij } V = ψ ij · n = 0 on ∂ ( K i ∪ K j ) � w i ( x ) in K i , ∇ · ψ ij = − w j ( x ) in K j . where w i is a weight function � with K i w i ( x ) dx =1. Applied Mathematics 68/89
Multiscale mixed finite element method Coarse grid Each coarse grid block is a connected set of cells from geomodel. Example: Coarse grid obtained with uniform coarsening in index space. Grid adaptivity at well locations: One block assigned to each cell in geomodel with well perforation. Applied Mathematics 69/89
Multiscale mixed finite element method Workflow At initial time Compute ψ for each domain Detect all adjacent blocks For each time-step: Assemble and solve coarse grid system. Recover fine grid velocity: v = � ij v ij ψ ij . Solve mass balance equations. Applied Mathematics 70/89
Multiscale mixed finite element method Layer 36 from SPE10 model 2 (Christie and Blunt, 2001). Example: Layer 36 from SPE10 (Christie and Blunt, 2001). Primary features Coarse pressure solution, subgrid resolution at well locations. Coarse velocity solution with subgrid resolution everywhere. Applied Mathematics 71/89
Numerical example Layer 85 from SPE10 model 2 (Christie and Blunt, 2001). Water injected at center of Pressure computed on fine grid an oil-filled reservoir. 350 0.9 Producer at each corner. 300 0.8 250 0.7 0.6 Geomodel: 60 × 220 grid. 200 0.5 150 0.4 Figures: water fraction 0.3 100 0.2 when volume injected = 50 0.1 0 30% of total flow volume. 0 100 200 300 400 500 600 MsFVM: coarse grid = 10 × 44 MsMFEM: coarse grid = 10 × 44 350 350 0.9 0.9 300 300 0.8 0.8 0.7 0.7 250 250 0.6 0.6 200 200 0.5 0.5 150 150 0.4 0.4 0.3 0.3 100 100 0.2 0.2 50 50 0.1 0.1 0 0 0 100 200 300 400 500 600 0 100 200 300 400 500 600 Applied Mathematics 72/89
11.20 11.25 11.30 Applied Mathematics 73/89
Implementational issues for the MsMFEM Coarse grid generation Choice of weighting function in definition of basis functions. Incorporation of global information in basis functions. Choice of numerical method for computation of basis functions Assembly of linear system Applied Mathematics 74/89
Coarse grid generation MsMFEM (unique) grid flexibility: Given a method applicable to solving the local flow problems on the subgrid, the MsMFEM can be formulated on any coarse grid where each grid block consists of a connected collection fine-grid cells. Applied Mathematics 75/89
Coarse grid generation MsMFEM (unique) grid flexibility: Given a method applicable to solving the local flow problems on the subgrid, the MsMFEM can be formulated on any coarse grid where each grid block consists of a connected collection fine-grid cells. Applied Mathematics 75/89
Coarse grid generation MsMFEM (unique) grid flexibility: Given a method applicable to solving the local flow problems on the subgrid, the MsMFEM can be formulated on any coarse grid where each grid block consists of a connected collection fine-grid cells. Applied Mathematics 75/89
Coarse grid generation MsMFEM (unique) grid flexibility: Given a method applicable to solving the local flow problems on the subgrid, the MsMFEM can be formulated on any coarse grid where each grid block consists of a connected collection fine-grid cells. The process of generating a coarse simulation grid from a complex geological model can be greatly simplified, especially when the fine grid is fully unstructured or has geometrical complications. Applied Mathematics 75/89
Coarse grid generation Limitations: Limited ability to model flow around large scale flow barriers Applied Mathematics 76/89
Coarse grid generation Limitations: Limited ability to model flow around large scale flow barriers Limited ability to model bi-directional flow across interfaces. Applied Mathematics 76/89
Coarse grid generation Limitations: Limited ability to model flow around large scale flow barriers Limited ability to model bi-directional flow across interfaces. Applied Mathematics 76/89
Coarse grid generation Guidelines Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction 6 7 8 6 7 8 5 5 2 1 3 1 3 4 2 Guidelines: 1: The coarse grid should minimize the occurrence of bi-directional flow across interfaces. Grid structures that increase the likelihood for bidirectional flow are: Coarse-grid faces with (highly) irregular shapes. Blocks that have only one neighbor. Blocks with no faces transverse to the major flow directions. Applied Mathematics 77/89
Coarse grid generation Guidelines Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction Flow direction 6 7 8 6 7 8 5 5 2 1 3 1 3 4 2 Guidelines: 2: Blocks should follow geological layers when possible. Avoids disrupting flow patterns. Applied Mathematics 77/89
Coarse grid generation Guidelines Guidelines: 3: Blocks in the coarse-grid should adapt to flow barriers. Applied Mathematics 77/89
Recommend
More recommend
