CC BY
23Математические заметки СВФУ Октябрь—декабрь, 2019. Том 26, № 4
UDC 519.63
SPLITTING SCHEMES FOR THE
THERMOPOROELASTICITY PROBLEM IN FRACTURED MEDIA
D. A. Ammosov, M. V. Vasilyeva, M. Babaei, and E. T. Chung
Abstract. We consider the thermoporoelasticity problem in the fractured geothermal reservoir. We use a hierarchical fracture representation, where small-scale highly connected fractures are represented by the classical dual porosity model and large-scale dense fractures are represented with the use of a discrete fracture model. The mathematical model is described by a coupled system of equations for temperature and pressure in the coupled dual continuum porous media with discrete fractures, where deformations are considered based on the effective media approach. For the numerical solution, we construct unstructured grids that resolve large-scale fractures explicitly on the grid level for the mixed-dimensional formulation of the pressure and temperature equations. The discrete system is constructed based on the finite element method with an implicit scheme for approximation by time. For effective solution of the obtained coupled system of equations for pressures, temperatures, and displacements for multicontinuum media, we present and study the splitting schemes based on fixed stress splitting. The results of the numerical simulation for the two-dimensional problem and a numerical study of the splitting schemes for the model problems are presented for two sets of parameters to show stability of the proposed schemes.
DOI: 10.25587/SVFU.2019.73.18.009 Keywords: thermoporoelasticity, double porosity, double permeability, dual continuum, discrete fracture model, finite element method, mathematical modeling, splitting schemes, fixed stress splitting, mixed-dimensional problem
1. Introduction
Numerical simulation of the thermoporoelasticity problems has many applications, for example, heat recovery from Enhanced Geothermal Systems (EGS), nuclear waste disposal and other coupled subsurface problems [1—3]. The mathematical model is described by a coupled system of equations for fluid flow, heat transfer and geomechanical deformations of the porous media [4-6]. Solving the systems of equations obtained from such multiphysical problems using the coupled method requires large computational resources. Splitting scheme is an alternative method that helps to solve the coupled system of equations sequentially, and thereby reduce
The work of D. A. Ammosov and M. V. Vasilyeva was supported by the mega-grant of the Russian Federation Government (No. 14.Y26.31.0013) and Russian Science Foundation (No. 17— 71-20055).
© 2019 D. A. Ammosov, M. V. Vasilyeva, M. Babaei, and E. T. Chung
computational costs. Splitting schemes for poroelasticity problems were considered in [7-10], the schemes for thermoporoelasticity problems were studied in [11,12]. Sequentional solution of the poroelasticity or thermoporoelasticity problems can lead to unstable solution in some cases of model parameters and, therefore, some stabilization or regularization should be applied. Stable splitting schemes were presented in [13-15].
Subsurface reservoirs usually have highly heteregeneous properties, moreover, they have fracture distribution in different scales. For simulation of the fractured porous media with dense highly connected fracture distribution, dual porosity models are used [16-18]. For simulation of the large-scale hydraulic fractures, we should construct unstructured grid that directly resolves complex fracture geometry. For the explicit consideration of the fracture's geometry, a mixed dimensional problem is used to describe a complex interaction between fracture and matrix [19-21]. For explicit fracture simulations, we have two main types of models (embedded fracture model and discrete fracture model), that differ by a mesh construction for porous matrix and fracture [22-24]. In the embedded fracture model, mesh for the fracture network and porous matrix is constructed separately, but for the discrete fracture model, we construct conforming fracture mesh with porous matrix mesh [25-27].
In the numerical simulation of the geothermal reservoirs, we have different scales for the fracture networks, one is for highly connected natural fracture networks and another for the large-scale hydraulic fractures. For describing the flow and heat transfer in the natural fracture networks, we use the dual continuum model and the discrete fracture for the large-scale hydraulic fractures. Mechanic deformation is described using multicontinuum approach, where effective stress contains parts for each continuum for pressure and temperature [28,29].
The structure of the paper is as follows. In Section 2, we present the mathematical model for the thermoporoelasticity problem in the fractured porous media. We introduce an approximation of the model using a finite element method in Section 3. The considered splitting schemes are presented in Section 4. The results of the numerical simulation for a two-dimensional problem and a numerical study of the splitting schemes for problems of filtration, poroelasticity and thermoporoelasticity are given in Section 5. Finally, we present conclusions in Section 6.
2. Mathematical model
The mathematical model of the thermoporoelasticity problem is described by the coupled system of equations for the displacements u, pressure p, and temperature T. For problems with the small-scale connected network of the fractures, we use the double porosity-double permeability (dual continuum) model for describing heat and mass transfer, where we introduce the two continua in the domain £ C (d = 2, 3) (Fig. 1) for the porous matrix and the small scale connected fracture network (pi, P2, T1, T2). The subscript 1 corresponds to the porous matrix, and the subscript 2 corresponds to the small-scale connected network of the fractures. We suppose that the small scale connected fracture network and the porous matrix continua exist in
each mesh cell. We introduce lower-dimensional domain 7 C Rd-i (Fig. 1), in which we define pressure and temperature of large- scale fractures (pf, Tf) with the heat and mass transfer terms. For the displacements u G £ C Rd, we use Biot's model with effective stress tensor that depends on pressure and temperature of the first continuum, for simplicity.
f2
Fig. 1. The domains £ C Rd and 7 C Rd-i, d = 2. Therefore, we have the following coupled system of equations for triple continuum
= 0, x G = 0, x G
= 0, x G
= 0, x G 7,
(1)
= 0, x G
= 0, x G
= 0, x G 7
j=1,2
model:
-V • a + aiVpi + ^iVTi
dV^ u dpi dTi , , , ,
+01-^-~ V • (fciVpi) + dij[Pi-Pj)
j=2,f
b-2^- - '72^ - V • (k2Vp2) + d2j(P2 - Pj)
j=1,f
C\ f)rT1
'»-¡jf- - >I>-Jf - V • (kfVpf) + Y. dfj(Pf -Pj)
j=1,2
wJ^+C.f+CTVin.Vt-n.T,,*.
-V • (AiVTi)+ rij(Ti - Tj) j=2,f
c2^ + C2№V(T2 • v2) - mT-2,0^
-V • (A2VT2) + £ r2j (T2 - Tj)
j=1,f
Cfd^ + crv(Trvf)-llfTf/-§i -V • (AfVTf) + £ rfj (Tf - Tj)
where a and e are the stress and strain tensors a = A(V • u)I + 2^e,
and A, ^ are the Lamé coefficients, bi = 1/M^ Mi is the Biot modulus, Vi = — kiVpi is the Darcy's velocity, ki = Ki/vf is the mobility of the ¿-th continuum , Ki is the permeability of the ¿-th continuum, Vf is the fluid viscosity, Ai is the thermal conductivity of the ¿-th continuum, ai is the Biot-Willis constant, pi = ^i(A + 2/3^), £i is thermal expansion of the solid phase, Ci is the equivalent volumetric heat capacity, ni is the thermal expansion coefficient, dij and rij are the mass and heat transfer coefficients (i = 1, 2, f ) [30,31]. We note that, in this model we ignore gravity forces and suppose diffusive transfer term for the heat equations.
System of equations (1) is considered with the following initial conditions
u
= Uo, Pi = PiiQ,
Ti = Ti
i,0,
(2)
and the boundary conditions:
a ■ n = 0, x G I\, u = 0, x G T2,
dp
dp2
dpf
-fci—— =—k2-z—= 0, x G diî, -kf—J- = 0, x G IV, pf = g, xGTd
dn
dTi
d n
d n
dTf
-Ai——- = = 0, X G diî, -\f—L = 0} x G IV, Tf = q, XGTD
dn
dn
(3)
where Ti U r2 = Oil, TD U TN = dy.
3. Finite-element approximation
We construct an unstructured grid, where the fractures are located on the interfaces between the cells. We use the continuous Galerkin finite element method with linear basis functions. For finite element approximation, we define spaces V = H(Q), VY = {v G H(7) : v|rD = g}, W = {w G [H(Q)]d : w|r2 = 0} and u G W, pi G V, Ti G V, P2 G V, T2 G V, pf G VY, Tf G VY. We have following variational formulation of the problem: find (un+1, p^+S p£+1, pnf+1, Tf+S T2n+1, Tfn+1) G W x V x V x VY x V x V x VY, such that
J a(un+1) : e(w) dx + J ai Vpi+1 • w dx + J foVTn+1 • w dx = 0,
' V- un+1 — V- un f p?+1 — pn , ct\-v\ dx + / b 1-v\dx
ni
^n+1 _ Tin
v1 dx + J k1 Vpl+1 ■ Vv1 dx
+ J d12 (p?+1 _ pi+1)v1 dx + J dff (pi+1 _ pi}+1)v1 dx = 0
C „"+1 _ pn r Tn+1 - Tn r
/ b/2 v2 d,x- t}2 ~ --2-v2 dx+ k2Vp™+1 • Vv2 dx
n n n
- J di2 (p?+1 - pn+1)v2 dx + J d2f (pn+1 - pn+1)v2 dx = 0, n Y
r pn+1 - p? f f1 - f r +_,
/ bf—-j-vf dx - i]f—-dx+ kfVpy ■ Vvf dx
Y Y Y
-/ dif (p?+1 - p?+1)vf dx -J d2f (p?+1 - p?+1)vf dx = 0, (4)
Y Y
/V u"+i - V un i T"+1 - T" /"
/3iTi,0-zid®+ / Ci—-—zi da:+ / C^V^^-Vi)*! d®
n n n
-J ViTi,oPl+ dx + J AiVT™+1 • Vzi da; + J r12 (T?+1 - T?+1)Zl dx
n n n
+ j rif (T1n+1 - Tfn+1)zi dx = 0,
Y
/Tn+i _ Tn /*
--z2 dx+ C2"'V(T2Tl+1 • V2)z2 dx
nn
-J V2T2,oP2+ ~ z2 dx + J A2VT™+1 • Vz2 da; n n
-J ri2 (Tin+i - T^1) Z2 dx + J r2f (T^1 - T/n+i) Z2 dx = 0,
nY
/Tn+1 _ Tn
Cf^—^zfdx + J Cjy■ Vf)zf dx
YY
-J VfTf,oP" dx + j\fVT™+1-Vzf dx
Y Y
-/rif (Tin+1 - T/n+1)zf dx -J r2f (T2n+1 - T/n+1)zf dx = 0,
YY
where (w, v1, v2, vf, z1, z2, zf) G W x V x V x V0Y x V x V x V0Y. We note that, we used an implicit scheme for approximation by time with time step t .
Let Th be a mesh for the domain <fh be the faces of the grid Th, and <fY will be the faces of large-scale fractures, where S'Y C S'h. We obtain the following discrete system in the matrix form for yn+1 = (un+1,pn+1,pn+1,pn+1,T1n+1,T2n+1,Tfn+1)T
yn+1 - yn
M--— + A yn+1 = 0,
t
where
(0 0
D Ci
M =
0 0 G
0 0
0 0
Ni
0 0
0 0
C2
0 0
N2
0
0 0 0
Cf
0 0
Nf
/ a d 0 0
0 B1 + Q+Z1 -Q -Zi
A:
0
Ni
0 0
Mi
0 0
G
0
0 0
N2
0 0
M2
0
0 0 0
0 0
-Q B2+Q + Z2 -z2
-zi -z2 bs + Z1+Z2
0 0
s1+h+r1
0 0 0
Nf
0 0
Mf
0 0 0 0 -H
0 0 0 0
r1
-H S2 + H+R2 -R2 -R1 -R2 Sf +R1+R2 )
For matrices, we have
A = {flij }
- I aij} 5 aij -
• DD - DT, D - {dij}, dij - / aiV* • dx,
n
• G - Ti,oGT, G - {gij}, gij - / ftV* • dx,
n
• C, - {cl}i3}, l -1, 2,f
Cii - / * dX, ^ -i Mj * dX Cfi - / f *fj f dX,
n n Y
• N - T,oNT, N, - {n,,ij}, l - 1, 2, f
ni'j - -J * dX, - -J V2*j * dX, nfi - -J Vf f f dX,
n n Y
• B, - {hj}, l -1, 2, f
bi,j - / kiV* •V*i dX, ^ -j W* •V*i dX, f -j kf Vf •V*f, dX,
nnY
• Q - {?ij} and H - {hij} qij - / di2*j * dX, hj - / ri2* * dX
n n
• Z, - {z^j}, l - 1, 2 with zi,ij - / dif *f,j*f,i dX, z2,j - / d2f*fj*f,i dX
YY
• Mi - {mi,ij}, l -1, 2, f
mi,j -j Ci** dX, m2,j -j C2** dX, mf,j -j Cf f f dX,
0
• S, = }, I =1, 2, f
sMj = J AiV^j • V^ dx + J CWV(^j • Vi)0i dx, n n
S2,ij = J A2V^j • V^i dx + J CW V(0j • V2)^i dx, n n nn
sf,ij = J Af V^f,j • V^f,idx + f CW V(f • Vf )^f,idx,
YY
• R, = {r,,ij}, I = 1, 2 with ri,ij = / rif-f,jf dx, = / r2f-f,j -f dx
YY
We use linear basis functions for domain £ and lower dimensional basis for domain 7. We use the artificial diffusion for stabilization of the temperature equation, A, = A, + Aad(V,) and take velocity from the previous time step pressure solution.
4. Splitting schemes
We consider two splitting schemes: fixed strain splitting scheme and fixed stress splitting scheme. In [32], using a von Neumann stability analysis, it was shown that fixed strain splitting scheme is conditionally stable and fixed stress splitting scheme is unconditionally stable. Contraction estimates and convergence rates for fixed stress splitting scheme were derived in [33, 34].
4.1. Fixed strain splitting scheme. In fixed strain splitting scheme, the equations for pressures and temeperatures are solved first, the change in volumetric strain, i.e. displacements divergence remains constant and is taken from the previous time step. Then the displacements are calculated using the previously obtained pressure and temperature of the first continuum.
• Fixed strain splitting scheme 1: p^1 ^ Tfl+1 ^ p^1 ^ T2n+1 ^ p^1 ^ T1n+1 ^ un+1.
In this scheme, we solve the equations, starting from the most permeable (heat-conducting) continuum f to the least permeable (heat-conducting) continuum 1. For each continuum, we first find the pressure and then the temperature. The divergence of the displacements is taken from the previous time step. At the end we find the distribution of displacements.
Aun+1 + Dpn+1 + GT1n+1 = 0,
- un-i pn+1 - pn Tn - Tn-i
D--— + cP-^ + JVi 1 1 + BlPl+1 + Q(pi + 1 -Pn2+^
+ Zi( pn+1 - pn+1)
0
pn+1 — pn Tn — Tn-i CV_2-P2 + N2h-±2_ + B2pn+1 + „+1 _p„x + z , „+1 _p„+lx = Qj
T T J
pn+i — „n Tn — Tn-i C/f T f + Nf f Tf + BfPnf+1 + Z^1 -p?) + Z2(p?+1 -p£) = 0,
~ un — ».n-1 Tn+i — Tn __ Pn+i — pn
GU-U- + Mi Ji-£i_ + ^Pi-+ + 1 + H{Tn+l _ T„+n
T T T
+ Ri (T"+i — Tfn+i) -0, (5)
T n+ i — T? ~ pn+1 — pn
f
T n+i — T? — pnt+i — p? Mf-£-i- + at/ —-^ + s'/t^1 + R1 (t™+1 - tf) + R2 - t2n) = 0,
I jyP/ ~P>
t f T
Fixed strain splitting scheme 2: (p^p^p^1) ^ (Tin+i, T2n+i, Tf+i)
un+1.
In this scheme, we solve all the equations for the pressures by the coupled method, and then we use the obtained pressures to solve the equations for the temperatures by the coupled method. The divergence of the displacements is taken from the previous time step. At the end we find the distribution of displacements.
Awn+i + Dpn+i + GTin+i - 0,
~ un — o.n-i pn+i — pn Tn — Tn-i
DU-V_ + CiEl-^ + , , + + 1 + +1
T T T
+ Zi(pn+i — pnf+V) -0,
pn+1 — pn Tn — Tn-i
C2P2-P2 + N/2 2 + B2pn+1 +Q{ n+1 _pn+l) + ^ { n+l „+!) = ^
T T J
pn+i — pn t? — T?-1
Cf f f +Nf f Tf + Bfpn+1 + Z1(pn+1 -p™+1)+Z2(pnf+1 -p™+1) = 0,
~ un — un-1 Tn+1 — Tn — nn+1 — pn
Q U-U- + MiJi-+ ^p,-^ + +1 + H{Tn+l _ T„+n
T T T (i i ) 2
+ ri (t?+1 — t;+1) -0, (6)
Tn+1 — T? ~ pn+1 — pn
T T 2 2 i 2 f
T n+1 — T? — «n+1 — pn Mf-L-—+Nf—-^+S,/T;+1+R1(T/ri+1-T1n+1)+R2(T/ri+1-T2ri+1) = 0,
4.2. Fixed stress splitting scheme. For simplicity, we suppose that the change in volumetric strain in the heat transfer equations remains constant, and we
take the displacements divergence from the previous time step. Then we assume that during filtering volumetric mean total stress is constant, and obtain:
2
Dun+1 = Dun + ^ (p™+1 - p™), (7)
After substitution (7) into the equation for the first continuum pressure in (5) we get:
un_un-l n+ l_pn 2pn+l_pn a2 pn_pn-1
L>--h Ci--h { —---—-)
t t K t K t
Tn _Tn— 1
+ JV! 1 T 1 + + + Z^1 - pf1) = 0, (8)
We can also write (8) in the form of a regularized scheme:
D----+ C^-^ + ST^-+Pl
t t t 2
Tn _ tn—1
+ Aii 1 t 1 + -BlPi+1 + QW+1 - Pn2+1) + Z1(pn1+1 - pnf+1) = 0, (9)
where iff, = j^, S = af, K = A + |/i is the drained bulk modulus.
For obtaining corresponding splitting schemes 1 and 2 we just replace the equation of the first continuum pressure in (5) and (6) by (9).
5. Numerical results
In this section, we present the numerical results for the thermoporoelasticity problem defined using the presented finite element approximation and study of the proposed splitting schemes. We consider two-dimensional problem in domain Q = [0,10] x [0,10] m2. Numerical implementation is based on the open source library FEniCS [35]. In order to resolve the fracture networks by the mesh, we construct unstructured grids using Gmsh software.
The computational grid for the domain Q has 15520 vertices and 30638 cells, the grid for y has 645 vertices and 640 cells. We depict the grids in Fig. 2, where we represent the unstructured grid with triangular cells for the two-dimensional domain Q on the left picture, the fracture domain y is represented on the left and right pictures.
Two sets of parameters are presented in Tab. 1 that have different elastic parameters y«, A, a1. We suppose that the first continuum represents the porous matrix, the second continuum represents the dense small-scale highly connected fracture networks, and the third continuum represents the large-scale fracture network. For transfer terms, we set <r12 = 2 • 10—13 and r12 = 10 for interaction between the first and the second continua, = 2 • 10—12 and r1f = 100 for interaction between the first and the third continua. Also we set <r2j = 10—11 and r2f = 1000 for interaction
Table 1. Problem properties
Parameter Unit Set 1 Set 2
Mi Pa 107 107
M2 Pa 107 107
Mf Pa 107 107
Kl 2 m 2 io-16 2 io-16
K2 2 io-15 io-15
Kf 2 io-12 io-12
V Pa-s 10"3 10"3
VI oc-l 3 • 10-4 3 • IO-4
m oc-l 3 • 10-4 3 • 10-4
rif oc-l 3 • 10-4 3 • IO-4
Ci J • m_3 _ oc-l 107 107
c2 J • m_3 _ oc-l 107 107
Cf J • m_3 _ oc-l 107 107
Ai W m-1• 10 10
a2 w m-1• 100 100
w m-1• 105 105
«1 - 1 0.1
a oC-l 3 • IO"5 3 • IO"5
A Pa 106 109
Pa 106 109
between the second and the third continua. We simulate for tmax = 86400 [s] with 100 time steps. We have following initial conditions: To = 200 [°C], po = 107 [Pa], «0 = 0 [m]. We set boundary conditions as follow: q = 20 [°C] and g = 1.2 • 107 [Pa] on the left boundary for the third continuum. On other boundaries we set zero Neumann boundary conditions, u • n = 0 [m] on the left and the bottom boundaries, on other boundaries we set a • n = 0 [Pa].
The results of the numerical simulations are presented in Fig. 3-5. In Fig. 3 we depict pressures for the first, the second and the third continua for different time steps, t^o, t^o and tioo. Temperature distributions for the first, the second and the third continua at times of t^o, t6o and t1oo are presented in Fig. 4. The displacements distribution at the final time is presented in Fig. 5. We observe how pressure increase from left to right, when we inject cold water. Pressure differences lead to the temperature propagation from the left boundary. In Fig. 5, we observe deformations that are caused due to pressure and temperature gradients (the results have been warped by the displacements vector multiplied by 120).
Fig. 2. Computational grids. Left: two-dimensional domain Q. Right: fracture networks y .
1,0e+07 1,2e+07 I.Oe+07 1.2e+07 l,0e+07 l,2e+07
Fig. 3. Pressure distribution (Pa) for different times tn, n = 30, 60 and 100 (from left to right). The first row: the first continuum. The second row: the second continuum. The third row: the third continuum.
The relative errors were calculated by the following formulas:
£Ti = Mk!.10o%, e„ = ■ 100%
\\Pc,i\\L2 \\Tc,i\\L2 \\uc\\L2
where sub index c is related to the coupled solution, i - 1, 2, f, \\ . \\L2 is the L2 norm.
In Fig. 6, 7, we present the dynamics of the relative errors epi, i - 1, 2, f for the filtration problem. Here, we consider three splitting schemes: splitting scheme 1, splitting scheme 2 and splitting scheme 3. In splitting scheme 1, we solve the equations starting from the most permeable continuum f to the least permeable continuum 1 (p?+1 ^ P?+1 ^ Pi + ^. In splitting scheme 2, we solve the equations independently, taking pressures of other continua from the previous time step (p?+1 — P?+1 — Pf+^. In splitting scheme 3, we solve the equations starting from the least permeable continuum 1 to the most permeable continuum f (p?+1 ^ p?+1 ^ Pf+^. In this problem, we have one set of parameters obtained by removing the elastic and heat transfer parameters. As can be seen from the graphs, all splitting schemes are stable for both sets. The values of the relative errors at the final time are presented in Table 2. All splitting schemes have small errors at the final time, but splitting scheme 1 shows the best result.
The dynamics of the relative errors ePi, i - 1, 2, f and eu for the poroelasticity problem are presented in Figures 8—11. Here, we consider two splitting schemes: fixed strain splitting scheme and fixed stress splitting scheme. In fixed strain splitting scheme, we solve the equations starting from the most permeable continuum f to the least permeable continuum 1. The divergence of displacements is taken from the previous time step. At the end we find the distribution of displacements (p>n+1 ^ p?+1 ^ P?+1 ^ . In the second scheme we use fixed stress splitting
scheme approximation for the first continuum pressure equation. In this problem, we consider two sets of parameters obtained by removing the heat transfer parameters from sets 1 and 2. As can be seen from the graphs, fixed strain splitting scheme is unstable for set 1. Fixed stress splitting scheme is stable for both sets. The values of the relative errors at the final time are presented in Table 3.
The dynamics of the relative errors epi, eTi, i - 1, 2, f and eu for thermoporoelasticity problem are shown in Fig. 12-18. As can be seen from the graphs, fixed strain splitting schemes are unstable for set 1. Fixed stress splitting schemes are stable for both sets. The values of the relative errors at the final time are presented in Tables 4, 5. We can see that fixed stress splitting scheme 2 is more accurate than fixed stress splitting scheme 1.
54,11 100 150 200,00 54,11 100 150 200,00 54,11 100 150 200,00
Fig. 4. Temperature distribution (°C) for different times tn, n = 30, 60 and 100 (from left to right). The first row: the first continuum. The second row: the second continuum. The third row: the third continuum.
-0,0036 -0,002 -0,001 0,0000 -0,0043 -0,003 -0,002 -0,001 0,0002
Fig. 5. Displacements distribution (m) at the final time, ux and uy (from left to right).
Fig. 6. Comparison of the relative errors of pressure for the filtration problem. Left: the first continuum. Right: the second continuum.
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 t
Fig. 7. Comparison of the relative errors of the third continuum pressure for the filtration problem.
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
Fig. 8. Comparison of the relative errors of the first continuum pressure for the poroelasticity problem. Left: set 1. Right: set 2.
Fig. 9. Comparison of the relative errors of the second continuum pressure for the poroelasticity problem. Left: set 1. Right: set 2.
Fig. 10. Comparison of the relative errors of the third continuum pressure for the poroelasticity problem. Left: set 1. Right: set 2.
Fig. 11. Comparison of the relative errors of displacements for the poroelasticity problem. Left: set 1. Right: set 2.
Fig. 12. Comparison of the relative errors of the first continuum pressure for the thermoporoelasticity problem. Left: set 1. Right: set 2.
0 10000 20000 30000 40000 50000 60000 70000 80000 30000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
t t
Fig. 13. Comparison of the relative errors of the second continuum pressure for the thermoporoelasticity problem. Left: set 1. Right: set 2.
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
t t
Fig. 14. Comparison of the relative errors of the third continuum pressure for the thermoporoelasticity problem. Left: set 1. Right: set 2.
Fig. 15. Comparison of the relative errors of the first continuum temperature for the thermoporoelasticity problem. Left: set 1. Right: set 2.
F.stressspl. sch. 1 (f->2-sl) ——
Û 10000 20000 30000 40000 50000 60000 70000 80000 00000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
Fig. 16. Comparison of the relative errors of the second continuum temperature for the thermoporoelasticity problem. Left: set 1. Right: set 2.
F. stress spl. sch. 1 !i >2->\) —;— F. strain spl. sch. 1 (t->2->l) —•— F, stress spl, sch, 2 (p-=T-=u) —
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
Fig. 17. Comparison of the relative errors of the third continuum temperature for the thermoporoelasticity problem. Left: set 1. Right: set 2.
0 10000 20000 30000 40000 50000 60000 70000 80000 90000 0 10000 20000 30000 40000 50000 60000 70000 80000 90000
t I
Fig. 18. Comparison of the relative errors of displacements for the thermoporoelasticity problem. Left: set 1. Right: set 2.
6. Conclusions
In this work, we considered thermoporoelasticity problem in the fractured geo-thermal reservoir. The discrete system was constructed based on finite element method with implicit scheme for approximation by time. For effective solution of the obtained coupled system of equations for pressures, temperatures and displacements for multicontinuum media, we presented and studied the splitting schemes that based on fixed stress splitting. The results of the numerical simulation for the two-dimensional problem and the numerical study of the splitting schemes for the model problems were presented for two sets of parameters to show stability of the proposed schemes.
The results of the numerical study for the filtration problem show that splitting scheme 1 has the smallest relative errors in comparison with other schemes. That is why we recommend solving equations starting from the most permeable continuum to the least permeable continuum. The numerical study of the splitting schemes for the poroelasticity problem shows that fixed strain splitting scheme is unstable for set 1, fixed stress splitting scheme is stable for both sets of parameters. The results of the numerical study for the thermoporoelasticity problem show that fixed strain splitting schemes are unstable for set 1, fixed stress splitting schemes are stable for both sets of parameters. We detected that fixed stress splitting scheme 2 is more accurate than fixed stress splitting scheme 1. Therefore, we recommend the fixed stress splitting schemes for the thermoporoelasticity problem in the fractured geothermal reservoir.
Table 2. The filtration problem relative errors at the final time
eP,2 (%) eP,f (%)
Splitting scheme 1 0.053 0.073 0.072
Splitting scheme 2 0.155 0.185 0.145
Splitting scheme 3 0.103 0.112 0.072
Table 3. The poroelasticity problem relative errors at the final time
Set 1 Set 2
(%) Sp,l (%) Sp,2 (%) sP,f (%) (%) Sp,l (%) Sp,2 (%) tp,f (%)
F. stress F. strain 1.204 1.4 • 1050 0.033 2.8 • 1048 0.05 1.2 • 1045 0.051 4.4 • 1045 0.82 0.82 0.052 0.052 0.071 0.071 0.072 0.072
Table 4. The thermoporoelasticity problem relative errors at the final time for set 1
(%) eP,i (%) eP,2 (%) tp,f (%) eT,i (%) ST,2 (%) eT,f (%)
F. stress spl. sch. 1 1.77 0.04 0.062 0.057 0.207 0.519 0.916
F. strain spl. sch. 1 IX I I I I oo
F. stress spl. sch. 2 1.332 0.029 0.012 0.007 0.005 0.005 0.005
F. strain spl. sch. 2 1.5 • 1057 7.1 • 1055 1.3 • 1053 9.5 • 1053 7.9 • 1012 3.8 • 1012 1.4 • 106
Table 5. The thermoporoelasticity problem relative errors at the final time for set 2
(%) Sp,l (%) eP,2 (%) tp,f (%) ST,1 (%) £T,2 (%) eT,f (%)
F. stress spl. sch. 1 0.821 0.052 0.079 0.073 0.211 0.52 0.916
F. strain spl. sch. 1 0.821 0.052 0.079 0.073 0.211 0.52 0.916
F. stress spl. sch. 2 0.003 0.009 0.015 0.007 0.003 0.005 0.006
F. strain spl. sch. 2 0.003 0.009 0.015 0.007 0.003 0.005 0.006
REFERENCES
1. Ghassemi A. and Zhou X., "A three-dimensional thermo- poroelastic model for fracture response to injection/extraction in enhanced geothermal systems," Geothermics, 40, No. 1, 39-49 (2011).
2. Smith D. W. and Booker J. R., "Boundary element analysis of linear thermoelastic consolidation," Int. J. Numer. Anal. Meth. Geomech., 20, No. 7, 457-488 (1996).
3. McTigue D. F., "Flow to a heated borehole in porous, thermoelastic rock: Analysis," Water Resour. Res., 26, No. 8, 1763-1774 (1990).
4. McTigue D. F., "Thermoelastic response of fluid-saturated porous rock," J. Geophys. Res., 91, No. B9, 9533-9542 (1986).
5. Kurashige M., "A thermoelastic theory of fluid-filled porous materials," Int. J. Solids Struct., 25, No. 9, 1039-1052 (1989).
6. Bear J. and Corapcioglu M. Y., "A mathematical model for consolidation in a thermoelastic aquifer due to hot water injection or pumping," Water Resour. Res., 17, No. 3, 723-736 (1981).
7. Girault V., Kumar K., and Wheeler M. F., "Convergence of iterative coupling of geomechanics with flow in a fractured poroelastic medium," Comput. Geosci., 20, No. 5, 997—1011 (2016).
8. Jha B. and Juanes R., "A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics," Acta Geotechnica, 2, No. 3, 139—153 (2007).
9. Kolesov A. E., Vabishchevich P. N., and Vasilyeva M. V., "Splitting schemes for poroelasticity and thermoelasticity problems," Computers Math. Appl., 67, No. 12, 2185-2198 (2014).
10. Wheeler M. F. and Gai X., "Iteratively coupled mixed and Galerkin finite element methods for poro-elasticity," Numer. Methods Partial Differ. Equ., 23, No. 4, 785-797 (2007).
11. Kim J., Unconditionally Stable Sequential Schemes for Thermoporomechanics: Undrained-Adiabatic and Extended Fixed-Stress Splits, Soc. Petrol. Eng. (2015).
12. Brun M. K., Ahmed E., Berre I., Nordbotten J. M., and Radu F. A., "Monolithic and splitting based solution schemes for fully coupled quasi-static thermo-poroelasticity with nonlinear convective transport," arXiv preprint arXiv:1902.05783 (2019).
13. Kim J., Tchelepi H. A., and Juanes R., "Stability and convergence of sequential methods for coupled flow and geomechanics: Fixed-stress and fixed-strain splits," Comput. Methods Appl. Mech. Eng., 200, No. 13, 1591-1606 (2011).
14. Vabishchevich P. N., Vasilyeva M. V., and Kolesov A. E., "Splitting scheme for poroelasticity and thermoelasticity problems," Comput. Math. Math. Phys., 54, No. 8, 1305-1315 (2014).
15. Kolesov A. and Vabishchevich P., "Splitting schemes with respect to physical processes for double-porosity poroelasticity problems," Russ. J. Numer. Anal. Math. Model., 32, No. 2, 99-113 (2017).
16. Barenblatt G. I., Zheltov Iu. P., and Kochina I. N., "Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]," J. Appl. Math. Mech., 24, No. 5, 1286-1303 (1960).
17. Arbogast T., Douglas J. Jr., and Hornung U., "Derivation of the double porosity model of single phase flow via homogenization theory," SIAM J. Math. Anal., 21, No. 4, 823-836 (1990).
18. Warren J. E. and Root P. J., "The behavior of naturally fractured reservoirs," Soc. Petrol. Eng. J., 3, No. 3, 245-255 (1963).
19. Martin V., Jaffre J., and Roberts J. E., "Modeling fractures and barriers as interfaces for flow in porous media," SIAM J. Sci. Comput., 26, No. 5, 1667-1691 (2005).
20. Spiridonov D. and Vasilyeva M., "Simulation of filtration problems in fractured porous media with mixed finite element method (Embedded Fracture Model)," Mat. Zamet. SVFU, 24, No. 3, 100-110 (2017).
21. D'Angelo C. and Scotti A., "A mixed finite element method for Darcy flow in fractured porous media with non-matching grids," ESAIM: Math. Model. Numer. Anal., 46, No. 2, 465-489 (2012).
22. Vasilyeva M., Chung E. T., Leung W. T., and Alekseev V., "Nonlocal multicontinuum (NLMC) upscaling of mixed dimensional coupled flow problem for embedded and discrete fracture models," arXiv preprint arXiv:1805.09407 (2018).
23. Vasilyeva M., Chung E. T., Cheung S. W., Wang Y., and Prokopiev G., "Nonlocal multiconti-nua upscaling for multicontinua flow problems in fractured porous media," arXiv preprint arXiv: 1807.05656 (2018).
24. Vasilyeva M., Babaei M., Chung E. T., and Spiridonov D., "Multiscale modelling of heat and mass transfer in fractured media for enhanced geothermal systems applications," Appl. Math. Model., 67, 159-178 (2018).
25. Karimi-Fard M., Durlofsky L. J., and Aziz K., "An efficient discrete fracture model applicable for general purpose reservoir simulators," in: SPE Reservoir Simulation Symp., Soc. Petrol. Eng. (2003).
26. Gong B., Karimi-Fard M., and Durlofsky L. J., "Upscaling discrete fracture characterizations to dual-porosity, dual-permeability models for efficient simulation of flow with strong gravitational effects," SPE J., 13, No. 1, 58-67 (2008).
27. Chung E. T., Efendiev Y., Leung T., and Vasilyeva M., "Coupling of multiscale and multi-continuum approaches," GEM-Int. J. Geomath., 8, No. 1, 9-41 (2017).
28. Zhang J., Roegiers J.-C., and Bai M., "Dual-porosity elastoplastic analyses of non-isothermal one-dimensional consolidation," Geotech. Geolog. Eng., 22, No. 4, 589-610 (2004).
29. Bai M. and Roegiers J.-C., "Fluid flow and heat flow in deformable fractured porous media," Int. J. Eng. Sci., 32, No. 10, 1615-1633 (1994).
30. Coussy O., Poromechanics, John Wiley & Sons, Chichester (2004).
31. Ammosov D., Vasilyeva M., Babaei M., and Chung E., "A coupled dual continuum and discrete fracture model for subsurface heat recovery with thermoporoelastic effects," Mat. Zamet. SVFU, 26, No. 1, 93-105 (2019).
32. Kim J., Tchelepi H. A., and Juanes R., Stability, Accuracy and Efficiency of Sequential Methods for Coupled Flow and Geomechanics, Soc. Petrol. Eng. (2009).
33. Mikelic A. and Wheeler M. F., "Convergence of iterative coupling for coupled flow and geomechanics," Comput. Geosci., 17, No. 3, 455-461 (2013).
34. Mikelic A., Wang B., and Wheeler M. F., "Numerical convergence study of iterative coupling for coupled flow and geomechanics," Comput. Geosci., 18, No. 3, 325-341 (2014).
35. Logg A., Mardal K.-A., and Wells G. (ed.)., Automated solution of differential equations by the finite element method: The FEniCS book, Springer Science & Business Media, New York (2012).
Submitted August 8, 2019 Revised November 22, 2019 Accepted Novemver 27, 2019
Dmitry A. Ammosov Multiscale Model Reduction Laboratory, Ammosov North-Eastern Federal University, 42 Kulakovsky Street, Yakutsk 677000, Russia dmitryammosov@gmail.com
Maria V. Vasilyeva
Department of Computational Technologies, Ammosov North-Eastern Federal University, 42 Kulakovsky Street, Yakutsk 677000, Russia; Institute for Scientific Computation, Texas A&M University, College Station, TX 77843-3368 vasilyevadotmdotv@gmail.com
Masoud Babaei
The University of Manchester,
School of Chemical Engineering and Analytical Science, Manchester, M13 9PL, UK masoud.babaei@manchester. ac.uk
Eric T. Chung
Department of Mathematics,
The Chinese University of Hong Kong (CUHK),
Hong Kong SAR
tschung@math.cuhk.edu.hk
