2 The problem and some preliminaries

The equation system that simulates wormhole propagation in porous media includes the momentum balance equation (2.1), the mass conservation equation (2.2), the solute transport equation (2.3) and a number of semiempirical equations that bridge the scale differences, as described earlier (2.4) and (2.5) (see, for example, Panga et al. (Reference Panga, Ziauddin and Balakotaiah2005), Wu et al. (Reference Wu, Salama and Sun2015), Kou et al. (Reference Kou, Sun and Wu2016), Li & Rui (Reference Li and Rui2017a )). They can be listed as follows:

with $F(\unicode[STIX]{x1D719})=1.75/\sqrt{150\unicode[STIX]{x1D719}^{3}}$ being the Forchheimer coefficient,

In the above equations, $\unicode[STIX]{x1D6FA}=(L_{lx},L_{rx})\times (L_{ly},L_{ry})$ . $J=(0,T]$ and $T$ denotes the final time. $\boldsymbol{u}$ is the velocity vector, $\unicode[STIX]{x1D70C}$ is the fluid density, $\unicode[STIX]{x1D719}$ is the porosity and $\unicode[STIX]{x1D707}$ is the fluid viscosity. $K$ is the permeability, $\boldsymbol{g}$ is the gravity vector and $\unicode[STIX]{x1D700}$ is a pseudo-compressibility parameter that results in a slight change of the density of the fluid phase in the dissolution process. As stated in Wu et al. (Reference Wu, Salama and Sun2015), in some cases of $\unicode[STIX]{x1D700}=0$ , the linear system may be singular and cannot be solved by the linear solver, and as a result, it is suggested that $\unicode[STIX]{x1D700}$ be a very small positive number to ensure an invertible coefficient matrix. $f=f_{I}+f_{P}$ , where $f_{P}$ and $f_{I}$ are the production and injection rates, respectively. $p$ is the pressure and $c_{f}$ is the cup-mixing concentration of the acid in the fluid phase. For simplicity we assume that the diffusion coefficient $\unicode[STIX]{x1D63F}=d_{mol}\unicode[STIX]{x1D644}=\text{diag}(D_{ll}),~l=1,2$ is a diagonal matrix in the following. $c_{I}$ is the injected concentration. The variable $c_{s}$ is the concentration of the acid at the fluid–solid interface. $\unicode[STIX]{x1D6FC}$ is the dissolving power of the acid and $\unicode[STIX]{x1D70C}_{s}$ is the density of the solid phase. $k_{c}$ is the local mass-transfer coefficient. $k_{s}$ is the surface reaction rate constant.

The above system is solved simultaneously to obtain the pressure, velocity and concentration of the acid. The term $-\unicode[STIX]{x1D735}p-(\unicode[STIX]{x1D707}/K)\boldsymbol{u}+\unicode[STIX]{x1D70C}\boldsymbol{g}$ is called the Darcy term. The term $(\unicode[STIX]{x1D707}/\unicode[STIX]{x1D719})\unicode[STIX]{x0394}\boldsymbol{u}$ is the Brinkman term, which is used to account for transitional flow between boundaries. The inertial term $(\unicode[STIX]{x1D70C}F/\sqrt{K})|\boldsymbol{u}|\boldsymbol{u}$ , which is called the Forchheimer term, is included because the inertial effects could become significant when the velocity is high in (2.1). In (2.3), we assume a first-order surface reaction (see, for example, Panga et al. (Reference Panga, Ziauddin and Balakotaiah2005), Kou et al. (Reference Kou, Sun and Wu2016)). The three terms on left-hand side represent the accumulation, convection and dispersion of the acid, respectively. The first term on the right-hand side denotes the transfer of acid species from the fluid phase to the fluid–solid interface (see, for example, Panga et al. (Reference Panga, Ziauddin and Balakotaiah2005), Kou et al. (Reference Kou, Sun and Wu2016)). It should be noted that temperature effects in Zhao, Hobbs & Ord (Reference Zhao, Hobbs and Ord2014) have been ignored in (2.1)–(2.5), so that only an isothermal acid dissolution system is studied in this model.

In the pore-scale model, the following two equations are given:

