Multiscale Extended Finite Element Method for Deformable Fractured Porous Media
MMultiscale Extended Finite Element Methodfor Deformable Fractured Porous Media
Fanxiang Xu a , ∗ Hadi Hajibeygi b Lambertus J. Sluys a a Materials, Mechanics, Management and Design, Faculty of Civil Engineering andGeosciences, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, TheNetherlands b Department of Geoscience and Engineering, Faculty of Civil Engineering andGeosciences, Delft University of Technology, P.O. Box 5048, 2600 GA Delft, TheNetherlands
Abstract
Deformable fractured porous media appear in many geoscience applications. Whilethe extended finite element (XFEM) has been successfully developed within thecomputational mechanics community for accurate modeling of the deformation,its application in natural geoscientific applications is not straightforward. This ismainly due to the fact that subsurface formations are heterogeneous and span largelength scales with many fractures at different scales. In this work, we propose anovel multiscale formulation for XFEM, based on locally computed basis functions.The local multiscale basis functions capture the heterogeneity and discontinuitiesintroduced by fractures. Local boundary conditions are set to follow a reduced-dimensional system, in order to preserve the accuracy of the basis functions. Us-ing these multiscale bases, a multiscale coarse-scale system is then governed al-gebraically and solved, in which no enrichment due to the fractures exist. Suchformulation allows for significant computational cost reduction, at the same time, itpreserves the accuracy of the discrete displacement vector space. The coarse-scalesolution is finally interpolated back to the fine scale system, using the same multi-scale basis functions. The proposed multiscale XFEM (MS-XFEM) is also integratedwithin a two-stage algebraic iterative solver, through which error reduction to anydesired level can be achieved. Several proof-of-concept numerical tests are presentedto assess the performance of the developed method. It is shown that the MS-XFEMis accurate, when compared with the fine-scale reference XFEM solutions. At thesame time, it is significantly more efficient than the XFEM on fine-scale resolution.As such, it develops the first scalable XFEM method for large-scale heavily fracturedporous media.
Key words:
Fractured porous media, extended finite element, multiscale,geomechanics, scalable iterative solver,
Preprint submitted to Elsevier Science August 7, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] A ug Introduction
Subsurface geological formations are often highly heterogeneous and heavilyfractured at multiple scales. Heterogeneity of the deformation properties (e.g.elasticity coefficients) can be of several orders of magnitude which occurs atfine scale (cm) resolution. The reservoirs also span large scales, in the order ofkms. Numerical simulation of mechanical deformation for such complex sys-tems is necessary to optimise the geo-engineering operations [1, 2], and assesstheir safety and manage the associated risks (e.g. fracture propagation, faultslip and induced seismicity). Though being crucially important, simulation ofthese systems are beyond the scope of classical numerical schemes.Presence of highly heterogeneous coefficients with high resolution within large-scale domains has been systematically addressed in the computational geo-science community through the development of multiscale finite element andfinite volume methods [3–6]. Recent developments also include mechanicaldeformation coupled with fluid pore pressure dynamics [7–11]. In presenceof many fractures, however, the complexity of the computational model in-creases significantly. As such, development of a robust multiscale strategy fordeformation of heavily fractured porous media, which also allows for conver-gent systematic error reduction [12–14], is of high interest in the geosciencecommunity.The presence of fractures within the computational domain can be includedexplicitly by two approaches of (1) unstructured grid and (2) immersed orembedded methods.The unstructured grid approach [15–17] generates a discrete computationaldomain in which fractures are always at the interfaces of elements. This al-lows for convenient treatment of their effect, however, at the cost of complexmeshing. The complex mesh generation for three-dimensional (3D) large scaledomains with many fractures is challenging, specially when fractures dynami-cally extend their geometries. On the other hand, the immersed or embeddedapproach allows for independent grids for matrix block and fractures, by intro-ducing enrichment of the discrete connectivity (for flow) and shape functions(for mechanics) [18–22]. These enriched formulations are aimed at representingdiscontinuities within the overlapping matrix cell, without any adjustment norrefinement of the grid [23]. The enrichment strategy for modeling deformationusing finite-element schemes in presence of fractured media are referred to as’extended finite element (XFEM)’ methods. ∗ Corresponding author.
Email addresses: [email protected] (Fanxiang Xu),
[email protected] (Hadi Hajibeygi),
[email protected] (Lambertus J. Sluys).
Consider the domain Ω bounded by Γ as shown in figure 1. Prescribed dis-placements or Dirichlet boundary condition are imposed on Γ u , while tractionsare imposed on Γ t . The crack surface Γ c (lines in 2-D and surfaces in 3-D) isassumed to be traction-free.Figure 1: An illustration of fractured domain setupThe momentum balance equations and boundary conditions read ∇ · σ + f = 0 in Ω (1) σ · −→ n = ¯ t on Γ t (2) σ · −→ n = 0 on Γ c (3) u = ¯ u on Γ u , (4)where σ is the stress tensor and u is the displacement field over the whole4omain. −→ n is the normal vector pointing outside the domain [34, 35].The constitutive law with linear elasticity assumption reads σ = C : ε = C : ∇ s u (5)where, ∇ s denotes the symmetrical operator and C is the property tensordefined as C = λ + 2 µ µ µ λ + 2 µ
00 0 µ , with λ and µ denoting the Lame’s parameters [36, 37].The strain tensor ε is expressed as ε = ∇ s u = 12 ( ∇ u + ∇ T u ) (6)where, ∇ denotes the gradient operator.Substituting Eqs. (5) and (6) in the governing equation (1) results in a 2ndorder Partial Differential Equation (PDE) for displacement field u ∇ · ( C : ∇ s u ) + f = 0 . (7)Equation (7) is then solved for computational domains with cracks (represent-ing faults and fractures). This is done by the extended finite element (XFEM)method, which is briefly revisited in the next section. The FEM with smooth shape functions N i provides an approximate numericalsolution to Eq. (7) for displacement unknowns, i.e., u = (cid:88) i ∈ I u i N i . (8)This formula can be used for computational domains without discontinuity.The FEM approximation is insufficient to capture discontinuities imposed bythe existence of the fractures and faults. As such, the XFEM method intro-duces two sets of enrichment to the original FEM in order to allow it to capturethe discontinuities without adapting the grid. These enrichment sets are asso-ciated with the body and tip of the fractures and faults. The body is enriched5y jump functions, and the tip by tip enrichment functions [25]. Below briefdescriptions of these two enrichment functions are provided. The jump enrichment represents the discontinuity involved in the displacementfield across the fracture and fault main body. The jump enrichment is oftenchosen as the step or Heaviside function, which can be expressed as H ( x ) = + − − .Note that Ω + and Ω − zones are determined based on the normal vector point-ing out of the fracture curve. For line fractures, the direction can be any side,as long as all discrete elements use the same + and - sides for a fracture. The tip enrichment represents the discontinuity of the displacement field nearthe fracture tip. This type of enrichment function, denoted by F l , is basedon the auxiliary displacement field near the fracture tip and contains fourfunctions, i.e., F l ( r, θ ) = {√ rsin ( θ , √ rcos ( θ , √ rsin ( θ sin ( θ ) , √ rcos ( θ sin ( θ } . (9)These four functions around the fracture tip inside the element are plottedin Figure 2. The red segment, shown on the base of the plots, represents thefracture which ends in the element. Note that only the √ rsin ( θ ) contains adiscontinuity around the fracture tip, while other functions are smooth.6igure 2: Four types of tip enrichment functions inside the element. The redsegment represents the crack with its tip located at center point (0,0). Thediscontinuity can be seen clearly in the top left function, √ rsin ( θ ). To decide whether the node is enriched or not, the node location related to thefracture is the key factor. The sketch of the enrichment mechanism is shownin Figure 3. More precisely, in this figure, the tip and jump enriched nodes arehighlighted in red and black, respectively.
Figure 3: Enrichment mechanism: node I and J will be enriched using tip andjump functions. 7 .2 XFEM linear system
The XFEM approximates the continuum displacement field u at fine-scalemesh resolution h by u h which is defined as u ≈ u h = (cid:88) i ∈ Ω h u i N i + (cid:88) j ∈ J a j N j H ( x ) + (cid:88) k ∈ K N k , (cid:104) (cid:88) l =1 F l ( x ) b lk (cid:105) (10)where N , H and F l represent, respectively, the classical FEM shape functions,the Heaviside function and the tip enrichment functions. The fine-scale meshhas Ω h nodes. Moreover, u denotes the standard degrees of freedom (DOFs)associated to the classical finite element method. a denotes the extra DOFsassociated to the jump enriched node. For the jump enriched nodes, in the2D domains, each node would contain 2 extra DOFs. Furthermore, b indicatesthe extra DOFs associated to the tip enrichment, which adds four extra DOFsper direction (8 in total in a 2D domain) for each tip inside an element.The first term in the right-hand-side (RHS) of Eq. (10) is the contributionof the classical finite element method. This term captures the smooth defor-mation, using classical shape functions. The second term, however, representsthe contribution of the jump enrichment. Note that the jump enrichment ismodeled by the weighted Heaviside functions, with weights being the classi-cal shape functions. There will be as many jump enrichment functions as thenumber of fractures inside an element. Finally, the third term in the RHS isthe contribution of the fracture tips. Note that if several fracture tips end upin an element, there will be 4 additional DOFs per tip per direction in thatelement.The resulting linear system entails the nodal displacement unknowns u , as wellas the jump level a and tip weight b per fracture (and fault). The augmentedXFEM linear system K h d h = f h , therefore, reads K uu K ua K ub K au K aa K ab K bu K ba K bb (cid:124) (cid:123)(cid:122) (cid:125) K h uab (cid:124) (cid:123)(cid:122) (cid:125) d h = f u f a f b (cid:124) (cid:123)(cid:122) (cid:125) f h . (11)Compared to the classical FEM, there exist several additional blocks involvedin the stiffness matrix, due to the existence of the discontinuities. The advan-tage of XFEM is that it does not rely on complex mesh geometry, instead, itallows fractures to overlap with the matrix elements. On the other hand, forgeoscientific fractured systems, the additional DOFs due to the enrichmentprocedure results in excessive computational costs. This imposes a significant8hallenge for the XFEM application in geoscience applications. In this paper,we develop a scalable multiscale procedure which constructs a coarse-scalesystem based on locally supported basis functions. The method is describedin the next section. A multiscale formulation provides an approximate solution u (cid:48) h to the fine-scaleXFEM deformation u h through u h ≈ u (cid:48) h = (cid:88) i ∈ Ω H N Hi u Hi , (12)where N Hi are the coarse-scale (multiscale) basis functions and u Hi are thecoarse-scale nodal displacements at coarse mesh Ω H . Note that this multiscaleformulation does not include any enrichment functions. Instead, all enrichmentfunctions are incorporated in the construction of accurate coarse-scale basisfunctions N H . This allows for significant computational complexity reduction,and makes the entire formulation attractive for field-scale geoscientific appli-cations.Next, construction of the coarse-scale system and the basis functions will bepresented. MS-XFEM solves the linear deformation system on a coarse mesh, imposed ona given fine-scale mesh, as shown in figure 4. The coarsening ratio is definedas the ratio between the coarse mesh size and fine-scale mesh size. × u h ≈ u (cid:48) h = P d H , (13)where P is the matrix of basis functions (i.e., prolongation operator) and d H is the coarse-scale deformation vector for u H unknowns. Algebraic formulationallows for convenient implementation of the proposed MS-XFEM, and its in-tegration as a black-box tool for any given classical XFEM solver. Therefore,the remainder of the article will be devoted to the formulation.The coarse-scale solution d H needs to be found by solving a coarse-scale sys-tem. To construct the coarse-scale system and solve for d H , one has to restrict(map) the fine-scale linear system( K h d h = f h ) to the coarse-scale, i.e.,( R K h P ) (cid:124) (cid:123)(cid:122) (cid:125) K H d H = R f h . (14)Here, R is the restriction operator with the size of Ω H × Ω h + j + t , where Ω h + j + t is the size of the fine-scale enriched XFEM system including jump and tipenrichment. Prolongation operator P has the dimension of Ω h + j + t × Ω H . Thisresults in the coarse-scale system matrix K H size of Ω H × Ω H .The finite-element-based restriction function is introduced as the transpose ofthe prolongation matrix, i.e., R = P T . (15)Therefore, the coarse-scale matrix K H is symmetric-positive-definite (SPD),if K h is SPD. Once the coarse-scale system is solved on Ω H space for d H ,one can find the approximate fine-scale solution using Eq. (13). Overall, themultiscale procedure can be summarised as finding an approximate solution d (cid:48) h according to d h ≈ d (cid:48) h = P d H = P ( R K h P ) − R f h . (16)Next, the prolongation operator P , i.e., the basis functions are explained indetail. Once P is known, all terms in Eq. (16) are defined. To obtain the basis functions, the governing equation (7) without any sourceterm using XFEM, i.e., (11) needs to be solved in each coarse element ω H .10his can be expressed as solving ∇ · ( C : ( ∇ S N Hi )) = 0 in Ω H , (17)subject to local boundary conditions. Here, we develop a reduced-dimensionalequilibrium equation to solve for the boundary cells [5, 38], i.e., ∇ (cid:107) · ( C r : ( ∇ S (cid:107) N Hi )) = 0 on Γ H . (18)Here, Γ H denotes the boundary cells of the coarse element Ω H . In addition, ∇ S (cid:107) denotes the reduced dimensional divergence and symmetrical gradient op-erators, which act parallel to the direction of the local domain boundary. For2D geometries, the reduced-dimensional boundary condition represents the 1D(rod) deformation model along the coarse element edges. Note that the localbasis functions involve transverse equilibrium, i.e., therefore the prolongationmatrix P reads P = P xx P xy P yx P yy . (19) (a) (b)Figure 5. Illustration of the multiscale local basis functions, constructed usingXFEM for the node H in x direction (a) and y direction (b). Figure 5 shows an example of a local system to be solved for basis functionsbelonging to the highlighted node H in x and y directions. Note that theDirichlet value of 1 is set at H for each directional basis functions, while allother 3 coarse mesh nodes are set to 0.Note that, as shown in Fig. 5, the boundary problem is solved for both edges11hich have the node H at one of their end values. More precisely, e.g., to findthe basis function in x-direction for node H , we set the value of u x ( H ) = 1at the location H . This causes extension of the horizontal boundary cells andbending of the vertical boundary.Once the boundary values are found, the internal cells are solved subjectedto Dirichlet values for the boundary cells. An illustration of a basis functionobtained using this algorithm is presented in figure 6. Figure 6: Illustration of a basis function that captures the discontinuity of afracture. Yellow segment represents the fracture.Note that the illustrated basis function captures the fractures, because ofthe XFEM enrichment procedure. The basis function N Hi will be stored inthe column i of the prolongation operator P . Once all basis functions arefound, the operator P is also known and one can proceed with the multiscaleprocedure as explained before.Next, we explain how the basis functions can be algebraically computed basedon the given XFEM fine-scale system. This crucial step allows for convenientintegration of our multiscale method into a given XFEM simulator. The basis function formulation (17) subjected to the local boundary condition(18) can be constructed and solved purely algebraically. This is important,since it allows for convenient integration of the devised multiscale methodinto any existing XFEM simulator. 12onsider the coarse cell (local domain) as shown in figure 7. The cells aresplit into 3 categories of internal, edge and vertex (node), depending on theirlocations [13].
Figure 7: Illustration of the 3 categories of Internal, Edge, and Vertex cells,corresponding to the position of each fine cell within a coarse element.Note that the vertex nodes are in fact the coarse mesh nodes, where the coarse-scale solution will be computed. The basis functions are needed to interpolatethe solution between the vertex cells through the edge and internal cells.To develop the basis functions, first the fine-scale stiffness matrix K h is per-muted, such that the terms for vertex, then edge and finally the internal cellsappear. The permutation operator T as such reorders K h into K v such that˘ K v = T K h T T = T K uu K ua K ub K au K aa K ab K bu K ba K bb T T = K II K IE K IV K EI K EE K EV K V I K V E K V V . (20)Here, I represents the internal nodes, E represents the edge nodes and V represents the vertex nodes. The permuted linear system, therefore, reads K II K IE K IV K EI K EE K EV K V I K V E K V V d I d E d V = f I f E f V (21)Note that the permuted system collects all entries of the XFEM discrete sys-tem belonging to I, E, and V cells. Therefore, the XFEM enrichment entries13ue to tips and jumps are within their corresponding I, E, and V entries.The reduced-dimensional boundary condition is now being imposed by replac-ing the 2D equation for E by a 1D XFEM discrete system. This leads the entry¯ K V I to vanish, as there will be no connectivity between the edge and internalcells for the edge cells. This 1D edge equations can then be expressed as K REE d E + K REV d V = 0 . (22)Knowing that the solution at vertex cells will be obtained from the coarse-scalesystem, the reorder fine-scale system matrix can now be reduced to K II K IE K IV K REE K REV I V V d (cid:48) I d (cid:48) E d (cid:48) V = . (23)Note that the equations for basis functions do not have any source terms intheir right-hand-side. The upper-triangular matrix of Eq. (23) can be easilyinverted to give the prolongation operator, i.e., given the coarse nodes solutions d (cid:48) V , one can obtain the solution at the edge via d (cid:48) E = − ( K REE ) − K REV d (cid:48) V = P E d (cid:48) V . (24)similarly, the solution at the internal cells reads d (cid:48) I = − K − II ( K IE d (cid:48) E + K IV d (cid:48) V )= − K − II ( − K IE ( K REE ) − K REV + K IV ) d (cid:48) V = P I d (cid:48) V . (25)Note that P E and P I are the sub-matrices of the prolongation operator, i.e., d (cid:48) = d (cid:48) I d (cid:48) E d (cid:48) V = − K − II ( − K IE ( K REE ) − K REV + K IV ) − ( K REE ) − K REV I V V (cid:124) (cid:123)(cid:122) (cid:125) P d (cid:48) V . (26)Here, I V V is the diagonal identity matrix equal to the size of the vertex nodes.After defining the prolongation operator algebraically, based on the entries ofthe 2D XFEM (for internal cells) and 1D XFEM (for edge cells), one can findthe multiscale solution. 14 .4 Iterative multiscale procedure (iMS-XFEM)
The multiscale solutions with the accurate XFEM basis functions can be usedto provide an approximate efficient solution for many practical applications.However, it is important to control the error and reduce it to any desiredtolerance [12] if needed. As such, the MS-XFEM is paired with a fine-scalesmoother (here, ILU(0) [32]) to allow for error reduction. Note that this itera-tive procedure can also be used within a GMRES iterative loop [39] to enhanceconvergence rates. The study of the most efficient iterative strategy to reducethe error is outside the scope of the current paper. The iterative procedurethen reads • Construct the P and R operators • Iterate until || r ν +1 || = || f h − K h d (cid:48) ν +1 || < e r
1. MS-XFEM stage: δd (cid:48) ν +1 / = P ( R K h P ) − R r ν
2. Smoothing stage (apply n s times ILU(0)): δd (cid:48) ν +1 = ( M ns ILU(0) ) − r ν +1 / Note that n s is defined by user. In this section several test cases are considered to investigate the performanceof MS-XFEM both as approximate solver and integrated within the iterativeerror reduction procedure.
In the first test, a square 2D domain of L × L with L = 10[m] is considered,which contains a single horizontal fracture in its centre, as shown in figure 8a.The fine-scale mesh consists of 40 ×
40 cells, while the MS-XFEM containsonly 5 × q = 5 × [N/m] magnitude. 15 (a) (b)Figure 8. Test case 1 :(a) illustration of the model setup, (b) heterogeneous Young’smodulus distribution. Note the units are SI. Results are shown in figure 9. The black lines on figure 9 (b) and (c) representthe coarse scale mesh. It is clear that the results of MS-XFEM on only 5 gridcells is in reasonable agreement with that of the fine-scale XFEM solver usinga 40 ×
40 mesh. Note that no enrichment for the MS-XFEM is used, and thebasis functions are computed using the XFEM method on local domains. (a)(b) || e y || = 2 . × − (c) || e y || = 1 . × − Figure 9. Test Case 1: displacement field for a heterogeneous fractured reservoirusing (a) fine scale XFEM and (b) MS-XFEM (c) iMS-XFEM after 3 iterations.
16 basis function for a fractured local domain is illustrated in figure 10. Notethat the discontinuity is captured by the basis functions, since XFEM is usedto solve for it. (a) (b)Figure 10. Basis functions of single fracture test case. Single discontinuity is capturedby axial equilibrium and transverse equilibrium solutions
The effect of the coarsening ratio is shown in figure 11. The error e in thisfigure is computed using e i = || u i,MS − u i,f || N , ∀ i ∈ x, y, (27)where N is the number of fine-scale mesh nodes. u i,MS and u i,f denote theproloned MS-XFEM solution displacement field and fine-scale solution dis-placement field, respectively. coarsening ratio -3.9-3.8-3.7-3.6-3.5-3.4-3.3-3.2-3.1-3-2.9 l og10 ( e ) e x e y Figure 11: Change of errors with different coarsening ratiosThe MS-XFEM errors are due to the local boundary conditions used to calcu-late the basis functions, and also because no additional enrichment functionsare imposed at coarse scale. Note this means that for heterogeneous domainsthere can be a finer resolution for coarse cells at which the local boundaryconditions impose more errors compared with coarser resolutions. In spite of17his, figure 11 clearly shows, for this example, a decaying trend of the errorwith respect to the finer coarse mesh is observed.
As discussed in section 3.4, one can pair the MS-XFEM in an iterative strat-egy in which the error is reduced to any desired level [13]. From figure 9 (c)that with 3 times of fine scale smoothers applied in the second stage after 3iterations the MS-XFEM result has been improved compared to the fine-scaleresult and the error is decreased significantly. Results of the iterative MS-XFEM procedure (iMS-XFEM) are shown in figure 12. Different smoothingsteps per iteration values n s are used. Note than neither GMRES [39] nor anyother iterative convergence enhancing procedure is used here. Clearly, one canreduce the multiscale errors to the machine accuracy by applying iMS-XFEMiterations. In particular, for practical applications, one can stop iterations af-ter a few counts, once the error norm is below the level of uncertainty τ withinthe parameters of the problem. e i (cid:54) τ, i = x, y (28)In which τ is chosen as 10 − in here.In figure 12 it is shown that convergence is achieved with n s rounds of the finescale smoother involved in the second stage. iteration -11-10-9-8-7-6-5-4-3 l og10 ( e x ) n s =1n s =3n s =5n s =9 (a) iteration -11-10-9-8-7-6-5-4-3 l og10 ( e y ) n s =1n s =3n s =5n s =9 (b)Figure 12. Iteration history of iMS-XFEM procedure with different number ofsmoothings per step. Errors for displacement in x (a) and y (b) direction are shown. .2 Test case 2: heterogeneous reservoir with multiple fractures The second test case is set to model deformation in a heterogeneous reservoirwith more fractures. The size and the heterogeneous properties of this testcase are the same as those in test case 1. Here, more fractures are considered.In addition, compared to test case 1, the east and west boundaries are alsoobserving distributed loads, as shown in figure 13.
Figure 13: Test case 2: Multiple fractures within a heterogeneous reservoirunder tension stress across three boundaries.Simulation results for both fine-scale XFEM and MS-XFEM are shown infigure 14. The black lines on figure 14 (b) and (c) represent the coarse scalemesh. It is clear that the MS-XFEM (without using iterations) results in arelatively accurate representation of the deformation field, compared with thefine-scale XFEM, using 8 × a)(b) || e y || = 1 . × − (c) || e y || = 4 . × − Figure 14. Test Case 2: displacement field for a heterogeneous media with multiplefractures for (a) fine scale XFEM and (b) MS-XFEM without iterative strategyapplied (c) iMS-XFEM result after 3 iterations.
An example of two basis functions for this test case is shown in figure 15. Inthe local plot, it is illustrated how 1 (15a) and 2 (15b) fractures are capturedby the basis functions. (a) (b)Figure 15. Illustration of the basis functions P xx for two different part of the domain,where two (a) and one (b) discontinuities are captured. The iMS-XFEM procedure, as explained before, is now applied to reduce themultiscale errors to machine precision. Still in figure 14 (c) the result qualityis improved a lot after 3 iterations with 3 times fine-scale smoothers appliedin the second stage. Note that since no GMRES nor a complete smoother isused, but the incomplete smoother ILU(0) for its efficiency. Results are shownin figure 16. 20
20 40 60 80 100 120 iteration -11-10-9-8-7-6-5-4-3-2 l og10 ( e x ) n s =1n s =3n s =5n s =9 (a) iteration -11-10-9-8-7-6-5-4-3 l og10 ( e y ) n s =1n s =3n s =5n s =9 (b)Figure 16. Iteration history of iMS-XFEM procedure with different number ofsmoothing per step. Errors for displacement in x (a) and y (b) direction are shown. A multiscale procedure for XFEM is proposed to model deformation of geo-logical heterogeneous fractured fields. The method resolves the discontinuitiesthrough local multiscale basis functions, which are computed using XFEMsubjected to local boundary conditions. The coarse-scale system is obtainedby using the basis functions, algebraically, which does not have any additionalenrichment functions, in contrast to the local basis function systems. Thisprocedure makes the MS-XFEM very efficient. Also, by combining it with afine-scale smoother, an iterative MS-XFEM (iMS-XFEM) procedure is devel-oped, which allows to reduce the error to any desired level of accuracy.Two heterogeneous test cases were studied as proof-of-concept, to investigatethe performance of the MS-XFEM. It was shown that MS-XFEM results inacceptable solutions, when no iterations are imposed. By applying iterations,one can further improve the results. For practical applications, when param-eters are uncertain, only a few iterations can be applied to maintain (andcontrol) the MS-XFEM quality of the solution.
Acknowledgements
Fanxiang Xu is sponsored by the Chinese Scholarship Council (CSC). Authorsacknowledge Yaolu Liu of TU Delft and all members of the Delft AdvancedReservoir Simulation (DARSim) and ADMIRE research group for fruitful dis-cussions. 21 eferences [1] C. A. Barton, M. D. Zoback, and D. Moos, “Fluid flow along potentially activefaults in crystalline rock,”
Geology , vol. 23, pp. 683–686, 08 1995.[2] E. L. Majer, R. Baria, M. Stark, S. Oates, J. Bommer, B. Smith, andH. Asanuma, “Induced seismicity associated with enhanced geothermalsystems,”
Geothermics , vol. 36, no. 3, pp. 185 – 222, 2007.[3] H. Hajibeygi and P. Jenny, “Multiscale finite-volume method for parabolicproblems arising from compressible multiphase flow in porous media,”
Journalof Computational Physics , vol. 228, pp. 5129–5147, aug 2009.[4] N. Castelletto, H. Hajibeygi, and H. A. Tchelepi, “Multiscale finite-elementmethod for linear elastic geomechanics,”
Journal of Computational Physics ,vol. 331, pp. 337–356, feb 2017.[5] I. Sokolova, M. G. Bastisya, and H. Hajibeygi, “Multiscale finite volume methodfor finite-volume-based simulation of poroelasticity,”
Journal of ComputationalPhysics , vol. 379, pp. 309–324, feb 2019.[6] P. Jenny, S. Lee, and H. Tchelepi, “Multi-scale finite-volume method for ellipticproblems in subsurface flow simulation,”
Journal of Computational Physics ,vol. 187, pp. 47–67, may 2003.[7] R. Deb and P. Jenny, “Finite volume-based modeling of flow-induced shearfailure along fracture manifolds,”
International Journal for Numerical andAnalytical Methods in Geomechanics , vol. 41, pp. 1922–1942, jul 2017.[8] G. Ren, J. Jiang, and R. M. Younis, “A fully coupled XFEM-EDFM model formultiphase flow and geomechanics in fractured tight gas reservoirs,”
ProcediaComputer Science , vol. 80, pp. 1404–1415, 2016.[9] N. Castelletto, H. Hajibeygi, and H. Tchelepi, “Hybrid multiscale formulationfor coupled flow and geomechanics,” in
ECMOR XV - 15th European Conferenceon the Mathematics of Oil Recovery , EAGE Publications BV, aug 2016.[10] A. Fumagalli and A. Scotti, “An efficient XFEM approximation of darcy flowsin arbitrarily fractured porous media,”
Oil & Gas Science and Technology –Revue d’IFP Energies nouvelles , vol. 69, no. 4, pp. 555–564, 2014.[11] B. Giovanardi, L. Formaggia, A. Scotti, and P. Zunino, “Unfitted FEM formodelling the interaction of multiple fractures in a poroelastic medium,” in
Lecture Notes in Computational Science and Engineering , pp. 331–352, SpringerInternational Publishing, 2017.[12] H. Hajibeygi, G. Bonfigli, M. A. Hesse, and P. Jenny, “Iterative multiscale finite-volume method,”
Journal of Computational Physics , vol. 227, pp. 8604–8621,oct 2008.
13] Y. Wang, H. Hajibeygi, and H. A. Tchelepi, “Algebraic multiscale solver for flowin heterogeneous porous media,”
Journal of Computational Physics , vol. 259,pp. 284 – 303, 2014.[14] E. T. Chung, Y. Efendiev, and G. Li, “An adaptive gmsfem for high-contrastflow problems,”
Journal of Computational Physics , vol. 273, pp. 54 – 76, 2014.[15] M. Rashid, “The arbitrary local mesh replacement method: An alternativeto remeshing for crack propagation analysis,”
Computer Methods in AppliedMechanics and Engineering , vol. 154, pp. 133–150, feb 1998.[16] T. Bittencourt, P. Wawrzynek, A. Ingraffea, and J. Sousa, “Quasi-automaticsimulation of crack progation for 2d lefm problems,”
Engineering FractureMechanics , vol. 55, pp. 321–334, 09 1996.[17] R. Cook,
Finite Element Modeling for Stress Analysis . Wiley, 1995.[18] G. Wells and L. Sluys, “Three-dimensional embedded discontinuity modelfor brittle fracture,”
International Journal of Solids and Structures , vol. 38,pp. 897–913, feb 2001.[19] M. T¸ ene, S. B. Bosma, M. S. A. Kobaisi, and H. Hajibeygi, “Projection-basedembedded discrete fracture model (pEDFM),”
Advances in Water Resources ,vol. 105, pp. 205–216, jul 2017.[20] A. R. Khoei, M. Vahab, E. Haghighat, and S. Moallemi, “A mesh-independentfinite element formulation for modeling crack growth in saturated porousmedia based on an enriched-FEM technique,”
International Journal of Fracture ,vol. 188, pp. 79–108, jun 2014.[21] Y. Efendiev, J. Galvis, G. Li, and M. Presho, “GENERALIZED MULTISCALEFINITE ELEMENT METHODS: OVERSAMPLING STRATEGIES,”
International Journal for Multiscale Computational Engineering , vol. 12, no. 6,pp. 465–484, 2014.[22] J.-Y. Wu and F.-B. Li, “An improved stable XFEM (is-XFEM) with anovel enrichment function for the computational modeling of cohesive cracks,”
Computer Methods in Applied Mechanics and Engineering , vol. 295, pp. 77–107,oct 2015.[23] T. Belytschko, Y. Y. Lu, and L. Gu, “Element-free galerkin methods,”
International Journal for Numerical Methods in Engineering , vol. 37, pp. 229–256, jan 1994.[24] J. Melenk and I. Babuˇska, “The partition of unity finite element method:Basic theory and applications,”
Computer Methods in Applied Mechanics andEngineering , vol. 139, no. 1, pp. 289 – 314, 1996.[25] N. Moes, J. Dolbow, and T. Belytschko, “A finite element method for crackgrowth without remeshing,”
International Journal for Numerical Methods inEngineering , vol. 46, pp. 131–150, sep 1999.
26] T. Belytschko and T. Black, “Elastic crack growth in finite elementswith minimal remeshing,”
International Journal for Numerical Methods inEngineering , vol. 45, no. 5, pp. 601–620, 1999.[27] A. M. Arag´on and A. Simone, “The discontinuity-enriched finite elementmethod,”
International Journal for Numerical Methods in Engineering , vol. 112,pp. 1589–1613, sep 2017.[28] F. P. Meer and L. J. Sluys, “A phantom node formulation with mixed modecohesive law for splitting in laminates,”
International Journal of Fracture ,vol. 158, pp. 107–124, may 2009.[29] G. N. Wells and L. J. Sluys, “A new method for modelling cohesive cracks usingfinite elements,”
International Journal for Numerical Methods in Engineering ,vol. 50, no. 12, pp. 2667–2682, 2001.[30] M. HosseiniMehr, C. Vuik, and H. Hajibeygi, “Adaptive dynamic multilevelsimulation of fractured geothermal reservoirs,”
Journal of ComputationalPhysics: X , p. 100061, may 2020.[31] M. HosseiniMehr, M. Cusini, C. Vuik, and H. Hajibeygi, “Algebraic dynamicmultilevel method for embedded discrete fracture model (f-ADM),”
Journal ofComputational Physics , vol. 373, pp. 324–345, nov 2018.[32] E. Chow and Y. Saad, “Experimental study of ILU preconditioners for indefinitematrices,”
Journal of Computational and Applied Mathematics , vol. 86, pp. 387–414, dec 1997.[33] H. Zhou and H. A. Tchelepi, “Two-stage algebraic multiscale linear solver forhighly heterogeneous reservoir models,”
SPE Journal , vol. 17, pp. 523–539, jun2012.[34] J. T. Camargo, J. A. White, and R. I. Borja, “A macroelementstabilization for mixed finite element/finite volume discretizations of multiphaseporomechanics,”
Computational Geosciences , 2020.[35] K. M. Terekhov, “Cell-centered finite-volume method for heterogeneousanisotropic poromechanics problem,”
Journal of Computational and AppliedMathematics , vol. 365, p. 112357, 2020.[36] H. Wang,
Theory of Linear Poroelasticity with Applications to Geomechanicsand Hydrogeology . Princeton Series in Geophysics, Princeton University Press,2017.[37] F. Gaspar, F. Lisbona, and P. Vabishchevich, “A finite difference analysis ofbiot’s consolidation model,”
Appl. Numer. Math. , vol. 44, no. 4, pp. 487 – 506,2003.[38] N. Castelletto, S. Klevtsov, H. Hajibeygi, and H. A. Tchelepi, “Multiscaletwo-stage solver for biot’s poroelasticity equations in subsurface media,”
Computational Geosciences , vol. 23, no. 2, pp. 207–224, 2019.
39] Y. Saad and M. H. Schultz, “GMRES: A Generalized Minimal ResidualAlgorithm for Solving Nonsymmetric Linear Systems,”
SIAM J. Sci. Stat.Comput. , vol. 7, no. 3, pp. 856–869, 1986., vol. 7, no. 3, pp. 856–869, 1986.