A Decoupled Scheme to Solve the Mass and Momentum Conservation Equations of the Improved Darcy-Brinkman-Forchheimer Framework in Matrix Acidization
** This work is supported by Peacock Plan Foundation of Shenzhen (No. 000255), National Natural Science Foundation of China (No. 11601345) andNatural Science Foundation of SZU (No. 2017059).
Corresponding Author. E-mail address: [email protected].
A Decoupled Scheme to Solve the Mass and Momentum ConservationEquations of the Improved Darcy-Brinkman-Forchheimer Framework inMatrix Acidization*
Yuanqing Wu , Jisheng Kou , Yu-Shu Wu , Shuyu Sun , Yilin Xia College of Mathematics and Statistics, Shenzhen University, Shenzhen, Guangdong 518048, China School of Civil Engineering, Shaoxing University, Shaoxing, Zhejiang 312000, China Department of Petroleum Engineering, Colorado School of Mines, 1600 Arapahoe Street, Goden, CO 80401, USA
Abstract
Matrix acidization simulation is a challenging task in the study of flows in porous media, due to thechanging porosity in the procedure. The improved DBF framework is one model to do this simulation,and its numerical scheme discretises the mass and momentum conservation equations together to form apressure-velocity linear system. However, this linear system can only be solved by direct solvers tosolve for pressure and velocity simultaneously, since zeros appear in the diagonal of the coefficientmatrix. Considering the large-scale attribute of matrix acidization simulation, the solving time of directsolvers is not intolerant. Thus, a decoupled scheme is proposed in this work to decouple the coupledpressure-velocity linear system into two independent linear systems: one is to solve for pressure, and theother one is to solve for velocity. Both of the new linear systems can be solved by parallel and iterativesolvers, which guarantees the large-scale simulation can be finished in a reasonable time period. Anumerical experiment is carried out to demonstrate the correctness of the decoupled scheme and itshigher computing efficiency.
Keywords : matrix acidization, improved Darcy-Brinkman-Forchheimer framework, decoupled scheme,parallel simulation 3D parallel code.
1. Introduction
In petroleum engineering, matrix acidization is a kind of operation, in which the oil productionchannel can be opened up again by the injection of acid flows, after the acid flows dissolve thedepositions such as mud and rocks. However, the numerical simulation of matrix acidization is not aneasy task in the flow study of porous media. Its key characteristic is the change of porosity with theprogress of acidization procedure, which brings a challenge to the simulation task. In order to reflect thechange of porosity accurately, the Darcy-Brinkman-Forchheimer (DBF) framework and its derivationsare provided by Wu et al ([1]-[5]). In this framework, the behaviour of acid flows in the matrix isdecided not only by Darcy’s law, but also by the viscous effect and inertial force, which is representedby the adding of Brinkman term and Forchheimer term to the momentum conservation equation of thetwo-scale model [6], respectively. After this, the DBF framework is named. The meaning of all kinds ofnotations in the following discussions is given in Table 1.Table 1 Nomenclature
Notation Meaning Notation Meaning p pressure t time μ fluid viscosity g gravity vector K permeability value C f cup-mixing concentration of the acid u velocity vector D e effective dispersion tensor ϕ porosity k c local mass-transfer coefficient ρ f mass density of the fluid a v interfacial surface area per unit volume F Forchheimer coefficient C s concentration of the acid at the fluid-solid interface ε pseudo-compressibility factor In the former literature, the mass conservation equation and momentum conservation equation ofthe improved DBF framework [5] (one of the derivations of the DBF framework) are combined togetherto form an equation system to solve for the pressure field and velocity field simultaneously. With thefinite difference method in space and semi-implicit method in time, this equation system can bediscretized to a linear system. Furthermore, such kind of linear system can be solved appropriately withthe serial and direct solvers such as LAPACK [7] and UMFPACK [8], if the linear system is not toolarge. However, as mentioned above, the key characteristic of matrix acidization is the change ofporosity, which requires a fine enough simulation grid to capture the structure of porosity well, i.e., thelinear system will become very large. Unfortunately, the computation of the serial and direct solvers isvery time wasting in such condition, which leads to the intolerable running time of matrix acidizationsimulation. One plausible solution is appealing to the parallel and iterative solvers such as HYPRE [9].However, such kind of solvers meet one main trouble when solving the linear system, which can be seenas below. The momentum conservation equation of the improved DBF framework is ρ f ∂∂t ( uϕ )+ ρ fuϕ ∙∇ uϕ =−∇p− μK u+∇ ∙ μ∇ uϕ − ρ f FK uu+ ρ f g. (1) The mass conservation equation of the improved DBF framework is ∂ϕ∂t +∇ ∙u=0. (2)
After discretization, the linear system
Ax=b is shown in Figure 1. If the 2D domain is considered, thediscretized x -momentum conservation equations, the discretized y -momentum conservation equationsand the discretized mass conservation equations are ranked one by one from the top to the bottom of thelinear system. Accordingly, the x -velocity field U , the y -velocity field V and the pressure field P are putinto the unknown vector x one by one. It is easy to see that the entries of the diagonal at the right-bottomcorner of the coefficient matrix A are zeros, which is the reason that the iterative solvers cannot solvethe linear system. Some works [2] suggest adding a compressibility term to the mass conservationequation, and then Equation (2) is changed to ε ∂p∂t + ∂ϕ∂t +∇ ∙u=0. (3) After such processing, the entries of the diagonal of the coefficient matrix A are non-zeros, and thereforethe linear system can be solved by iterative solvers. However, the introduction of the compressibilityterm changes the attribute of the flows, which should be incompressibility in nature. A trade-off is to setthe pseudo-compressibility factor ε to be a very small positive number.Figure 1 The pressure-velocity linear system after discretization.
2. Decoupled Scheme
Although the trade-off strategy fixes the solvability issue of the linear system with iterative solvers,it cannot manifest the incompressibility attribute of the flows in matrix acidization. Thus, a decoupledscheme is suggested to solve the problem. In this scheme, the Darcy-Forchheimer equation is solvedfirstly, which can be shown as below μK u Dn+1 + ρ f FK u n u Dn+1 =−∇p n+1 + ρ f g, (4) in which u D is an intermediate velocity, and the superscripts represent the time step. By Equation (4), u D can be solved as u Dn+1 = μK + ρ f FK u n −1 −∇p n+1 + ρ f g . (5) In this condition, the mass conservation equation can be written as ∂ϕ∂t +∇ ∙u
Dn+1 =0. (6)
Then, by substituting Equation (5) to (6), the pressure equation can be obtained as ∂ϕ∂t +∇ ∙ μK + ρ f FK u n −1 −∇p n+1 + ρ f g =0. (7) This equation is an elliptic equation with the unknown p n+1 , and therefore all the entries of the diagonalof the coefficient matrix are non-zeros, which guarantees its solvability by iterative solvers. Not onlythat, the incompressibility attribute of the flows is also guaranteed. However, the intermediate velocity u D is computed only from the Darcy term and Forchheimer term, and therefore the Brinkman correctionshould be included as ρ f ∆t u n+1 ϕ n+1 − u n ϕ n + μK + ρ f FK u n u n+1 −u Dn+1 + ρ f ∇ ∙ u n ϕ n ⨂ u n+1 ϕ n+1 − μ∇ n+1 ϕ n+1 =0. (8) Equation (8) is another elliptic equation with the unknown u n+1 . With the same reason as above, it canbe solved by iterative solvers. Moreover, it is easy to see that if Equation (4) is substituted into (8), thediscretized Equation (1) can be achieved, which means that u n+1 in Equation (8) also satisfy (1). Thisis a demonstration that the decoupled scheme has the same velocity solution as the improved DBFframework. In the decoupled scheme, the pressure computed from Equation (7) is an intermediatepressure instead of the true pressure, since it is derived from the Darcy-Forchheimer equation (4). Theintermediate pressure is enough to compute the velocity, and the true pressure will not influence thecomputing of the other variables. Thus, the true pressure will not be computed. In sum, the decoupledscheme has the same solution as the improved DBF framework, except for the pressure.Thus, by the mathematical derivation above, the coupled pressure-velocity linear system of theimproved DBF framework is decoupled into two linear systems: one is to solve for the pressure field,and the other one is to solve for the velocity field. Since both of the two linear systems are derived fromthe elliptic equations, they can be solved by iterative solvers. Furthermore, if parallel and iterativesolvers are used, the large-scale simulation of matrix acidization can be finished in a reasonable timeperiod by a parallel simulator. Meanwhile, since the decoupled scheme has the same solution (exceptpressure) as the improved DBF framework, the incompressibility attribute of flows in matrix acidizationis maintained.Figure 2 The flowchart of the decoupled scheme.
3. Experiment
A numerical experiment is carried out in this section to demonstrate the reliability and applicabilityof the decoupled scheme. For the purpose of comparison with the improved DBF framework, theconcentration equation (9) is introduced to finish the description of the improved DBF framework. ∂(ϕC f )∂t +∇ ∙uC f =∇ ∙ϕD e ∙∇C f −k c a v C f −C s . (9) The other equations of the improved DBF framework will not be given here again, and they can bereferenced in ([1][5]). With the same philosophy as before, one linear system with the unknown C fn+1 isformed, after Equation (9) is discretized. Since Equation (9) is an elliptic equation, the linear system canbe solved by iterative solvers. Until now, there are three linear systems in the decoupled scheme, and allof them can be solved by iterative solvers. The flowchart of the decoupled scheme can be shown inFigure 2.Although the parallel simulator of matrix acidization was realized in the existing works, and theparallel solver such as MUMPS [10] was applied in the simulator, the merits of parallelism are not fullyexploited considering MUMPS is a direct solver. With the decoupled scheme, parallel and iterativesolvers such as HYPRE can be applied in the simulator, and much more acceleration can be gained. Inthe following experiment, three tasks will be done:1. Two single processors with the solver UMFPACK run on the improved DBF framework and thedecoupled scheme, respectively, to investigate the correctness of the decoupled scheme.2. The solver time of both the improved DBF framework and decoupled scheme is given tocompare the computing efficiency of them.3. The parallel and iterative solver HYPRE is applied to the decoupled scheme to evaluate theparallel performance of the simulator.In this experiment, matrix acidization is done in a rectangular domain with the size of 0.1-m by0.04-m. Acid flow is injected into the domain from the left side with the velocity of −6 m/s. Apressure of Pa is imposed on the right side of the domain. The up and down boundaries ofthe domain are closed. The concentration of the injected acid flow is mol/m . The averageinitial porosity of the matrix is 0.18, and the average initial permeability of the matrix is −16 m . The simulator is written in FORTRAN 90 and MPI, and runs on macOS Catalina with four cores. Ifa grid of 180 ×
72 cells is imposed on the domain, and the time step is 10 s, the simulation results areshown as below. Table 2 Results comparison
Improved DBF framework Decoupled schemePore volume to breakthrough
Average porosity
Average concentration Solver time
For Task 1, the key parameters at breakthrough are compared in Table 2. It can be seen that the porevolume to breakthrough, average porosity and average concentration in the domain are more or less thesame for the improved DBF framework and decoupled scheme. Moreover, the porosity profile,concentration profile and the streamlines are also compared in Figure 3. From the table and the figure, itcan be said that both the improved DBF framework and the decoupled scheme simulate very similarresults. In such condition, the correctness of the decoupled scheme is demonstrated by the experiment.Improved DBF framework Decoupled schemeFigure 3 The porosity profile, concentration profile and streamlines in the domain for the improved DBFframework (left column) and decoupled scheme (right column).For Task 2, it can be known from the last row of Table 2 that the solver time of the decoupledscheme is much less than the counterpart of the improved DBF framework, which demonstrates that thedecoupled scheme has much higher computing efficiency than the improved DBF framework.For Task 3, since there are only four cores in our Mac, the simulator will use one, two, and fourprocessors to run, respectively. GMRES [11] in HYPRE is used to solve all the three linear systems. Thesolver time is 2267.0 s, 1141.4 s and 600.7 s, respectively, for the processors above. From the result, itcan be known that the parallel and iterative solver HYPRE has a good scalability on the decoupledscheme.
4. Conclusion
Accurate simulation of matrix acidization requires the change of porosity be well captured by thesimulator, which can be guaranteed by the use of fine enough grids. However, such kind of grids makesthis simulation fall in the field of large-scale simulation, and the application of parallel and iterativesolvers becomes mandatory. The former numerical scheme of the improved DBF framework generated acoupled pressure-velocity linear system, and both pressure and velocity were solved together.Considering that zeros appear in the diagonal of the coefficient matrix, the coupled linear system canonly be solved by direct solvers, which leads to huge running time of the large-scale simulation asmentioned above. Thus, in this work, a decoupled scheme is provided to decouple the coupledpressure-velocity linear system into two linear systems: one is to solve for pressure, and the other one isto solve for velocity. By this way, both of the new linear systems can be solved by parallel and iterativesolvers such as HYPRE. As a result of that, the large-scale simulation can be finished in a reasonabletime period. We also demonstrate mathematically that the decoupled scheme has the same simulationresults as the coupled scheme of the improved DBF framework, except for the pressure, which is alsoproved by the numerical experiment. Besides that, the numerical experiment also shows the highercomputing efficiency of the decoupled scheme, and the good scalability of HYPRE on the decoupledscheme.The next step of this work is to derive the decoupled scheme for the thermal DBF framework [12],which is an expansion of the improved DBF framework by adding the energy conservation equation.However, this kind of derivation is not straightforward. As mentioned above, the pressure p derivedfrom the decoupled scheme is not the true pressure, which leads to that the pressure p cannot be useddirectly in the energy conservation equation of the thermal DBF framework. Some kind of correctionform of the pressure p should be considered. References [1]. Wu, Yuanqing, Amgad Salama, and Shuyu Sun. "Parallel simulation of wormhole propagation withthe Darcy–Brinkman–Forchheimer framework."
Computers and Geotechnics
69 (2015): 564-577.[2]. Wu, Yuanqing.
Parallel Reservoir Simulations with Sparse Grid Techniques and Applications toWormhole Propagation . Diss. 2015.[3]. Kou, Jisheng, Shuyu Sun, and Yuanqing Wu. "Mixed finite element-based fully conservativemethods for simulating wormhole propagation."
Computer Methods in Applied Mechanics andEngineering
298 (2016): 279-302.[4]. Kou, Jisheng, Shuyu Sun, and Yuanqing Wu. "A semi-analytic porosity evolution scheme forsimulating wormhole propagation with the Darcy–Brinkman–Forchheimer model."
Journal ofComputational and Applied Mathematics
348 (2019): 401-420.[5]. Wu, Yuanqing, and Maoqing Ye. "A Newton’s Second Law Abided Darcy-Brinkman-ForchheimerFramework in Matrix Acidization Simulation."
International Conference on Computational &Experimental Engineering and Sciences . Springer, Cham, 2019.[6]. Panga, Mohan KR, Murtaza Ziauddin, and Vemuri Balakotaiah. "Two ‐ scale continuum model forsimulation of wormholes in carbonate acidization." AIChE journal ‐ performancecomputers." Concurrency: Practice and Experience
ACM Transactions on Mathematical Software (TOMS)
International Conference on Computational Science . Springer, Berlin, Heidelberg,2002.[10]. Amestoy, Patrick R., Iain S. Duff, and J-Y. L'excellent. "Multifrontal parallel distributedsymmetric and unsymmetric solvers."
Computer methods in applied mechanics andengineering