where $a_{v}$ is the interfacial area and the subscript $0$ stands for the initial value. Equation (2.6) is the Carman–Kozeny correction, which establishes the relationship between the permeability and the porosity. Equation (2.6) shows that the interfacial area can be calculated from porosity and permeability.

The boundary and initial conditions are as follows.

where $\boldsymbol{n}$ is the unit outward normal vector of the domain $\unicode[STIX]{x1D6FA}$ . For the physical quantities, we assume that $0<\unicode[STIX]{x1D719}_{0\ast }\leqslant \unicode[STIX]{x1D719}_{0}\leqslant \unicode[STIX]{x1D719}_{0}^{\ast }<1$ , $0<K_{0\ast }\leqslant K_{0}\leqslant K_{0}^{\ast }$ and $0\leqslant a_{0\ast }\leqslant a_{0}\leqslant a_{0}^{\ast }$ , where $g_{\ast }$ and $g^{\ast }$ denote the upper and lower bounds of $g$ , respectively. In addition, we assume that $0\leqslant c_{f0},~c_{I}\leqslant 1$ .

We consider the MAC scheme for the model problem. We use the partitions and notations as follows.

The two-dimensional domain $\unicode[STIX]{x1D6FA}$ is partitioned as $\unicode[STIX]{x1D6FF}_{x}\times \unicode[STIX]{x1D6FF}_{y}$ , where

For simplicity we also use the following notations.

For possible integers $i,j$ , $0\leqslant i\leqslant N_{x},~0\leqslant j\leqslant N_{y}$ , define

It is clear that

For a given partition, let $C_{0}$ be a positive number such that

For a function $f(x,y)$ , let $f_{l,m}$ denote $f(x_{l},y_{m})$ , where $l$ may take values $i,~i+1/2$ for integer $i$ , and $m$ may take values $j,~j+1/2$ for integer $j$ . For discrete functions with values at appropriate nodal points, define

For functions $f$ and $g$ , define some discrete $l^{2}$ inner products and norms as follows.

Further define discrete $l^{2}$ inner products and norms as follows.

For vector-valued functions $\boldsymbol{u}=(u^{x},u^{y})$ , it is clear that

and $\Vert d_{y}u^{y}\Vert _{l^{2},M},\Vert D_{x}u^{y}\Vert _{l^{2},T_{x}}$ can be represented similarly. Finally, define the discrete $H^{1}$ -norm and discrete $l^{2}$ -norm of a vectored-valued function $\boldsymbol{u}$ ,

Next we give some definitions and lemmas which will be used in error estimates. For $0\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$ we define the discrete function

The next lemma can be proved as in Weiser & Wheeler (Reference Weiser and Wheeler1988).

Define the discrete interpolations of $u^{x}$ and $u^{y}$ in $\unicode[STIX]{x1D6FA}$ and on the boundary $\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}$ as follows.

Extending the value of $\tilde{u} ^{x}$ and $\tilde{u} ^{y}$ in the the following way,

Now we estimate $(d_{x}\tilde{u} _{i+1/2,j+1/2}^{x}-d_{x}\tilde{u} _{i-1/2,j+1/2}^{x})/h_{i}$ . For $(x_{i+1/2},y_{j+1/2})\in \unicode[STIX]{x1D6FA}$ set

Next we estimate $(D_{y}\tilde{u} _{i,j+1}^{x}-D_{y}\tilde{u} _{i,j}^{x})/k_{j+1/2}$ . Set

According to (2.10) we have

Similarly to the definitions of (2.44) and (2.47) we set

Then there hold

Then similarly to Lemmas 4 and 5 we can obtain the next two lemmas.

Finally, we make the following assumptions which will be used in the rest of the paper.

1. (H1) $p\in W_{\infty }^{2}(J;W_{\infty }^{2}(\unicode[STIX]{x1D6FA}))$ , $\boldsymbol{u}\in W_{\infty }^{2}(J;W_{\infty }^{4}(\unicode[STIX]{x1D6FA}))^{2}$ .

2. (H2) $c_{f}\in W_{\infty }^{2}(J;W_{\infty }^{4}(\unicode[STIX]{x1D6FA}))$ , $\unicode[STIX]{x1D719}\in W_{\infty }^{2}(J;L^{\infty }(\unicode[STIX]{x1D6FA}))$ .

