# An optimization approach for flow simulations in poro-fractured media with complex geometries

A new discretization approach is presented for the simulation of flow in complex poro-fractured media described by means of the Discrete Fracture and Matrix Model. The method is based on the numerical optimization of a properly defined cost-functional and allows to solve the problem without any constraint on mesh generation, thus overcoming one of the main complexities related to efficient and effective simulations in realistic DFMs.

## Authors

• 11 publications
• 4 publications
• 7 publications
03/24/2021

### A Finite-Volume Moving-Mesh Method for Two-phase Flow in Fracturing Porous Media

Flow in fractured porous media is modeled frequently by discrete fractur...
03/15/2021

### A novel approach for the efficient modeling of material dissolution in electrochemical machining

This work presents a novel approach to efficiently model anodic dissolut...
05/03/2012

### Discretization of a matrix in the problem of quadratic functional binary minimization

The capability of discretization of matrix elements in the problem of qu...
07/23/2018

### Conformal Mesh Parameterization Using Discrete Calabi Flow

In this paper, we introduce discrete Calabi flow to the graphics researc...
12/20/2021

### Flow and Transport in Three-Dimensional Discrete Fracture Matrix Models using Mimetic Finite Difference on a Conforming Multi-Dimensional Mesh

We present a comprehensive workflow to simulate single-phase flow and tr...
02/16/2017

### High Accuracy Mantle Convection Simulation through Modern Numerical Methods. II: Realistic Models and Problems

Computations have helped elucidate the dynamics of Earth's mantle for se...
12/20/2019

### A Three-Field based optimization formulation for flow simulations in networks of fractures on non conforming meshes

A new numerical scheme is proposed for flow computation in complex discr...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The present work deals with the simulation of the flow in the subsoil, modelled by means of the Discrete Fracture and Matrix (DFM) model. According to this model, underground fractures are represented as planar polygons arbitrarily oriented in a three dimensional porous matrix. The flows considered here are governed by the Darcy law in the three dimensional matrix and by an averaged Darcy law on each fracture plane, with suitable matching conditions at fracture-matrix interfaces and at fracture intersections. The quantity of interest is the hydraulic head, given by the sum of the pressure head and elevation. Single phase stationary flow is considered, with the assumption of continuity of the hydraulic head at both fracture-matrix interfaces and at fracture-fracture intersections and no longitudinal flow is allowed along fracture intersections. This is a simplified model with respect to other DFM approaches, described, for example in MJR2005 or, more recently, in Boon2018 , but still representative of realistic configurations, characterized, e.g., by highly permeable fractures. The main focus of the present work is on geometrical complexity aspects, proposing a problem formulation and a numerical approach suitable for complex and randomly generated networks. The described approach can however be extended to different flow models and different coupling conditions. The geometrical complexity of DFM models mainly arises from the multi-scale nature of the resulting domains and from the presence of multiple intersecting interfaces, where the solution displays an irregular behavior. DFM models are proposed as an alternative to homogenization techniques Qi2005 , dual and multy-porosity models DPbook , or embedded discrete fracture matrix (EDFM) models Li2008 ; Moinfar2014 ; 2015WR017729 , and are characterized by the explicit representation of the underground fractures, dimensionally reduced to planar interfaces into the porous matrix. As a consequence of the random orientation, fractures usually form an intricate system of intersections, with the presence of fractures with very different sizes spanning several orders of magnitude that generate intersections with huge geometrical complexities as, for example, 2D and 3D geometrical objects with very different dimension and objects with enormous aspect ratios. The research on effective numerical tools for DFM simulations is particularly active, see e.g. Angot2009 ; Ahmed2015 ; Brenner2016 ; Antonietti2016 ; BPSf ; Odsaeter2019 ; Antonietti2019 ; Chernyshenko2020 . One of the key aspects is the meshing of the domain, with a mesh conforming to the interfaces, suitable for standard approaches for the imposition of interface conditions. The generation of a conforming mesh for realistic fracture networks might, in fact, result in an impossible task, for the extremely high number of geometrical constraints. The mesh conformity constraint at the interfaces can be relaxed by using extended finite elements as suggested, e.g., by FS13 ; formaggia2014 . Different approaches are based on the Mimetic Finite Difference method (MFD) Lipnikov2014 , as described, for example, in Wheeler2015 ; Antonietti2016 , or on Hybrid High Order (HHO) methods as proposed by Chave2018 , where a partial non-conformity is allowed between the mesh of the porous medium and of the fractures, or also on Discontinuous Galerkin discretizations, as in Antonietti2019 . Two or multi-point flux approximation based techniques are described in Sandve2012 ; Faille2016 and gradient schemes in Brenner2016 . Virtual Element (VEM) based discretizations have also been recently investigated to ease the mesh generation process in complex DFMs, as in BBFPSV where the VEM is coupled to the Boundary Element method, and in Fumagalli2019 , in Coulet2019 for poro-elasticity problems, or in Benedetto2020 where an arbitrary order mixed VEM formulation is proposed.

This work presents a development of an optimization-based approach, first proposed for Discrete Fracture Networks BPSa ; BPSb ; BBoS ; BSV ; BBV2019 and recently extended to DFM problems in BPSf . This approach avoids any mesh conformity requirement for the imposition of interface conditions, which are instead enforced through the minimization of a properly defined cost functional. The computation of the quantities involved in functional definition does not require any constraint on the mesh. Further, the resolution of the optimization problem via a gradient-based scheme allows to de-couple the problems on each fracture and the problem on the porous matrix, thus paving the way for an efficient parallel implementation of the numerical scheme, similarly to what done in BSV ; BDV . The discretization scheme described in BPSf relies on the Boundary Element Method for the discretization of the problem on three dimensional matrix blocks, thus requiring the splitting of the original three dimensional domain into sub-domains not crossing the fractures, and thus implying a partial mesh conformity at the fracture-matrix interfaces. Here, the three dimensional domain is not split into sub-domains and Finite Elements are used for the discretization of the matrix, on tetrahedral elements that can arbitrarily cross the fractures. Finite elements on triangular meshes are used for the fractures, with elements not conforming to the tetrahedral mesh and also arbitrarily placed with respect to fracture-fracture intersections. The proposed discretization approach thus greatly improves the usability of the method to general DFM geometries, allowing a trivial meshing process of extremely complex domains, thanks to the complete independence of the mesh from all the interfaces.

The structure of the manuscript is the following: Section 2 describes both the classical and the optimization based formulation of the flow problem in a DFM; the following Section 3 describes the derivation of the discrete problem and the proof of its well posedness; Section 4 shows how an equivalent unconstrained optimization problem is derived, and the gradient based scheme used for problem resolution; Section 5 reports some numerical results and finally some conclusions are proposed in Section 6.

## 2 Problem description

This section is devoted to a brief description of the problem of interest, referring to BPSf for a more detailed exposition and well posedness results. Let us consider a polyhedral block of porous material, denoted as , crossed by a fracture network given by the union of planar polygonal fractures , in the three-dimensional space, i.e. . We further denote by the set of all fracture indexes. Fractures might intersect, and fracture intersections, also called traces, are indicated as , . We assume, for simplicity, that each trace is given by the intersection of exactly two fractures, such that an injective map can be defined between a trace index and a couple of fracture indexes, as being . Further, is the set of indexes of all the traces on fracture and the set of indexes of all the traces in the network. Let us introduce the domain , thus given by the original block without the internal fractures. Calling the boundary of , let us denote by the portion of that matches fracture , for , the superscript “” or “” referring to one of the two sides of the boundary “around” the fracture (see Figure 1

); the unit normal vector to

is , always pointing outward from . A jump operator is introduced for any sufficiently regular vector function on , defined as the jump of along the normal direction to the faces :

 [[v⋅n]]Γi:=(v|Γ+i⋅nΓ+i)−(v|Γ−i⋅nΓ+i).

Similarly, for we denote by the fracture without traces, i.e. , and for each trace , , for any sufficiently regular vector function on , the jump of the normal component of across trace on is denoted as:

with the two sides of the portion of the boundary of lying on and the normal unit vector to with a fixed orientation on . These jump operators are easily extended to functions defined on the whole 3D domain and on the whole fractures , , with the superscripts still denoting the two sides of the interface , , or , .

The portion of not matching any fracture is split in a Dirichlet part and a Neumann part , , where, for simplicity of exposition, we assume homogeneous Dirichlet and Neumann boundary conditions are enforced. Similarly, the boundary of each fracture , , is split in a Dirichlet and Neumann part, and , respectively. If fracture lies in the interior of , then we set , and homogeneous Neumann boundary conditions are prescribed on . If , we assume that , whereas, if there is more than one fracture in the network, we allow for . The problem of the equilibrium distribution of the hydraulic head in can be then stated in strong formulation as:

 −∇⋅(KD∇HD) = fin ~D (1) −∇i⋅(Ki∇iHi) = −[[KD∇HD⋅n]]Γiin ~Fi, i=1,…,NF (2) (HD)|Γ±i = Hi, i=1,…,NF (3) Hi = Hj, on Sm, m=1,…,NS, i,j=σ(m) (4) [[Ki∇iHi⋅n]]Sm = −[[Kj∇jHj⋅n]]Sm, m=1,…,NS, i,j=σ(m) (5) HD = 0on ΓD (6) KD∇HD⋅nΓN = 0on ΓN (7) Hi = 0on γiD, i=1,…,NF (8) Ki∇iHi⋅nγiN = 0on γiN, i=1,…,NF (9)

where is the hydraulic head in , the hydraulic head on , and is a volumetric source term. The operator represents the three-dimensional gradient in , is the two-dimensional gradient on the plane containing fracture , whereas , for is a symmetric positive definite matrix representing the transmissivity of the porous matrix and , is a symmetric positive definite matrix representing the tangential transmissivity of the fracture on its tangential plane. Finally, is the outward unit normal vector to , and for a given index , the outward unit normal vector on the plane of fracture .

Here, for simplicity, we have considered only the source term on the fractures deriving from the exchange with the porous matrix and homogeneous boundary conditions, but the extension to a more general case is immediate. Conditions (3) and (4) express the continuity of the solution at fracture-matrix interfaces and at fracture intersections, respectively, whereas Equation (5) enforces the balance of fluxes at the traces.

Let us now introduce the following functional spaces: first, on each fracture , , we define the function space as ; then on the whole three dimensional domain , the space is defined as the space of functions in whose trace on each interface is a function in , i.e.:

 H1Ω(D)={v∈H10(D):v|ΓD=0, v|Γ±i∈Vi, i=1,…,NF }.

Also, on each trace , we set the spaces and . We introduce the following variables: defined on trace of fracture as

 Umi=[[Ki∇iHi⋅n]]Sm+αHi|Sm,∀i=1,…,NF, ∀m∈Si, (10)

thus representing a sort of internal Robin boundary condition on the traces; and, for all , , with

 Qi:=[[KD∇HD⋅n]]Γi+βHD|Γi, (11)

thus again being a linear combination of the jump of the co-normal derivative of across interface and the trace of on , and the dual of . We remark that, as the hydraulic head is continuous across interfaces .

We also define the bilinear forms: ,

for all , bilinear forms ,

 aFi(vi,wi)=∫FiKi∇ivi∇iwidFi+α∑m∈Si∫Smvi|Smwi|SmdS;

 bi(q,v)=⟨q,v⟩V′i,Vi

and, for , form ,

 cm(um,v)=⟨um,v⟩Um,Hm.

Then, problem (3)-(9) can be written in weak formulation as: find , , , , , such that, for all , for all , :

being the scalar product in . The coupling conditions in weak form are given by: for all , and for all

 bi(HD|Fi−Hi,μi) = 0,∀μi∈V′i, (14) cm(ηm,Hi|Sm−Hj|Sm) = 0,∀ηm∈Um, i,j=σ(m) (15) cm(Umi+Umj,νm) = 0,∀νm∈Hm, i,j=σ(m). (16)

Parameters and ensure stability of the problems written independently on each fracture and on the three dimensional domain. This is required to obtain a discrete formulation suitable for parallel computing. If , or equivalently , well posedness of the problem on the unique fracture is guaranteed by the assumption on .

Problem (12)-(16) is well posed. To show this, let us introduce the function space defined as:

 H1Ω+(D)={v∈H10(D):vi:=v|Fi∈Vi, ∀i=1,…,NF, vi|Sm=vj|Sm, ∀m=1,…,NS,i,j=σ(m)}

and thus incorporating the matching conditions at the interfaces. Let us then write the following problem: find such that, for all

 (17)

Problem (17) is well posed, as it can be easily seen that is an Hilbert space with the scalar product, BPSf :

 (v,w)H1Ω+:=(K∇v,∇w)~D+NF∑i=1(Ki∇iwi,∇ivi)Fi.

Problem (12)-(16) is equivalent to problem (17); indeed, recalling that, for , conditions (14)-(15) are satisfied by construction. Moreover summing (13) for and (12), using (16) and the definition of and , for , , we get (17). We propose a reformulation of problem (12)-(16) well suited for discretization on non conforming meshes and parallel computing, based on a PDE constrained optimization approach. To this end, we introduce a cost functional expressing the error in the fulfilment of the interface conditions as continuity and flux conservation:

 J(HD,HF,US):=NF∑i=1(∥HD|Fi−Hi∥2Vi)+ NS∑m=1(∥Hi|Sm−Hj|Sm∥2Hm+∥Umi+Umj−α(Hi|Sm+Hj|Sm)∥2Um),

being and . The solution to problem (12)-(16) is obtained as the minimum of functional constrained by the PDE equations on the 3D domain and on the fractures:

 minJ(HD,HF,US) (18) constrained by (12) (19) and by (13) ∀i=1,…,NF. (20)

## 3 Discrete formulation

The PDE constrained optimization formulation is specifically designed to allow for an easy discretization of the problem using non conforming meshes and to obtain a discrete problem suitable for effective resolution using parallel computing resources. The imposition of the interface constraints expressed by equations (14)-(16

) with a standard approach requires some sort of mesh conformity at the interfaces: either a perfect matching of the nodes on the meshes to enforce conditions by means of degrees of freedom equality constraints, or the weaker condition of alignment of mesh edges with the interfaces, to use mortaring techniques. In contrast, the imposition of interface conditions through the functional only requires the computation of integrals on the traces, as shown below, and thus meshes can be arbitrarily placed with respect to the interfaces, see Figure

2 for an example of non-conforming meshes in the rock matrix and on the fractures. Further, the minimization process allows to decouple the problems on the fractures and on the three-dimensional domain, for parallel computing.

The discretization strategy proposed in this work is based on the use of standard finite elements on tetrahedra for the three dimensional domain and finite elements on triangles for the fractures. Let us then denote by the tetrahedral mesh on , characterized by a mesh parameter , by a triangular mesh on , , with mesh parameter , and by a possibly different triangular mesh on , with mesh parameter . We further introduce a discretization of the one-dimensional traces, different on each fracture, denoted by , with mesh parameter , , . We denote by the finite dimensional approximation of on , , with the number of degrees of freedom (dofs) and a finite element basis function in 3D; for , we further call the approximation of on , , with the number of dofs and a 2D basis function; the approximation of on , having dofs, and one basis function; the approximation of on , , with dofs and a basis function. Tables 1-2 summarize the labels used for the dimensions of the discrete variables, the name used to denote the basis functions and the notation used, in the following, for the matrices collecting integrals of these basis functions. We build arrays of dofs by collecting column-wise the dofs of each discrete function and with abuse of notation we denote the dof array with the same symbol of the corresponding function, thus having arrays , , , , and , . We define arrays , , as for with , and we further collect column wise arrays , , and forming:

 hF=⎡⎢ ⎢⎣h1⋮hNF⎤⎥ ⎥⎦, q=⎡⎢ ⎢⎣q1⋮qNF⎤⎥ ⎥⎦, ui=⎡⎢ ⎢ ⎢⎣um1i⋮um♯Sii⎤⎥ ⎥ ⎥⎦, u+i=⎡⎢ ⎢⎣um1⋮um♯Si⎤⎥ ⎥⎦, u=⎡⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢⎣u1⋮um⋮uNS⎤⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥⎦,

where are the indexes in ordered increasingly.

The discrete version of functional is the following:

 J(hD,hF,u)=NF∑i=1∥hD|Fi−hi∥2L2(Fi)+ NS∑m=1(∥hi−hj∥2L2(Sm)+∥umi+umj−α(hi+hj)∥2L2(Sm)) (21)

obtained replacing the discretized variables and using norms. The discrete functional can be written in matrix form, computing the integrals of the basis functions and collecting the values into matrices. Considering the first norm in , we have:

 ∥hD|Fi−hi∥2L2(Fi)=∫Fi⎛⎜⎝NhD∑k=1hD,kϕk|Fi−Nhi∑j=1hi,jψi,j⎞⎟⎠2dFi

and we can define three matrices as follows, for each , , , :

 (GiD)k,ℓ=∫Fiϕk|Fiϕk|Fi,(GiDF)k,ℓ=∫Fiϕk|Fiψi,ℓ,(GiF)k,ℓ=∫Fiψi,kψi,ℓ

such that

The computation of matrix is not straightforward, as the two involved variables are defined on different meshes. In particular, the intersection of the three dimensional tetrahedral mesh with the fracture plane needs to be computed. This operation defines a polygonal tessellation of which is then sub-triangulated, thus generating a triangular interface mesh. This sub-triangulation process can be performed without any mesh quality requirement, as the resulting mesh is used only for quadrature purposes. The computation of the elements in is finally performed first computing the intersection of the elements of the interface mesh with the elements in , and subsequently the required integral on the intersection region. Element neighbourhood information is used to efficiently perform the task. The computation of the interface mesh is a quite complex and expensive task. Also in this case element neighbourhood information is used for efficiency, and further can be performed independently fracture by fracture and thus in parallel, which is of paramount importance for the applicability of the method to complex geometries.

We can proceed similarly with the remaining terms of the functional ; to this end, for , and all the possible couples of indexes such that , we define matrices , , :

 (Gmij)k,ℓ=∫Smψi,k|Smψj,ℓ|Sm, (Bmij)k,ℓ=∫Smψi,k|Smϱmj,ℓ, (Cmij)k,ℓ=∫Smϱmi,kϱmj,ℓ

such that

 ∥hi−hj∥2L2(Sm)=[hTihTj][Gmii−Gmij−GmjiGmjj][hihj],

and

 ∥umi+umj−α(hi+hj)∥2L2(Sm)= ([(umi)T(umj)T][CmiiCmijCmjiCmjj]−2α[hTihTj][BmiiBmijBmjiBmjj])[umiumj]

We can collect these matrices, defined locally at the various interfaces into global matrices to derive a compact form of the functional. Let us define matrix , , as a block matrix, with diagonal blocks in positions - are given by , . Extra diagonal blocks in positions - () are instead equal to , if , or a zero block otherwise. Further let us define matrix , and matrices and respectively as

Matrix , , is finally set as

 G=[GFDGFDF(GFDF)TGFF+GS].

For all , let us assemble matrices , , collecting row-wise matrices , for increasing values of and , , i.e.:

 B+i=[Bm1ijBm1ik⋯Bm♯SiipBm♯Siiz]

with , . Let us introduce matrices , with , defined such that . Matrix is finally obtained collecting column wise matrices . Matrix , is a block diagonal matrix with diagonal blocks, each diagonal block in position - being equal to

 [CmiiCmijCmjiCmjj],i,j=σ(m), i

The matrix formulation of then reads as:

We can re-write also the discrete constraint equations in matrix form. We follow a standard procedure and we define matrix