3. (H3) The coefficients $K,\unicode[STIX]{x1D70C},F$ are smooth functions that satisfy $\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})/\sqrt{K(\unicode[STIX]{x1D719})}\in C^{1}(\unicode[STIX]{x1D6FA})$ and there exist positive constants $C_{1}$ and $C_{2}$ such that

   (2.58a,b ) $$\begin{eqnarray}C_{1}\leqslant \frac{\unicode[STIX]{x1D707}}{K}\leqslant C_{2},\quad 0\leqslant \frac{\unicode[STIX]{x1D70C}F}{\sqrt{K}}\leqslant C_{2}.\end{eqnarray}$$

4. (H4) There exist two positive constants $D_{\ast }$ and $D^{\ast }$ such that

   (2.59) $$\begin{eqnarray}D_{\ast }\leqslant D_{ll}(\boldsymbol{x})\leqslant D^{\ast },\quad \text{for }\boldsymbol{x}\in \unicode[STIX]{x1D6FA},l=1,2.\end{eqnarray}$$

5. (H5) $\unicode[STIX]{x1D700},\unicode[STIX]{x1D6FC},\unicode[STIX]{x1D70C}_{s},\unicode[STIX]{x1D707},k_{c}$ and $k_{s}$ are all given positive constants.

3 A fully conservative finite difference algorithm on non-uniform staggered grids

In this section, in order to solve the nonlinear system (2.1)–(2.6) efficiently, the fully conservative finite difference algorithm on non-uniform staggered grids is considered.

Define $\boldsymbol{w}=(w^{x},w^{y})=\boldsymbol{u}c_{f}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D63F}\unicode[STIX]{x1D735}c_{f}$ . Then using (2.5)–(2.6), (2.3) and (2.4) can be transformed into

where $r_{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D6FC}a_{0}k_{c}k_{s}/\unicode[STIX]{x1D70C}_{s}(k_{c}+k_{s}),$ and noting the assumption of $a_{0}$ , we obtain $0\leqslant r_{\unicode[STIX]{x1D719}\ast }\leqslant r_{\unicode[STIX]{x1D719}}\leqslant r_{\unicode[STIX]{x1D719}}^{\ast }.$

For the definition of the scheme, we make some preparations. For a discrete function $q_{i+1/2,j+1/2}$ , define a piecewise-constant function on $\unicode[STIX]{x1D6FA}$ such that

Furthermore, for a pair of discrete functions $\{V_{i,j+1/2}^{x}\}$ and $\{V_{i+1/2,j}^{y}\}$ , define the interpolant operator $\unicode[STIX]{x1D644}_{1}$ as follows:

where

Let $N(U,V)$ denote the norm function for a vector $(U,V)$ ,

Then define the square root average operators $S_{x}$ and $S_{y}$ as

Direct calculation shows that

Define the cutoff operator $\unicode[STIX]{x1D712}$ for velocity as

where $M$ is a given large positive constant.

Denote by $\{Z^{n},\boldsymbol{W}^{n},P^{n},\boldsymbol{U}^{n}\}_{n=1}^{N}$ , the finite difference approximations to $\{\!c_{f}^{n},\boldsymbol{w}^{n},p^{n}$ , $\boldsymbol{u}^{n}\!\}_{n=1}^{N}$ , respectively. Their values are defined as follows.

where ${\mathcal{P}}_{h}$ is a bilinear interpolant operator which is defined as follows: for points $(x,y)$ , assuming $x\in [x_{i-1/2},x_{i+1/2}],y\in [y_{j-1/2},y_{j+1/2}]$ , then ${\mathcal{P}}_{h}\unicode[STIX]{x1D713}(x,y)$ can be evaluated by

For the calculation of the discrete porosity $\unicode[STIX]{x1D713}$ , we use the following scheme.

where $\overline{Z}^{n-1}=\max \{0,\min \{Z^{n-1},1\}\}.$

Set the boundary and initial approximations as follows:

The difference method will consist of three parts: firstly, if the approximate concentration $Z_{i+1/2,j+1/2}^{n-1}$ and porosity $\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n-1},~n=1,\ldots ,N$ are known, equation (3.19) will be used to obtain a new porosity $\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}$ . Since (3.15)–(3.17) form a nonlinear finite difference scheme. By using Picard iteration, an approximation $P_{i+1/2,j+1/2}^{n}$ to the pressure and $U_{i,j+1/2}^{x,n}$ , $U_{i+1/2,j}^{y,n}$ to the velocity will be calculated using $\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}$ . Finally, a new concentration $Z_{i+1/2,j+1/2}^{n}$ will be calculated using $U_{i,j+1/2}^{x,n}$ , $U_{i+1/2,j}^{y,n}$ and $\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}$ , and we get the approximations $W_{i,j+1/2}^{x,n}$ and $W_{i+1/2,j}^{y,n}$ . It is easy to see that at each time level, the difference scheme has an explicit solution or is a linear system with a strictly diagonally dominant coefficient matrix – thus the approximate solutions exist uniquely.

4 Error analysis for the discrete scheme

In this section, we give the error estimates of (3.12)–(3.19). Set

Next we will prove a priori bounds for the discrete solution $\unicode[STIX]{x1D713}$ which will be used in what follows.

Based on the above lemmas, we can obtain

Since $(U^{x}-\tilde{u} ^{x})_{i,j+1/2}=0$ for $i=0,N_{x}$ and $(U^{y}-\tilde{u} ^{y})_{i+1/2,j}=0$ for $j=0,N_{y}$ , we know that

Recalling the definition of $\tilde{\boldsymbol{u}}$ , we have the main error estimates below.

5 Numerical examples

In this section, some numerical experiments using the fully conservative finite difference method on non-uniform staggered grids have been carried out.

We test examples 1 and 2 to verify the convergence rates. The time step is refined as $\unicode[STIX]{x0394}t=1/N_{x}^{2}=1/N_{y}^{2}$ to show the convergence and

where $\unicode[STIX]{x1D644}$ is an identity matrix. In Example 1, set $\unicode[STIX]{x1D6FA}=(0,1)\times (0,1)$ . In Example 2, set $\unicode[STIX]{x1D6FA}=(-1,1)\times (-1,1)$ . The initial spatial partition is a $5\times 5$ grid. Then the grid is refined by dividing each edge into two equal parts and non-uniform meshes are used which are generated from the corresponding uniform mesh by adding a random displacement in both $x$ and $y$ directions (see figure 1). The numerical results of Example 1 are listed in tables 1–8, where tables 1–3 present the results with uniform grids and tables 4–8 demonstrate the results with non-uniform grids. The numerical results of Example 2 are listed in tables 9–16, where tables 9–11 present the results with uniform grids and tables 12–16 demonstrate the results with non-uniform grids. Here on uniform grids we do not present the results for $u^{y}$ since its results are similar to $u^{x}$ . In tables 4–8, 12–16, the ratios between the maximum and minimum mesh sizes in each dimension (that is, $h/h_{min}$ and $k/k_{min}$ ) are listed to illustrate the non-uniformity of the grid, where $h_{min}=\min _{1\leqslant i\leqslant N_{x}}\{h_{i-1/2}\},~k_{min}=\min _{1\leqslant j\leqslant N_{y}}\{k_{j-1/2}\}$ . For simplicity, in the following tables, we set $E_{\boldsymbol{u}}=\boldsymbol{U}-\boldsymbol{u}$ and define

where $\unicode[STIX]{x1D709}=M,T,T_{x},T_{y}$ .

Figure 1. The non-uniform mesh generated from the $20\times 20$ uniform mesh for Example 2.

Figure 2. Numerical porosity of Example 3.

Figure 3. Numerical porosity of Example 4.

Figure 4. Numerical concentration of Example 4.

Figure 5. Numerical pressure of Example 4.

Figure 6. Numerical velocity of Example 4.

In the following examples 3 and 4, we give the dissolution simulations, which are similar to Kou et al. (Reference Kou, Sun and Wu2016), Wu et al. (Reference Wu, Salama and Sun2015). $\unicode[STIX]{x1D6FA}=(0,0.2)\times (0,0.2)$ , $J=[0,1\times 10^{6}]$ , $\unicode[STIX]{x0394}t=1\times 10^{4}$  s. The initial concentration of acid flow is zero. The initial pressure is $10^{5}$  Pa. The parameter $\unicode[STIX]{x1D700}$ is set to 0.01. In addition, some properties of acid flow in porous media, the injection and production flow rates are listed as follows: