1 Introduction
As a result of the injection of acids into a acid dissolution system, a worm-like hole is generated and developed in the subsurface formations, which is called a wormhole. Recently, researchers have paid increasing attention to wormhole formation during reactive dissolution of carbonates, which plays a significant role in enhancing oil production rates when muds and fines deposit at the perforated well pore pipe (see, for example, Fredd & Fogler (Reference Fredd and Fogler1998), Golfier et al. (Reference Golfier, Zarcone, Bazin, Lenormand, Lasseux and Quintard2002), Panga, Ziauddin & Balakotaiah (Reference Panga, Ziauddin and Balakotaiah2005), Liu, Zhang, Mou & Zhou (Reference Liu, Zhang, Mou and Zhou2013), Wu, Salama & Sun (Reference Wu, Salama and Sun2015), Kou, Sun & Wu (Reference Kou, Sun and Wu2016), Li & Rui (Reference Li and Rui2017a )). By increasing the permeability of the damaged zone near the well bore, the primary issue of this process is to increase the production rate. Wormholes are channels with high porosity which establish an effective connectivity between the reservoir and the well. These channels are created by the injected acid, which dissolves the material near the well bore. To simulate wormhole generation and propagation, a combination of the finite-element method and the finite difference method was developed in Zhao et al. (Reference Zhao, Hobbs, Hornby, Ord, Peng and Liu2008). Kou et al. (Reference Kou, Sun and Wu2016) researched mixed finite-element-based fully conservative methods to simulate wormhole propagation. To simulate the wormhole formation procedure, the new Darcy–Brinkman–Forchheimer (DBF) framework was presented in Wu et al. (Reference Wu, Salama and Sun2015). Recently, Li & Rui (Reference Li and Rui2017a , Reference Li and Rui2018a ) have proposed and analysed a block-centred finite difference method to solve both incompressible and compressible wormhole propagation.
Although a large number of numerical methods have been applied to research wormhole generation and propagation in porous media, most of them are based on the Darcy framework. Although this may be adequate when there is no dramatic change in porosity, for scenarios where there is a significant change in porosity (and consequently the permeability) it can be considerably unreliable. Actually, the fluid velocity becomes high in the high-porosity area in cases where the porosity is large (but less than one) and its spatial distribution also becomes very heterogeneous. Also there are scenarios where the porosity equals unity, which means an area of clear fluid within the porous media. Thus, an alternative practical framework that accounts for both the porous media and clear fluid area should be proposed. This is the Darcy–Brinkman–Forchheimer (DBF) framework (see, for example, Vafai & Kim (Reference Vafai and Kim1995), Rees (Reference Rees1999), Alazmi & Vafai (Reference Alazmi and Vafai2000), Wu et al. (Reference Wu, Salama and Sun2015)). In this work, we propose the DBF framework to describe matrix acidization for the purpose of obtaining more reliable results. In comparisons with experimental data under conditions of variable porosity in Kladias & Prasad (Reference Kladias and Prasad1991), the DBF framework predicts the convective flow reasonably well. Besides, the DBF framework accounts for boundary-layer development, macroscopic and microscopic shear stress, and the initial force. A porous medium and clear fluid interface have been extensively studied using the DBF framework and continuity of velocities and stress at the interface in Vafai & Kim (Reference Vafai and Kim1995).
Among the many numerical methods that are available for solving the Stokes and Navier–Stokes equations, one of the best and simplest is the marker and cell method (MAC). The MAC scheme has the ability to enforce the incompressibility constraint of the velocity field in a pointwise fashion. Moreover, it has been shown to locally conserve the mass, momentum and kinetic energy (Perot Reference Perot2000, Reference Perot2011). The MAC method is a class of finite volume methods on rectangular cells, with pressure approximated at the cell centre, the
$x$
-component of velocity approximated at the midpoint of the vertical edges of the cell, and the
$y$
-component of velocity approximated at the midpoint of the horizontal edges of the cell.
The MAC method can also be interpreted as a mixed finite-element method coupled with a quadrature formula. This interpretation can be found in Girault & Raviart (Reference Girault and Raviart1979, Reference Girault and Raviart2012), where the mixed method is analysed. Error estimates for the MAC method were also demonstrated in Girault & Lopez (Reference Girault and Lopez1996). In 1998, Han & Wu (Reference Han and Wu1998) proposed that the MAC scheme can be obtained from a new mixed finite-element method. Kanschat (Reference Kanschat2008) presented that discontinuous Galerkin discretizations in their lowest-order version were very similar to the MAC finite difference scheme and proposed a way to interpolate solutions obtained by the MAC scheme to obtain a pointwise divergence-free function for the Stokes equations. Inspired by the work of Kanschat (Reference Kanschat2008), in Minev (Reference Minev2008) it was demonstrated that some well-known finite difference schemes can be interpreted within the framework of local discontinuous Galerkin (LDG) methods using low-order piecewise solenoidal discrete spaces. Recently, second-order superconvergence for both velocity and pressure of the MAC scheme for Stokes and Navier–Stokes equations on non-uniform grids was derived rigorously (Rui & Li Reference Rui and Li2017; Li & Rui Reference Li and Rui2018b ).
Block-centred finite differences, sometimes called cell-centred finite differences, can be thought as the lowest-order Raviart–Thomas mixed element method (Raviart & Thomas Reference Raviart and Thomas1977), with a suitable quadrature formulation. The mixed finite elements for elliptic problems with tensor coefficients as cell-centred finite differences was proposed in Arbogast, Wheeler & Yotov (Reference Arbogast, Wheeler and Yotov1997). Rui & Pan (Reference Rui and Pan2012) demonstrated a block-centred finite difference method for the Darcy–Forchheimer model. Block-centred finite difference methods have been developed (Rui & Pan Reference Rui and Pan2013; Li & Rui Reference Li and Rui2016; Liu, Cui & Xin Reference Liu, Cui and Xin2018) and it has been applied to the nonlinear time-fractional parabolic equation by Li & Rui (Reference Li and Rui2017b ). Applications of block-centred finite difference methods enable us to approximate pressure, velocity, porosity, concentration and its flux in different discrete norms with second-order accuracy on non-uniform grids, which is superconvergence. Also block-centred finite difference methods can guarantee local mass conservation. Moreover, applications of block-centred finite difference methods enable us to transfer the saddle point problem to symmetric positive definite problems.
Both the MAC and block-centred finite difference methods are based on finite difference methods on non-uniform staggered grids. Thus, in this paper, the finite difference method on non-uniform staggered grids is derived for wormhole propagation with the Darcy–Brinkman–Forchheimer framework in porous media. Actually, in our proposed scheme, there are at least three issues in the theoretical analysis for wormhole propagation with the Darcy–Brinkman–Forchheimer framework: estimating and bounding the porosity, which can change with time; complications resulting from the introduction of an auxiliary flux variable to guarantee full mass conservation; and the fully coupling relation of multivariables. To resolve these issues, we introduce some useful lemmas and consider the coupled analysis method. To guarantee full mass conservation, we introduce an auxiliary flux variable. Based on the cutoff operator of velocity and the coupled analysis approach, and inspired by the analysis technique in Rui, Zhao & Pan (Reference Rui, Zhao and Pan2015), Rui & Li (Reference Rui and Li2017) for the Darcy–Forchheimer and Stokes equations, we obtain second-order superconvergence of the pressure, velocity, porosity, concentration and its flux in different discrete norms on non-uniform grids. We also obtain second-order superconvergence for some terms of the
$H^{1}$
norm of the velocity on non-uniform grids. Finally, some numerical experiments are presented to verify the theoretical analysis and effectiveness of the given scheme.
The paper is organized as follows. In § 2 we give the problem and some preliminaries. In § 3 we present the fully conservative finite difference algorithm on non-uniform staggered grids. In § 4 we demonstrate error analysis for the discrete scheme. In § 5 we carry out some numerical experiments using the finite difference scheme on non-uniform staggered grids to verify the theoretical analysis and effectiveness of the proposed scheme.
Throughout the paper we use
$C$
, with or without a subscript, to denote a positive constant, which can have different values depending on the circumstances in which it is used.
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:
$$\begin{eqnarray}\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}\boldsymbol{u}}{\unicode[STIX]{x2202}t}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}\unicode[STIX]{x0394}\boldsymbol{u}+\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}|\boldsymbol{u}|\boldsymbol{u}=-\unicode[STIX]{x1D735}p-\frac{\unicode[STIX]{x1D707}}{K(\unicode[STIX]{x1D719})}\boldsymbol{u}+\unicode[STIX]{x1D70C}\boldsymbol{g},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},t\in J\end{eqnarray}$$
with
$F(\unicode[STIX]{x1D719})=1.75/\sqrt{150\unicode[STIX]{x1D719}^{3}}$
being the Forchheimer coefficient,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{u}=f,\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},t\in J & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D719}c_{f})}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\boldsymbol{u}c_{f})-\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\unicode[STIX]{x1D63F}\boldsymbol{\cdot }\unicode[STIX]{x1D735}c_{f})=k_{c}a_{v}(c_{s}-c_{f})+f_{P}c_{f}+f_{I}c_{I},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},t\in J & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\frac{\unicode[STIX]{x1D6FC}k_{c}a_{v}(c_{f}-c_{s})}{\unicode[STIX]{x1D70C}_{s}},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},t\in J & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle c_{s}=\frac{c_{f}}{1+k_{s}/k_{c}},\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},t\in J. & \displaystyle\end{eqnarray}$$
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:
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}l@{}}\displaystyle \frac{K}{K_{0}}=\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{0}}\left(\frac{\unicode[STIX]{x1D719}(1-\unicode[STIX]{x1D719}_{0})}{\unicode[STIX]{x1D719}_{0}(1-\unicode[STIX]{x1D719})}\right)^{2},\\ \displaystyle \frac{a_{v}}{a_{0}}=\frac{\unicode[STIX]{x1D719}}{\unicode[STIX]{x1D719}_{0}}\sqrt{\frac{K_{0}\unicode[STIX]{x1D719}}{K\unicode[STIX]{x1D719}_{0}}},\end{array}\right\} & \displaystyle\end{eqnarray}$$
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.
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\boldsymbol{u}=0,~(\boldsymbol{u}c_{f}-\unicode[STIX]{x1D719}\unicode[STIX]{x1D63F}\unicode[STIX]{x1D735}c_{f})\boldsymbol{\cdot }\boldsymbol{n}=0,\quad \boldsymbol{x}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA},t\in J,\displaystyle \\ p(\boldsymbol{x},0)=p_{0}(\boldsymbol{x}),\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},\displaystyle \\ \boldsymbol{u}(\boldsymbol{x},0)=\boldsymbol{u}_{0}(\boldsymbol{x}),\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},\displaystyle \\ c_{f}(\boldsymbol{x},0)=c_{f0}(\boldsymbol{x}),\quad \boldsymbol{x}\in \unicode[STIX]{x1D6FA},\end{array}\right\}\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}_{x}:L_{lx}=x_{0}<x_{1}<\cdots <x_{N_{x}-1}<x_{N_{x}}=L_{rx}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}_{y}:L_{ly}=y_{0}<y_{1}<\cdots <y_{N_{y}-1}<y_{N_{y}}=L_{ry}. & \displaystyle\end{eqnarray}$$
For simplicity we also use the following notations.
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}x_{-1/2}=x_{0}=L_{lx},\quad x_{N_{x}+1/2}=x_{N_{x}}=L_{rx},\\ y_{-1/2}=y_{0}=L_{ly},\quad y_{N_{y}+1/2}=y_{N_{y}}=L_{ry}.\end{array}\right\}\end{eqnarray}$$
For possible integers
$i,j$
,
$0\leqslant i\leqslant N_{x},~0\leqslant j\leqslant N_{y}$
, define
$$\begin{eqnarray}\displaystyle & \displaystyle x_{i+1/2}=\frac{x_{i}+x_{i+1}}{2},\quad h_{i+1/2}=x_{i+1}-x_{i},\quad h=\max _{i}h_{i+1/2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle h_{i}=x_{i+1/2}-x_{i-1/2}=\frac{h_{i+1/2}+h_{i-1/2}}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle y_{j+1/2}=\frac{y_{j}+y_{j+1}}{2},\quad k_{j+1/2}=y_{j+1}-y_{j},\quad k=\max _{j}k_{j+1/2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle k_{j}=y_{j+1/2}-y_{j-1/2}=\frac{k_{j+1/2}+k_{j-1/2}}{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FA}_{i+1/2,j+1/2}=(x_{i},x_{i+1})\times (y_{j},y_{j+1}). & \displaystyle\end{eqnarray}$$
It is clear that
$$\begin{eqnarray}\displaystyle & \displaystyle h_{0}=\frac{h_{1/2}}{2},\quad h_{N_{x}}=\frac{h_{N_{x}-1/2}}{2},\quad k_{0}=\frac{k_{1/2}}{2},\quad k_{N_{y}}=\frac{k_{N_{y}-1/2}}{2}. & \displaystyle\end{eqnarray}$$
For a given partition, let
$C_{0}$
be a positive number such that
$$\begin{eqnarray}\min _{i,j}\{h_{i+1/2},k_{j+1/2}\}\geqslant C_{0}\max _{i,j}\{h_{i+1/2},k_{j+1/2}\}.\end{eqnarray}$$
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
$$\begin{eqnarray}\left.\begin{array}{@{}cc@{}}\displaystyle [d_{x}\,f]_{i+1/2,m}=\frac{f_{i+1,m}-f_{i,m}}{h_{i+1/2}}, & \displaystyle [D_{y}\,f]_{l,j+1}=\frac{f_{l,j+3/2}-f_{l,j+1/2}}{k_{j+1}},\\ \displaystyle [D_{x}\,f]_{i+1,m}=\frac{f_{i+3/2,m}-f_{i+1/2,m}}{h_{i+1}}, & \displaystyle [d_{y}\,f]_{l,j+1/2}=\frac{f_{l,j+1}-f_{l,j}}{k_{j+1/2}},\\ \displaystyle [d_{t}\,f]^{n}=\frac{f^{n}-f^{n-1}}{\unicode[STIX]{x0394}t}.\end{array}\right\}\end{eqnarray}$$
For functions
$f$
and
$g$
, define some discrete
$l^{2}$
inner products and norms as follows.
$$\begin{eqnarray}\displaystyle (f,g)_{l^{2},M} & \equiv & \displaystyle \mathop{\sum }_{i=0}^{N_{x}-1}\mathop{\sum }_{j=0}^{N_{y}-1}h_{i+1/2}k_{j+1/2}\,f_{i+1/2,j+1/2}g_{i+1/2,j+1/2},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle (f,g)_{l^{2},T} & \equiv & \displaystyle \mathop{\sum }_{i=1}^{N_{x}-1}\mathop{\sum }_{j=1}^{N_{y}-1}h_{i}k_{j}\,f_{i,j}g_{i,j},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle (f,g)_{l^{2},T_{x}} & \equiv & \displaystyle \mathop{\sum }_{i=0}^{N_{x}}\mathop{\sum }_{j=1}^{N_{y}-1}h_{i}k_{j}\,f_{i,j}g_{i,j},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle (f,g)_{l^{2},T_{y}} & \equiv & \displaystyle \mathop{\sum }_{i=1}^{N_{x}-1}\mathop{\sum }_{j=0}^{N_{y}}h_{i}k_{j}\,f_{i,j}g_{i,j},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \Vert \,f\Vert _{l^{2},\unicode[STIX]{x1D709}}^{2} & \equiv & \displaystyle (f,f)_{l^{2},\unicode[STIX]{x1D709}},\quad \unicode[STIX]{x1D709}=M,T,T_{x},T_{y}.\end{eqnarray}$$
Further define discrete
$l^{2}$
inner products and norms as follows.
$$\begin{eqnarray}\displaystyle & \displaystyle (f,g)_{l^{2},T,M}\equiv \mathop{\sum }_{i=1}^{N_{x}-1}\mathop{\sum }_{j=0}^{N_{y}-1}h_{i}k_{j+1/2}\,f_{i,j+1/2}g_{i,j+1/2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (f,g)_{l^{2},M,T}\equiv \mathop{\sum }_{i=0}^{N_{x}-1}\mathop{\sum }_{j=1}^{N_{y}-1}h_{i+1/2}k_{j}\,f_{i+1/2,j}g_{i+1/2,j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \Vert \,f\Vert _{l^{2},T,M}^{2}\equiv (f,f)_{l^{2},T,M},\quad \Vert \,f\Vert _{l^{2},M,T}^{2}\equiv (f,f)_{l^{2},M,T}. & \displaystyle\end{eqnarray}$$
For vector-valued functions
$\boldsymbol{u}=(u^{x},u^{y})$
, it is clear that
$$\begin{eqnarray}\displaystyle \Vert d_{x}u^{x}\Vert _{l^{2},M}^{2} & \equiv & \displaystyle \mathop{\sum }_{i=0}^{N_{x}-1}\mathop{\sum }_{j=0}^{N_{y}-1}h_{i+1/2}k_{j+1/2}|d_{x}u_{i+1/2,j+1/2}^{x}|^{2},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \Vert D_{y}u^{x}\Vert _{l^{2},T_{y}}^{2} & \equiv & \displaystyle \mathop{\sum }_{i=1}^{N_{x}-1}\mathop{\sum }_{j=0}^{N_{y}}h_{i}k_{j}|D_{y}u_{i,j}^{x}|^{2},\end{eqnarray}$$
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}$
,
$$\begin{eqnarray}\displaystyle \Vert D\boldsymbol{u}\Vert ^{2} & \equiv & \displaystyle \Vert d_{x}u^{x}\Vert _{l^{2},M}^{2}+\Vert D_{y}u^{x}\Vert _{l^{2},T_{y}}^{2}+\Vert D_{x}u^{y}\Vert _{l^{2},T_{x}}^{2}+\Vert d_{y}u^{y}\Vert _{l^{2},M}^{2}.\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \Vert \boldsymbol{u}\Vert _{l^{2}}^{2} & \equiv & \displaystyle \Vert u^{x}\Vert _{l^{2},T,M}^{2}+\Vert u^{y}\Vert _{l^{2},M,T}^{2}.\end{eqnarray}$$
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
$$\begin{eqnarray}\unicode[STIX]{x1D6FF}_{c_{f},i+1/2,j+1/2}=\frac{h_{i+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}c_{f,i+1/2,j+1/2}}{\unicode[STIX]{x2202}x^{2}}+\frac{k_{j+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}c_{f,i+1/2,j+1/2}}{\unicode[STIX]{x2202}y^{2}}.\end{eqnarray}$$
Lemma 1 (see, for example, Rui & Li (Reference Rui and Li2017)).
If
$c_{f}$
is sufficiently smooth then it holds that
$$\begin{eqnarray}\left.\begin{array}{@{}lll@{}}\displaystyle \frac{\unicode[STIX]{x2202}c_{f,i,j+1/2}}{\unicode[STIX]{x2202}x} & = & \displaystyle D_{x}(c_{f}-\unicode[STIX]{x1D6FF}_{c_{f}})_{i,j+1/2}+\unicode[STIX]{x1D716}_{i,j+1/2}^{x}(c_{f}),\quad (x_{i},y_{j+1/2})\in \unicode[STIX]{x1D6FA},\\ \displaystyle \frac{\unicode[STIX]{x2202}c_{f,i+1/2,j}}{\unicode[STIX]{x2202}y} & = & \displaystyle D_{y}(c_{f}-\unicode[STIX]{x1D6FF}_{c_{f}})_{i+1/2,j}+\unicode[STIX]{x1D716}_{i+1/2,j}^{y}(c_{f}),\quad (x_{i+1/2},y_{j})\in \unicode[STIX]{x1D6FA},\end{array}\right\}\end{eqnarray}$$
with the following approximate properties
$$\begin{eqnarray}\unicode[STIX]{x1D716}_{i,j+1/2}^{x}(c_{f})=O((h^{2}+k^{2})\Vert c_{f}\Vert _{3,\infty }),\quad \unicode[STIX]{x1D716}_{i+1/2,j}^{y}(c_{f})=O((h^{2}+k^{2})\Vert c_{f}\Vert _{3,\infty }).\end{eqnarray}$$
The next lemma can be proved as in Weiser & Wheeler (Reference Weiser and Wheeler1988).
Lemma 2. Let
$\{V_{i,j+1/2}^{x}\},\{V_{i+1/2,j}^{y}\}$
and
$\{q_{i+1/2,j+1/2}^{x}\},\{q_{i+1/2,j+1/2}^{y}\}$
be discrete functions with
$V_{0,j+1/2}^{x}=V_{N_{x},j+1/2}^{x}=V_{i+1/2,0}^{y}=V_{i+1/2,N_{y}}^{y}=0$
, with appropriate integers
$i$
and
$j$
. Then there holds
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}(D_{x}q^{x},V^{x})_{l^{2},T,M}=-(q^{x},d_{x}V^{x})_{l^{2},M},\\ (D_{y}q^{y},V^{y})_{l^{2},M,T}=-(q^{y},d_{y}V^{y})_{l^{2},M}.\end{array}\right\}\end{eqnarray}$$
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.
$$\begin{eqnarray}\displaystyle \displaystyle \tilde{u} _{i,j+1/2}^{x} & = & \displaystyle u_{i,j+1/2}^{x}-\frac{k_{j+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}u_{i,j+1/2}^{x}}{\unicode[STIX]{x2202}y^{2}},\quad 1\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \tilde{u} _{i+1/2,j}^{y} & = & \displaystyle u_{i+1/2,j}^{y}-\frac{h_{i+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}u_{i+1/2,j}^{y}}{\unicode[STIX]{x2202}x^{2}},\quad 0\leqslant i\leqslant N_{x}-1,1\leqslant j\leqslant N_{y}-1.\end{eqnarray}$$
Extending the value of
$\tilde{u} ^{x}$
and
$\tilde{u} ^{y}$
in the the following way,
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{i,-1/2}^{x}=\tilde{u} _{i,0}^{x}=u_{i,0}^{x},\quad \tilde{u} _{i,N_{y}+1/2}^{x}=\tilde{u} _{i,N_{y}}^{x}=u_{i,N_{y}}^{x},\quad 1\leqslant i\leqslant N_{x}-1. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \tilde{u} _{-1/2,j}^{y}=\tilde{u} _{0,j}^{y}=u_{0,j}^{y},\quad \tilde{u} _{N_{x}+1/2,j}^{y}=\tilde{u} _{N_{x},j}^{y}=u_{N_{x},j}^{y},\quad 1\leqslant j\leqslant N_{y}-1. & \displaystyle\end{eqnarray}$$
Lemma 3. Let
$\boldsymbol{u}=(u^{x},u^{y})$
be the solution of (2.1) and (2.2), where
$\tilde{u} ^{x}$
,
$\tilde{u} ^{y}$
are defined by (2.35) and (2.36), then for
$0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1$
there holds
$$\begin{eqnarray}[d_{x}\tilde{u} ^{x}+d_{y}\tilde{u} ^{y}]_{i+1/2,j+1/2}=\frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}y}+r_{i+1/2,j+1/2}^{\boldsymbol{u}},\end{eqnarray}$$
where
$$\begin{eqnarray}r_{i+1/2,j+1/2}^{\boldsymbol{u}}=O((h^{2}+k^{2})(\Vert u^{x}\Vert _{3,\infty }+\Vert u^{y}\Vert _{3,\infty })).\end{eqnarray}$$
Proof. From the definition (2.35) we have that
$$\begin{eqnarray}\displaystyle d_{x}\tilde{u} _{i+1/2,j+1/2}^{x} & = & \displaystyle \displaystyle \frac{1}{h_{i+1/2}}(\tilde{u} _{i+1,j+1/2}^{x}-\tilde{u} _{i,j+1/2}^{x})\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{1}{h_{i+1/2}}\int _{x_{i}}^{x_{i+1}}\left(\frac{\unicode[STIX]{x2202}u^{x}}{\unicode[STIX]{x2202}x}(x,y_{j+1/2})-\frac{k_{j+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{3}u^{x}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y^{2}}(x,y_{j+1/2})\right)\text{d}x\nonumber\\ \displaystyle & = & \displaystyle \displaystyle \frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}x}+O((h^{2}+k^{2})\Vert u^{x}\Vert _{3,\infty }).\end{eqnarray}$$
Similarly
$$\begin{eqnarray}\displaystyle & d_{y}\tilde{u} _{i+1/2,j+1/2}^{y}=\displaystyle \frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}y}+O((h^{2}+k^{2})\Vert u^{y}\Vert _{3,\infty }). & \displaystyle\end{eqnarray}$$
Then adding (2.41) and (2.42) we obtain
$$\begin{eqnarray}[d_{x}\tilde{u} ^{x}+d_{y}\tilde{u} ^{y}]_{i+1/2,j+1/2}=\frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}u_{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}y}+O((h^{2}+k^{2})(\Vert u^{x}\Vert _{3,\infty }+\Vert u^{y}\Vert _{3,\infty })).\end{eqnarray}$$
Setting the last term on the right-hand side of (2.43) as
$r_{i+1/2,j+1/2}^{\boldsymbol{u}}$
completes the proof.
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
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{x,x}=\frac{h_{i+1/2}^{2}}{6}\frac{\unicode[STIX]{x2202}^{3}u_{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}x^{3}}-\frac{k_{j+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{3}u_{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}x\unicode[STIX]{x2202}y^{2}}. & & \displaystyle\end{eqnarray}$$
Lemma 4 (see, for example, Rui & Li (Reference Rui and Li2017)).
Let
$\boldsymbol{u}=(u^{x},u^{y})$
be the solution of (2.1) and (2.2),
$\tilde{u} ^{x}$
is defined by (2.35). When
$u^{x}$
is sufficiently smooth for
$1\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1$
there holds
$$\begin{eqnarray}\displaystyle \frac{d_{x}\tilde{u} _{i+1/2,j+1/2}^{x}-d_{x}\tilde{u} _{i-1/2,j+1/2}^{x}}{h_{i}}=\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{i,j+1/2}^{x}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{x,x}-\unicode[STIX]{x1D716}_{i-1/2,j+1/2}^{x,x}}{h_{i}}+r_{i,j+1/2}^{x,x},\qquad & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & r_{i,j+1/2}^{x,x}=O((h^{2}+k^{2})\Vert u^{x}\Vert _{4,\infty }). & \displaystyle\end{eqnarray}$$
Next we estimate
$(D_{y}\tilde{u} _{i,j+1}^{x}-D_{y}\tilde{u} _{i,j}^{x})/k_{j+1/2}$
. Set
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D716}_{i,j}^{x,y}=-\frac{\unicode[STIX]{x2202}^{3}u_{i,j}^{x}}{\unicode[STIX]{x2202}y^{3}}\frac{k_{j+1/2}^{3}+k_{j-1/2}^{3}}{24k_{j}},\quad 1\leqslant i\leqslant N_{x}-1,1\leqslant j\leqslant N_{y}-1,\\ \displaystyle \unicode[STIX]{x1D716}_{i,0}^{x,y}=-\frac{\unicode[STIX]{x2202}^{3}u_{i,1/2}^{x}}{\unicode[STIX]{x2202}y^{3}}\frac{k_{1/2}^{2}}{12},\quad \unicode[STIX]{x1D716}_{i,N_{y}}^{x,y}=-\frac{\unicode[STIX]{x2202}^{3}u_{i,N_{y}-1/2}^{x}}{\unicode[STIX]{x2202}y^{3}}\frac{k_{N_{y}-1/2}^{2}}{12},1\leqslant i\leqslant N_{x}-1.\end{array}\right\}\end{eqnarray}$$
According to (2.10) we have
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D716}_{i,-1/2}^{x,y}=\unicode[STIX]{x1D716}_{i,0}^{x,y},\quad \unicode[STIX]{x1D716}_{i,N_{y}+1/2}^{x,y}=\unicode[STIX]{x1D716}_{i,N_{y}}^{x,y}. & \displaystyle\end{eqnarray}$$
Lemma 5 (see, for example, Rui & Li (Reference Rui and Li2017)).
Under the condition of Lemma 4, for
$1\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1$
there holds
$$\begin{eqnarray}\displaystyle \frac{D_{y}\tilde{u} _{i,j+1}^{x}-D_{y}\tilde{u} _{i,j}^{x}}{k_{j+1/2}}=\frac{\unicode[STIX]{x2202}^{2}u_{i,j+1/2}^{x}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x1D716}_{i,j+1}^{x,y}-\unicode[STIX]{x1D716}_{i,j}^{x,y}}{k_{j+1/2}}+r_{i,j+1/2}^{x,y}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}r_{i,j+1/2}^{x,y}=O(k^{2}\Vert u^{x}\Vert _{4,\infty }).\end{eqnarray}$$
Similarly to the definitions of (2.44) and (2.47) we set
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{y,y}=-\frac{h_{i+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{3}u_{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}x^{2}\unicode[STIX]{x2202}y}+\frac{k_{j+1/2}^{2}}{6}\frac{\unicode[STIX]{x2202}^{3}u_{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}y^{3}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D716}_{i,j}^{y,x}=-\frac{\unicode[STIX]{x2202}^{3}u_{i,j}^{y}}{\unicode[STIX]{x2202}x^{3}}\frac{h_{i+1/2}^{3}+h_{i-1/2}^{3}}{24h_{i}},\quad 1\leqslant i\leqslant N_{x}-1,~1\leqslant j\leqslant N_{y}-1,\\ \displaystyle \unicode[STIX]{x1D716}_{0,j}^{y,x}=-\frac{\unicode[STIX]{x2202}^{3}u_{1/2,j}^{y}}{\unicode[STIX]{x2202}x^{3}}\frac{h_{1/2}^{2}}{12},~\unicode[STIX]{x1D716}_{N_{x},j}^{y,x}=-\frac{\unicode[STIX]{x2202}^{3}u_{N_{x}-1/2,j}^{y}}{\unicode[STIX]{x2202}x^{3}}\frac{h_{N_{x}-1/2}^{2}}{12},~1\leqslant j\leqslant N_{y}-1.\end{array}\right\} & \displaystyle\end{eqnarray}$$
Then there hold
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D716}_{-1/2,j}^{y,x}=\unicode[STIX]{x1D716}_{0,j}^{y,x},\quad \unicode[STIX]{x1D716}_{N_{x}+1/2,j}^{y,x}=\unicode[STIX]{x1D716}_{N_{x},j}^{y,x}. & \displaystyle\end{eqnarray}$$
Then similarly to Lemmas 4 and 5 we can obtain the next two lemmas.
Lemma 6. Let
$\boldsymbol{u}=(u^{x},u^{y})$
be the solution of (2.1) and (2.2),
$\tilde{u} ^{y}$
is defined by (2.36). When
$u^{y}$
is sufficiently smooth, for
$0\leqslant i\leqslant N_{x}-1,1\leqslant j\leqslant N_{y}-1$
, there holds
$$\begin{eqnarray}\frac{d_{y}\tilde{u} _{i+1/2,j+1/2}^{y}-d_{y}\tilde{u} _{i+1/2,j-1/2}^{y}}{k_{j}}=\displaystyle \frac{\unicode[STIX]{x2202}^{2}u_{i+1/2,j}^{y}}{\unicode[STIX]{x2202}y^{2}}+\frac{\unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{y,y}-\unicode[STIX]{x1D716}_{i+1/2,j-1/2}^{y,y}}{k_{j}}+r_{i+1/2,j}^{y,y},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & r_{i+1/2,j}^{y,y}=O((h^{2}+k^{2})\Vert u^{y}\Vert _{4,\infty }). & \displaystyle\end{eqnarray}$$
Lemma 7. Under the condition of Lemma 6, for
$0\leqslant i\leqslant N_{x}-1,~1\leqslant j\leqslant N_{y}-1$
, there holds
$$\begin{eqnarray}\displaystyle \frac{D_{x}\tilde{u} _{i+1,j}^{y}-D_{x}\tilde{u} _{i,j}^{y}}{h_{i+1/2}}=\frac{\unicode[STIX]{x2202}^{2}u_{i+1/2,j}^{y}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x1D716}_{i+1,j}^{y,x}-\unicode[STIX]{x1D716}_{i,j}^{y,x}}{h_{i+1/2}}+r_{i+1/2,j}^{y,x}, & & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & r_{i+1/2,j}^{y,x}=O(h^{2}\Vert u^{y}\Vert _{4,\infty }). & \displaystyle\end{eqnarray}$$
Finally, we make the following assumptions which will be used in the rest of the paper.
(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}$
.(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}))$
.(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}$$
(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}$$
(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
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}{\displaystyle \frac{\unicode[STIX]{x2202}c_{f}}{\unicode[STIX]{x2202}t}}+c_{f}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{w}+\frac{k_{c}k_{s}}{k_{c}+k_{s}}a_{v}c_{f}-f_{P}c_{f}=f_{I}c_{I}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}=\frac{r_{\unicode[STIX]{x1D719}}(1-\unicode[STIX]{x1D719})}{1-\unicode[STIX]{x1D719}_{0}}c_{f}, & \displaystyle\end{eqnarray}$$
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
$$\begin{eqnarray}\unicode[STIX]{x1D6F1}_{h}q(x,y)=q_{i+1/2,j+1/2},\quad (x,y)\in \unicode[STIX]{x1D6FA}_{i+1/2,j+1/2}.\end{eqnarray}$$
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:
$$\begin{eqnarray}\unicode[STIX]{x1D644}_{1}\boldsymbol{V}=(I_{1}^{x}V^{x},I_{1}^{y}V^{y}),\end{eqnarray}$$
where
$$\begin{eqnarray}\left.I_{1}^{x}V^{x}(x,y)\right|_{\unicode[STIX]{x1D6FA}_{i,j+1/2}}=V_{i,j+1/2}^{x},\quad \left.I_{1}^{y}V^{y}(x,y)\right|_{\unicode[STIX]{x1D6FA}_{i+1/2,j}}=V_{i+1/2,j}^{y}.\end{eqnarray}$$
Let
$N(U,V)$
denote the norm function for a vector
$(U,V)$
,
$$\begin{eqnarray}N(U,V)=(U^{2}+V^{2})^{1/2}.\end{eqnarray}$$
Then define the square root average operators
$S_{x}$
and
$S_{y}$
as
$$\begin{eqnarray}\displaystyle & \displaystyle [S_{x}\boldsymbol{V}]_{i,j+1/2}=\frac{1}{|\unicode[STIX]{x1D6FA}_{i,j+1/2}|}\int _{\unicode[STIX]{x1D6FA}_{i,j+1/2}}N(I_{1}^{x}V^{x},I_{1}^{y}V^{y})\,\text{d}x\,\text{d}y, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle [S_{y}\boldsymbol{V}]_{i+1/2,j}=\frac{1}{|\unicode[STIX]{x1D6FA}_{i+1/2,j}|}\int _{\unicode[STIX]{x1D6FA}_{i+1/2,j}}N(I_{1}^{x}V^{x},I_{1}^{y}V^{y})\,\text{d}x\,\text{d}y. & \displaystyle\end{eqnarray}$$
Direct calculation shows that
$$\begin{eqnarray}\displaystyle [S_{x}\boldsymbol{V}]_{i,j+1/2} & = & \displaystyle \frac{1}{4h_{i}}\{\!h_{i-1/2}[N(V_{i,j+1/2}^{x},V_{i-1/2,j+1}^{y})+N(V_{i,j+1/2}^{x},V_{i-1/2,j}^{y})]\nonumber\\ \displaystyle & & \displaystyle +\,h_{i+1/2}[N(V_{i,j+1/2}^{x},V_{i+1/2,j+1}^{y})+N(V_{i,j+1/2}^{x},V_{i+1/2,j}^{y})]\!\},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle [S_{y}\boldsymbol{V}]_{i+1/2,j} & = & \displaystyle \frac{1}{4k_{j}}\{\!k_{j-1/2}[N(V_{i+1,j-1/2}^{x},V_{i+1/2,j}^{y})+N(V_{i,j-1/2}^{x},V_{i+1/2,j}^{y})]\nonumber\\ \displaystyle & & \displaystyle +\,k_{j+1/2}[N(V_{i+1,j+1/2}^{x},V_{i+1/2,j}^{y})+N(V_{i,j+1/2}^{x},V_{i+1/2,j}^{y})]\!\}.\end{eqnarray}$$
Define the cutoff operator
$\unicode[STIX]{x1D712}$
for velocity as
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D712}(\boldsymbol{u})=\left\{\begin{array}{@{}ll@{}}\boldsymbol{u}, & \text{if }|\boldsymbol{u}|\leqslant M,\\ M\boldsymbol{u}/|\boldsymbol{u}|, & \text{if }|\boldsymbol{u}|>M,\end{array}\right. & \displaystyle\end{eqnarray}$$
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.
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}[d_{t}Z]_{i+1/2,j+1/2}^{n}+Z_{i+1/2,j+1/2}^{n}[d_{t}\unicode[STIX]{x1D713}]_{i+1/2,j+1/2}^{n}+[d_{x}W^{x}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad +\,[d_{y}W^{y}]_{i+1/2,j+1/2}^{n}+\frac{k_{c}k_{s}}{k_{c}+k_{s}}a_{v}(\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})Z_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad -\,[f_{P}Z]_{i+1/2,j+1/2}^{n}=[f_{I}c_{I}]_{i+1/2,j+1/2}^{n},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle W_{i,j+1/2}^{x,n}=\unicode[STIX]{x1D712}(U_{i,j+1/2}^{x,n})[{\mathcal{P}}_{h}Z]_{i,j+1/2}^{n}-D_{11}^{n}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}]_{i,j+1/2}^{n}[D_{x}Z]_{i,j+1/2}^{n},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle W_{i+1/2,j}^{y,n}=\unicode[STIX]{x1D712}(U_{i+1/2,j}^{y,n})[{\mathcal{P}}_{h}Z]_{i+1/2,j}^{n}-D_{22}^{n}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}]_{i+1/2,j}^{n}[D_{y}Z]_{i+1/2,j}^{n},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}[d_{t}P]_{i+1/2,j+1/2}^{n}+[d_{t}\unicode[STIX]{x1D713}]_{i+1/2,j+1/2}^{n}+[d_{x}U^{x}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,[d_{y}U^{y}]_{i+1/2,j+1/2}^{n}=f_{i+1/2,j+1/2}^{n},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D707}K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})+\unicode[STIX]{x1D70C}F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}])_{i,j+1/2}U_{i,j+1/2}^{x,n}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{d_{x}U_{i+1/2,j+1/2}^{x,n}-d_{x}U_{i-1/2,j+1/2}^{x,n}}{h_{i}}-\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{D_{y}U_{i,j+1}^{x,n}-D_{y}U_{i,j}^{x,n}}{k_{j+1/2}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D70C}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}[d_{t}U^{x}]_{i,j+1/2}^{n}+[D_{x}P]_{i,j+1/2}^{n}=[\unicode[STIX]{x1D70C}g^{x}]_{i,j+1/2}^{n},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D707}K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})+\unicode[STIX]{x1D70C}F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{y}\boldsymbol{U}^{n}])_{i,j+1/2}U_{i+1/2,j}^{y,n}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i+1/2,j}^{-1}\frac{D_{x}U_{i+1,j}^{y,n}-D_{x}U_{i,j}^{y,n}}{h_{i+1/2}}-\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i+1/2,j}^{-1}\frac{d_{x}U_{i+1/2,j+1/2}^{y,n}-d_{x}U_{i+1/2,j-1/2}^{y,n}}{k_{j}}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\unicode[STIX]{x1D70C}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i+1/2,j}^{-1}[d_{t}U^{y}]_{i+1/2,j}^{n}+[D_{y}P]_{i+1/2,j}^{n}=[\unicode[STIX]{x1D70C}g^{y}]_{i+1/2,j}^{n},\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle {\mathcal{P}}_{h}\unicode[STIX]{x1D713}(x,y) & = & \displaystyle \left(\unicode[STIX]{x1D713}_{i-1/2,j-1/2}\left(\frac{x_{i+1/2}-x}{x_{i+1/2}-x_{i-1/2}}\right)+\unicode[STIX]{x1D713}_{i+1/2,j-1/2}\left(\frac{x-x_{i-1/2}}{x_{i+1/2}-x_{i-1/2}}\right)\!\right)\nonumber\\ \displaystyle & & \displaystyle \times \left(\frac{y_{j+1/2}-y}{y_{j+1/2}-y_{j-1/2}}\right)+\left(\unicode[STIX]{x1D713}_{i-1/2,j+1/2}\left(\frac{x_{i+1/2}-x}{x_{i+1/2}-x_{i-1/2}}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\unicode[STIX]{x1D713}_{i+1/2,j+1/2}\left(\frac{x-x_{i-1/2}}{x_{i+1/2}-x_{i-1/2}}\right)\!\right)\left(\frac{y-y_{j-1/2}}{y_{j+1/2}-y_{j-1/2}}\right).\end{eqnarray}$$
For the calculation of the discrete porosity
$\unicode[STIX]{x1D713}$
, we use the following scheme.
$$\begin{eqnarray}[d_{t}\unicode[STIX]{x1D713}]_{i+1/2,j+1/2}^{n}=\frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}(1-\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}\overline{Z}_{i+1/2,j+1/2}^{n-1},\end{eqnarray}$$
where
$\overline{Z}^{n-1}=\max \{0,\min \{Z^{n-1},1\}\}.$
Set the boundary and initial approximations as follows:
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}ll@{}}W_{0,j+1/2}^{x,n}=W_{N_{x},j+1/2}^{x,n}=0, & 0\leqslant j\leqslant N_{y}-1,\\ W_{i+1/2,0}^{y,n}=W_{i+1/2,N_{y}}^{y,n}=0, & 0\leqslant i\leqslant N_{x}-1,\\ U_{0,j+1/2}^{x,n}=U_{N_{x},j+1/2}^{x,n}=0, & 0\leqslant j\leqslant N_{y}-1,\\ U_{i,0}^{x,n}=U_{i,N_{y}}^{x,n}=0, & 0\leqslant i\leqslant N_{x},\\ U_{0,j}^{y,n}=U_{N_{x},j}^{y,n}=0, & 0\leqslant j\leqslant N_{y},\\ U_{i+1/2,0}^{y,n}=U_{i+1/2,N_{y}}^{y,n}=0, & 0\leqslant i\leqslant N_{x}-1,\\ Z_{i+1/2,j+1/2}^{0}=c_{f0,i+1/2,j+1/2}, & 0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1,\\ P_{i+1/2,j+1/2}^{0}=p_{0,i+1/2,j+1/2}, & 0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1,\\ U_{i,j+1/2}^{x,0}=\tilde{u} _{0,i,j+1/2}^{x}, & 0\leqslant i\leqslant N_{x},0\leqslant j\leqslant N_{y},\\ U_{i+1/2,j}^{y,0}=\tilde{u} _{0,i+1/2,j}^{y}, & 0\leqslant i\leqslant N_{x},0\leqslant j\leqslant N_{y},\\ \unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{0}=\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}, & 0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1.\end{array}\right\} & \displaystyle\end{eqnarray}$$
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.
Remark.
In (3.13) and (3.14), we can replace
$\unicode[STIX]{x1D712}(U_{i,j+1/2}^{x,n})$
,
$\unicode[STIX]{x1D712}(U_{i+1/2,j}^{y,n})$
with
$U_{i,j+1/2}^{x,n}$
,
$U_{i+1/2,j}^{y,n}$
. The boundedness of
$\Vert \boldsymbol{U}\Vert _{\infty }^{n}$
can be proved similar to the hypothesis of equation (23) in Li & Rui (Reference Li and Rui2017a
).
4 Error analysis for the discrete scheme
In this section, we give the error estimates of (3.12)–(3.19). Set
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle E_{p}=P-p,\quad E_{\tilde{\boldsymbol{u}}}=(E_{\tilde{\boldsymbol{u}}}^{x},E_{\tilde{\boldsymbol{u}}}^{y})=\boldsymbol{U}-\tilde{\boldsymbol{u}},\\ E_{c_{f}}=Z-c_{f},\quad E_{\boldsymbol{w}}=(E_{\boldsymbol{w}}^{x},E_{\boldsymbol{w}}^{y})=\boldsymbol{W}-\boldsymbol{w},\\ E_{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D713}-\unicode[STIX]{x1D719}.\end{array}\right\}\end{eqnarray}$$
Next we will prove a priori bounds for the discrete solution
$\unicode[STIX]{x1D713}$
which will be used in what follows.
Lemma 8. Assuming that
$0<\unicode[STIX]{x1D719}_{0\ast }\leqslant \unicode[STIX]{x1D719}_{0}\leqslant \unicode[STIX]{x1D719}_{0}^{\ast }<1$
, then the discrete porosity
$\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}$
is bounded, i.e.,
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{0\ast }\leqslant \unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}<1,\quad 0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1,n\leqslant N.\end{eqnarray}$$
It also holds that
$$\begin{eqnarray}0\leqslant [d_{t}\unicode[STIX]{x1D713}]_{i+1/2,j+1/2}^{n}<r_{\unicode[STIX]{x1D719}}^{\ast },\quad 0\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1,n\leqslant N.\end{eqnarray}$$
Proof. The proof is given by induction. It is trivial that
$\unicode[STIX]{x1D719}_{0\ast }\leqslant \unicode[STIX]{x1D713}_{i,j}^{0}<1$
. Suppose that
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{0\ast }\leqslant \unicode[STIX]{x1D713}_{i,j}^{n-1}<1,\quad n\leqslant N;\end{eqnarray}$$
next we prove that
$\unicode[STIX]{x1D713}_{i,j}^{n}$
also does.
For simplicity, set
$\unicode[STIX]{x1D6FD}^{n-1}=r_{\unicode[STIX]{x1D719}}\overline{Z}^{n-1}/(1-\unicode[STIX]{x1D719}_{0})\unicode[STIX]{x0394}t,$
where
$\overline{Z}^{n-1}=\max \{0,\min \{Z^{n-1},1\}\}.$
Then (3.19) can be transformed into the following.
$$\begin{eqnarray}\unicode[STIX]{x1D713}_{i,j}^{n}=\frac{\unicode[STIX]{x1D6FD}_{i,j}^{n-1}}{1+\unicode[STIX]{x1D6FD}_{i,j}^{n-1}}+\frac{\unicode[STIX]{x1D719}_{i,j}^{n-1}}{1+\unicode[STIX]{x1D6FD}_{i,j}^{n-1}}.\end{eqnarray}$$
By (4.5), we can easily obtain that
$\unicode[STIX]{x1D719}_{i,j}^{n}<1$
. Equation (3.19) can also be calculated as
$$\begin{eqnarray}\unicode[STIX]{x1D719}_{i,j}^{n}-\unicode[STIX]{x1D719}_{i,j}^{n-1}=\unicode[STIX]{x1D6FD}_{i,j}^{n-1}(1-\unicode[STIX]{x1D719}_{i,j}^{n}),\end{eqnarray}$$
and thus we have that
$\unicode[STIX]{x1D719}_{i,j}^{n}>\unicode[STIX]{x1D719}_{i,j}^{n-1}$
. This completes the proof.
Lemma 9. The approximate errors of discrete porosity satisfy
$$\begin{eqnarray}\displaystyle & \displaystyle \Vert E_{\unicode[STIX]{x1D719}}^{m}\Vert _{l^{2},M}^{2}\leqslant C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert E_{c_{f}}^{n-1}\Vert _{l^{2},M}^{2})+C\unicode[STIX]{x0394}t^{2}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}\leqslant C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert E_{c_{f}}^{n-1}\Vert _{l^{2},M}^{2})+C\unicode[STIX]{x0394}t^{2},\quad m\leqslant N, & \displaystyle\end{eqnarray}$$
where the positive constant
$C$
is independent of
$h,~k$
and
$\unicode[STIX]{x0394}t$
.
Proof. Subtracting (3.2) from (3.19), we can obtain
$$\begin{eqnarray}\displaystyle d_{t}E_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n} & = & \displaystyle \frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}(1-\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}\overline{Z}_{i+1/2,j+1/2}^{n-1}\nonumber\\ \displaystyle & & \displaystyle -\,\frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}(1-\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n})}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}c_{f,i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle +\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n}}{\unicode[STIX]{x2202}t}-d_{t}\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & = & \displaystyle \frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}(1-\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}(\overline{Z}_{i+1/2,j+1/2}^{n-1}-c_{f,i+1/2,j+1/2}^{n-1})\nonumber\\ \displaystyle & & \displaystyle +\,\frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}(1-\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}(c_{f,i+1/2,j+1/2}^{n-1}-c_{f,i+1/2,j+1/2}^{n})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{r_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}c_{f,i+1/2,j+1/2}^{n}}{1-\unicode[STIX]{x1D719}_{0,i+1/2,j+1/2}}E_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}+R_{1,i+1/2,j+1/2}^{n},\end{eqnarray}$$
where
$R_{1}^{n}=\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{n}/\unicode[STIX]{x2202}t-d_{t}\unicode[STIX]{x1D719}^{n}.$
Multiplying (4.9) by
$E_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}h_{i+1/2}k_{j+1/2}$
and making a summation on
$i,j$
for
$0\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$
, we have that
$$\begin{eqnarray}\displaystyle (d_{t}E_{\unicode[STIX]{x1D719}}^{n},E_{\unicode[STIX]{x1D719}}^{n})_{l^{2},M} & = & \displaystyle \left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(\overline{Z}^{n-1}-c_{f}^{n-1}),E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle +\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(c_{f}^{n-1}-c_{f}^{n}),E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle -\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}c_{f}^{n}}{1-\unicode[STIX]{x1D719}_{0}}E_{\unicode[STIX]{x1D719}}^{n},E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}+(R_{1}^{n},E_{\unicode[STIX]{x1D719}}^{n})_{l^{2},M}.\end{eqnarray}$$
The term on the left-hand side of (4.10) can be transformed into
$$\begin{eqnarray}(d_{t}E_{\unicode[STIX]{x1D719}}^{n},E_{\unicode[STIX]{x1D719}}^{n})_{l^{2},M}=\frac{\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}-\Vert E_{\unicode[STIX]{x1D719}}^{n-1}\Vert _{l^{2},M}^{2}}{2\unicode[STIX]{x0394}t}+\frac{\unicode[STIX]{x0394}t}{2}\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}.\end{eqnarray}$$
Taking notice of Lemma 8, the first term on the right-hand side of (4.10) can be estimated by
$$\begin{eqnarray}\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(\overline{Z}^{n-1}-c_{f}^{n-1}),E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\leqslant C(\Vert E_{c_{f}}^{n-1}\Vert _{l^{2},M}^{2}+\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}),\end{eqnarray}$$
where we used the fact that
$|\overline{Z}^{n-1}-c_{f}^{n-1}|\leqslant |Z^{n-1}-c_{f}^{n-1}|$
.
The second term on the right-hand side of (4.10) can be bounded by
$$\begin{eqnarray}\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(c_{f}^{n-1}-c_{f}^{n}),E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\leqslant C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C\Vert \frac{\unicode[STIX]{x2202}c_{f}}{\unicode[STIX]{x2202}t}\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}^{2}\unicode[STIX]{x0394}t^{2}.\end{eqnarray}$$
Combining (4.10) with (4.11)–(4.13), multiplying by
$2\unicode[STIX]{x0394}t$
, and summing for
$n$
from 1 to
$m$
,
$m\leqslant N$
, we have
$$\begin{eqnarray}\Vert E_{\unicode[STIX]{x1D719}}^{M}\Vert _{l^{2},M}^{2}\leqslant \Vert E_{\unicode[STIX]{x1D719}}^{0}\Vert _{l^{2},M}^{2}+C\mathop{\sum }_{n=1}^{M}\unicode[STIX]{x0394}t\Vert E_{c_{f}}^{A,n-1}\Vert _{l^{2},M}^{2}+C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C\unicode[STIX]{x0394}t^{2}.\end{eqnarray}$$
Recalling the initial condition on
$\unicode[STIX]{x1D713}^{0}$
gives (4.7).
On the other hand, multiplying (4.9) by
$d_{t}E_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}h_{i+1/2}k_{j+1/2}$
and making a summation on
$i,j$
for
$0\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$
, we have that
$$\begin{eqnarray}\displaystyle \Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2} & = & \displaystyle \left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(\overline{Z}^{n-1}-c_{f}^{n-1}),d_{t}E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle +\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}(1-\unicode[STIX]{x1D713}^{n})}{1-\unicode[STIX]{x1D719}_{0}}(c_{f}^{n-1}-c_{f}^{n}),d_{t}E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle -\left(\frac{r_{\unicode[STIX]{x1D719}}^{n}c_{f}^{n}}{1-\unicode[STIX]{x1D719}_{0}}E_{\unicode[STIX]{x1D719}}^{n},d_{t}E_{\unicode[STIX]{x1D719}}^{n}\right)_{l^{2},M}+(R_{1}^{n},d_{t}E_{\unicode[STIX]{x1D719}}^{n})_{l^{2},M}.\end{eqnarray}$$
Similarly we can obtain that
$$\begin{eqnarray}\unicode[STIX]{x0394}t\mathop{\sum }_{n=1}^{m}\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}\leqslant C\unicode[STIX]{x0394}t\mathop{\sum }_{n=1}^{m}\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C\unicode[STIX]{x0394}t\mathop{\sum }_{n=1}^{m}\Vert E_{c_{f}}^{A,n-1}\Vert _{l^{2},M}^{2}+C\unicode[STIX]{x0394}t^{2}.\end{eqnarray}$$
This completes the proof.
Lemma 10. The approximate errors of discrete pressure and velocity satisfy
$$\begin{eqnarray}\displaystyle \Vert E_{p}^{m}\Vert _{l^{2},M}^{2}+\Vert E_{\tilde{\boldsymbol{u}}}^{m}\Vert _{l^{2}}^{2}+\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert \unicode[STIX]{x1D63F}E_{\tilde{\boldsymbol{u}}}^{n}\Vert ^{2} & {\leqslant} & \displaystyle C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2})\nonumber\\ \displaystyle & & \displaystyle +\,O(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}),\quad m\leqslant N,\end{eqnarray}$$
where the positive constant
$C$
is independent of
$h,~k$
and
$\unicode[STIX]{x0394}t$
.
Proof. Define
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D6FF}_{p,i+1/2,j+1/2}^{\tilde{\boldsymbol{u}}}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{h_{i+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}p_{i+1/2,j+1/2}^{n}}{\unicode[STIX]{x2202}x^{2}}+\frac{k_{j+1/2}^{2}}{8}\frac{\unicode[STIX]{x2202}^{2}p_{i+1/2,j+1/2}^{n}}{\unicode[STIX]{x2202}y^{2}}+\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{4\sqrt{K(\unicode[STIX]{x1D719})}}\right)_{i+1/2,j+1/2}\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,\frac{[\tilde{u} ^{x}\tilde{u} ^{y}]_{i+1/2,j+1/2}}{|\tilde{\boldsymbol{u}}_{i+1/2,j+1/2}|}\left(\frac{\unicode[STIX]{x2202}\tilde{u} _{i+1/2,j+1/2}^{y}}{\unicode[STIX]{x2202}x}h_{i+1/2}^{2}+\frac{\unicode[STIX]{x2202}\tilde{u} _{i+1/2,j+1/2}^{x}}{\unicode[STIX]{x2202}y}k_{j+1/2}^{2}\right).\end{eqnarray}$$
Define the residuals as follows
$$\begin{eqnarray}\displaystyle & r_{i,j+1/2}^{x}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}r_{i,j+1/2}^{x,x}+\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}r_{i,j+1/2}^{x,y},\quad 1\leqslant i\leqslant N_{x}-1,0\leqslant j\leqslant N_{y}-1, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & r_{i+1/2,j}^{y}=\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}r_{i+1/2,j}^{y,y}+\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}r_{i+1/2,j}^{y,x},\quad 0\leqslant i\leqslant N_{x}-1,1\leqslant j\leqslant N_{y}-1. & \displaystyle\end{eqnarray}$$
From Lemmas 4 and 5, we obtain that for
$i=1,\ldots ,N_{x}-1,~j=0,\ldots ,N_{y}-1$
,
$$\begin{eqnarray}\displaystyle & & \displaystyle \left.\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}d_{t}\tilde{u} _{i,j+1/2}^{x,n}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}\right|_{i,j+1/2}^{n}\left.\frac{d_{x}\tilde{u} _{i+1/2,j+1/2}^{x,n}-d_{x}\tilde{u} _{i-1/2,j+1/2}^{x,n}}{h_{i}}-\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}\right|_{i,j+1/2}^{n}\frac{D_{y}\tilde{u} _{i,j+1}^{x,n}-D_{y}\tilde{u} _{i,j}^{x,n}}{k_{j+1/2}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\frac{\unicode[STIX]{x1D707}}{K}\tilde{u} _{i,j+1/2}^{x,n}+\left.\frac{\unicode[STIX]{x1D70C}F}{\sqrt{K}}\right|_{i,j+1/2}^{n}[S_{x}\tilde{\boldsymbol{u}}]\tilde{u} _{i,j+1/2}^{x,n}+[D_{x}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70C}g_{i,j+1/2}^{x,n}-\left.\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}\right|_{i,j+1/2}^{n}\frac{\unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{x,x,n}-\unicode[STIX]{x1D716}_{i-1/2,j+1/2}^{x,x}}{h_{i}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\left.\frac{\unicode[STIX]{x1D707}}{\unicode[STIX]{x1D719}}\right|_{i,j+1/2}^{n}\frac{\unicode[STIX]{x1D716}_{i,j+1}^{x,y,n}-\unicode[STIX]{x1D716}_{i,j}^{x,y}}{k_{j+1/2}}-r_{i,j+1/2}^{x,n}+R_{i,j+1/2}^{x,n},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle R_{i,j+1/2}^{x,n} & = & \displaystyle [D_{x}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i,j+1/2}^{n}-\left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}\right|_{i,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \,+\left(\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\right)_{i,j+1/2}^{n}([S_{x}\tilde{\boldsymbol{u}}]_{i,j+1/2}^{n}-|\tilde{\boldsymbol{u}}|_{i,j+1/2}^{n})\tilde{u} _{i,j+1/2}^{x,n}\right)\nonumber\\ \displaystyle & & \displaystyle \,+\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}d_{t}\tilde{u} _{i,j+1/2}^{x,n}-\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}u_{i,j+1/2}^{x,n}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D707}}{K}(\tilde{u} ^{x,n}-u^{x,n})_{i,j+1/2}\nonumber\\ \displaystyle & & \displaystyle \,+\frac{\unicode[STIX]{x1D70C}F}{\sqrt{K}}|_{i,j+1/2}^{n}(|\tilde{\boldsymbol{u}}|\tilde{u} _{i,j+1/2}^{x,n}-|\boldsymbol{u}|u_{i,j+1/2}^{x,n}).\end{eqnarray}$$
Similarly to the derivation of Lemma 4.2 in Rui et al. (Reference Rui, Zhao and Pan2015), we can obtain that when
$\unicode[STIX]{x1D70C}F/\sqrt{K}\in C^{1}(\unicode[STIX]{x1D6FA})$
,
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\right)_{i,j+1/2}^{n}([S_{x}\tilde{\boldsymbol{u}}]_{i,j+1/2}^{n}-|\tilde{\boldsymbol{u}}|_{i,j+1/2}^{n})\tilde{u} _{i,j+1/2}^{x,n}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \quad =\frac{1}{4h_{i}}\left\{\!\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\frac{\tilde{u} ^{x}\tilde{u} ^{y}}{|\tilde{\boldsymbol{u}}|}\frac{\unicode[STIX]{x2202}\tilde{u} ^{y}}{\unicode[STIX]{x2202}x}\right)_{i+1/2,j+1/2}^{n}h_{i+1/2}^{2}-\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\frac{\tilde{u} ^{x}\tilde{u} ^{y}}{|\tilde{\boldsymbol{u}}|}\frac{\unicode[STIX]{x2202}\tilde{u} ^{y}}{\unicode[STIX]{x2202}x}\right)_{i-1/2,j+1/2}^{n}h_{i-1/2}^{2}\right\}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,O(h^{2}+k^{2}).\end{eqnarray}$$
Thus we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \left|[D_{x}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i,j+1/2}^{n}-\left.\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}x}\right|_{i,j+1/2}^{n}+\left(\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\right)_{i,j+1/2}^{n}([S_{x}\tilde{\boldsymbol{u}}]_{i,j+1/2}^{n}-|\tilde{\boldsymbol{u}}_{i,j+1/2}^{n}|)\tilde{u} _{i,j+1/2}^{x,n}\right)\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C(h^{2}+k^{2}).\end{eqnarray}$$
Recalling Lemma 3, (2.2) can be transformed into
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}[d_{t}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i+1/2,j+1/2}^{n}+[d_{t}\unicode[STIX]{x1D719}]_{i+1/2,j+1/2}^{n}+[d_{x}\tilde{u} ^{x}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad +\,[d_{y}\tilde{u} ^{y}]_{i+1/2,j+1/2}^{n}=f_{i+1/2,j+1/2}^{n}+R_{i+1/2,j+1/2}^{p,n},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle R_{i+1/2,j+1/2}^{p,n} & = & \displaystyle \unicode[STIX]{x1D700}[d_{t}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i+1/2,j+1/2}^{n}-\left.\unicode[STIX]{x1D700}\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}t}\right|_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle +\,[d_{t}\unicode[STIX]{x1D719}]_{i+1/2,j+1/2}^{n}-\left.\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t}\right|_{i+1/2,j+1/2}^{n}+r_{i+1/2,j+1/2}^{\boldsymbol{u},n}.\end{eqnarray}$$
Subtracting (4.21) from (3.16), we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}[d_{t}E_{\tilde{\boldsymbol{u}}}^{x}]_{i,j+1/2}^{n}+\unicode[STIX]{x1D707}[K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})E_{\tilde{\boldsymbol{u}}}^{x}]_{i,j+1/2}^{n}+[D_{x}(E_{p}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}])_{i,j+1/2}U_{i,j+1/2}^{x,n}-\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}[S_{x}\tilde{\boldsymbol{u}}]\tilde{u} ^{x}\right)_{i,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{d_{x}E_{\tilde{\boldsymbol{u}},i+1/2,j+1/2}^{x,n}-d_{x}E_{\tilde{\boldsymbol{u}},i-1/2,j+1/2}^{x,n}}{h_{i}}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{D_{y}E_{\tilde{\boldsymbol{u}},i,j+1}^{x,n}-D_{y}E_{\tilde{\boldsymbol{u}},i,j}^{x,n}}{k_{j+1/2}}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D707}\left([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}-[\unicode[STIX]{x1D719}^{n}]_{i,j+1/2}^{-1}\right)\left(\frac{\unicode[STIX]{x2202}^{2}u_{i,j+1/2}^{x,n}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}u_{i,j+1/2}^{x,n}}{\unicode[STIX]{x2202}y^{2}}+r_{i,j+1/2}^{x,x,n}+r_{i,j+1/2}^{x,y,n}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{\unicode[STIX]{x1D716}_{i+1/2,j+1/2}^{x,x,n}-\unicode[STIX]{x1D716}_{i-1/2,j+1/2}^{x,x,n}}{h_{i}}+\unicode[STIX]{x1D707}[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\frac{\unicode[STIX]{x1D716}_{i,j+1}^{x,y,n}-\unicode[STIX]{x1D716}_{i,j}^{x,y,n}}{k_{j+1/2}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}(K^{-1}(\unicode[STIX]{x1D719})-K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}))_{i,j+1/2}^{n}\tilde{u} _{i,j+1/2}^{x,n}+\unicode[STIX]{x1D70C}\left([\unicode[STIX]{x1D719}^{n}]_{i,j+1/2}^{-1}-[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]_{i,j+1/2}^{-1}\right)\nonumber\\ \displaystyle & & \displaystyle \qquad \times \,[d_{t}\tilde{u} ^{x}]_{i+1/2,j}^{n}+\,r_{i,j+1/2}^{x,n}-R_{i,j+1/2}^{x,n}.\end{eqnarray}$$
Multiplying (4.28) by
$E_{\tilde{\boldsymbol{u}},i,j+1/2}^{x,n}h_{i}k_{j+1/2}$
and making a summation on
$i,~j$
for
$1\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$
, we can obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}(({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},T,M}+\unicode[STIX]{x1D707}\Vert \sqrt{K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T,M}^{2}-(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n},d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}]U^{x,n}-\frac{F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{x}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\left(d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n},d_{x}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n})\right)_{l^{2},M}+\unicode[STIX]{x1D707}\left(D_{y}E_{\tilde{\boldsymbol{u}}}^{x,n},D_{y}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n})\right)_{l^{2},T_{y}}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D707}\left(\left([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}-[\unicode[STIX]{x1D719}^{n}]^{-1}\right)\left(\frac{\unicode[STIX]{x2202}^{2}u^{x,n}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}u^{x,n}}{\unicode[STIX]{x2202}y^{2}}+r^{x,x,n}+r^{x,y,n}\right)\!,E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{x,x,n},d_{x}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n}))_{l^{2},M}-\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{x,y,n},D_{y}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n}))_{l^{2},T_{y}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\left(\left(K^{-1}(\unicode[STIX]{x1D719}^{n})-K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})\right)\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(\left([\unicode[STIX]{x1D719}^{n}]^{-1}-[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}\right)d_{t}\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}+\left((r^{x,n}-R^{x,n}),E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\!.\end{eqnarray}$$
The fifth term on the left-hand side of (4.29) can be estimated by
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D707}\left(d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n},d_{x}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n})\right)_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \quad \geqslant \unicode[STIX]{x1D707}\Vert d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},M}^{2}-C\Vert d_{x}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}\Vert _{\infty }\Vert E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T,M},\end{eqnarray}$$
where taking notice of equation (3.19) and the boundedness of
$\boldsymbol{W}^{n-1}$
, we can obtain that
$\Vert d_{x}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}\Vert _{\infty }\leqslant ~C$
.
The last term on the left-hand side of (4.29) can be estimated similar to (4.30). Then noting Lemma 8, we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},T,M}-(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n},d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}]U^{x,n}-\frac{F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{x}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\Vert d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},M}^{2}+\unicode[STIX]{x1D707}\Vert D_{y}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T_{y}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \unicode[STIX]{x1D707}\left(\left([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}-[\unicode[STIX]{x1D719}^{n}]^{-1}\right)\left(\frac{\unicode[STIX]{x2202}^{2}u^{x,n}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}u^{x,n}}{\unicode[STIX]{x2202}y^{2}}+r^{x,x,n}+r^{x,y,n}\right)\!,E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{x,x,n},d_{x}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n}))_{l^{2},M}-\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{x,y,n},D_{y}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n}))_{l^{2},T_{y}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\left(\left(K^{-1}(\unicode[STIX]{x1D719}^{n})-K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})\right)\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(\left([\unicode[STIX]{x1D719}^{n}]^{-1}-[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}\right)d_{t}\tilde{u} ^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x1D707}}{4}\Vert d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},M}^{2}+\frac{\unicode[STIX]{x1D707}}{4}\Vert D_{y}E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T_{y}}^{2}+C\Vert E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T,M}^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\left((r^{x,n}-R^{x,n}),E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}+C(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
Similarly in the
$y$
direction, we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n})_{l^{2},M,T}-(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n},d_{y}E_{\tilde{\boldsymbol{u}}}^{y,n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{y}\boldsymbol{U}^{n}]U^{y,n}-\frac{F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{y}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\Vert d_{y}E_{\tilde{\boldsymbol{u}}}^{y,n}\Vert _{l^{2},M}^{2}+\unicode[STIX]{x1D707}\Vert D_{x}E_{\tilde{\boldsymbol{u}}}^{y,n}\Vert _{l^{2},T_{x}}^{2}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D707}\left(\left([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}-[\unicode[STIX]{x1D719}^{n}]^{-1}\right)\left(\frac{\unicode[STIX]{x2202}^{2}u^{y,n}}{\unicode[STIX]{x2202}x^{2}}+\frac{\unicode[STIX]{x2202}^{2}u^{y,n}}{\unicode[STIX]{x2202}y^{2}}+r^{y,x,n}+r^{y,y,n}\right)\!,E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \qquad -\,\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{y,y,n},d_{y}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{y,n}))_{l^{2},M}-\unicode[STIX]{x1D707}(\unicode[STIX]{x1D716}^{y,x,n},D_{x}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}E_{\tilde{\boldsymbol{u}}}^{y,n}))_{l^{2},T_{x}}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\left(\left(K^{-1}(\unicode[STIX]{x1D719}^{n})-K^{-1}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})\right)\tilde{u} ^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(\left([\unicode[STIX]{x1D719}^{n}]^{-1}-[{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}\right)d_{t}\tilde{u} ^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x1D707}}{4}\Vert d_{y}E_{\tilde{\boldsymbol{u}}}^{y,n}\Vert _{l^{2},M}^{2}+\frac{\unicode[STIX]{x1D707}}{4}\Vert D_{x}E_{\tilde{\boldsymbol{u}}}^{y,n}\Vert _{l^{2},T_{x}}^{2}+C\Vert E_{\tilde{\boldsymbol{u}}}^{y,n}\Vert _{l^{2},M,T}^{2}\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \qquad +\left((r^{y,n}-R^{y,n}),E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}+C(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}),\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle & & \displaystyle R_{i+1/2,j}^{y,n}\nonumber\\ \displaystyle & & \displaystyle \quad =[D_{y}(p-\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i+1/2,j}^{n}-\frac{\unicode[STIX]{x2202}p}{\unicode[STIX]{x2202}y}|_{i+1/2,j}^{n}+\left(\!\!\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719})}{\sqrt{K(\unicode[STIX]{x1D719})}}\right)_{i+1/2,j}^{n}\!\!([S_{y}\tilde{\boldsymbol{u}}]_{i+1/2,j}^{n}-|\tilde{\boldsymbol{u}}_{i+1/2,j}^{n}|)\tilde{u} _{i+1/2,j}^{y,n}\!\right)\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}d_{t}\tilde{u} _{i+1/2,j}^{y,n}-\frac{\unicode[STIX]{x1D70C}}{\unicode[STIX]{x1D719}}\frac{\unicode[STIX]{x2202}u_{i+1/2,j}^{y,n}}{\unicode[STIX]{x2202}t}+\frac{\unicode[STIX]{x1D707}}{K}(\tilde{u} ^{y,n}-u^{y,n})_{i+1/2,j}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\frac{\unicode[STIX]{x1D70C}F}{\sqrt{K}}|_{i+1/2,j}^{n}(|\tilde{\boldsymbol{u}}|\tilde{u} _{i+1/2,j}^{y,n}-|\boldsymbol{u}|u_{i+1/2,j}^{y,n}).\end{eqnarray}$$
Using the similar estimation in Pan & Rui (Reference Pan and Rui2012), we have
$$\begin{eqnarray}\left|\frac{\unicode[STIX]{x1D70C}F}{\sqrt{K}}|^{n}(|\tilde{\boldsymbol{u}}|\tilde{\boldsymbol{u}}^{n}-|\boldsymbol{u}|\boldsymbol{u}^{n})\right|\leqslant C|\tilde{\boldsymbol{u}}^{n}-\boldsymbol{u}^{n}|\leqslant C(h^{2}+k^{2}).\end{eqnarray}$$
Recalling the estimates of (4.25), we can easily obtain
$$\begin{eqnarray}\Vert \text{}\text{R}^{n}\Vert _{l^{2}}\leqslant C(h^{2}+k^{2}+\unicode[STIX]{x0394}t),\end{eqnarray}$$
where
$\boldsymbol{R}^{n}=(R^{x,n},R^{y,n})$
.
Subtracting (4.26) from (3.15), we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}[d_{t}(E_{p}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}})]_{i+1/2,j+1/2}^{n}+[d_{t}E_{\unicode[STIX]{x1D719}}]_{i+1/2,j+1/2}^{n}+[d_{x}E_{\tilde{\boldsymbol{u}}}^{x}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad +\,[d_{y}E_{\tilde{\boldsymbol{u}}}^{y}]_{i+1/2,j+1/2}^{n}=-R_{i+1/2,j+1/2}^{p,n}.\end{eqnarray}$$
Multiplying (4.41) by
$[E_{p}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}}}]_{i+1/2,j+1/2}^{n}h_{i+1/2}k_{j+1/2}$
and making a summation on
$i,~j$
for
$0\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$
, we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}\left(d_{t}(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}),E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}\right)_{l^{2},M}+(d_{x}E_{\tilde{\boldsymbol{u}}}^{x,n},E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}+(d_{y}E_{\tilde{\boldsymbol{u}}}^{y,n},E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \quad =-(d_{t}E_{\unicode[STIX]{x1D719}}^{n},E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}-(R^{p,n},E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}.\end{eqnarray}$$
Combining (4.42) with (4.31) and (4.32) and using the Cauchy–Schwartz inequality, we can obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D700}(d_{t}(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}),E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}+\unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n})_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\!\left(\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}]\right)U^{x,n}-\left(\!\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{x}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{x,n}\right),E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\!\left(\!\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{y}\boldsymbol{U}^{n}]\right)U^{y,n}-\left(\!\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{y}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{y,n}\right)\!,E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D707}\Vert \unicode[STIX]{x1D63F}E_{\tilde{\boldsymbol{u}}}^{n}\Vert ^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\Vert \unicode[STIX]{x1D719}^{n}-{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}\Vert _{l^{2}}^{2}+C\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+\frac{\unicode[STIX]{x1D707}}{2}\Vert \unicode[STIX]{x1D63F}E_{\tilde{\boldsymbol{u}}}^{n}\Vert ^{2}+C\Vert E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}\Vert _{l^{2},M}^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,C(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
The first term on the left-hand side of (4.43) can be estimated by
$$\begin{eqnarray}\unicode[STIX]{x1D700}(d_{t}(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}),E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})_{l^{2},M}=\frac{\unicode[STIX]{x1D700}}{2}d_{t}\Vert E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n}\Vert _{l^{2},M}^{2}+\frac{\unicode[STIX]{x1D700}\unicode[STIX]{x0394}t}{2}\Vert d_{t}(E_{p}^{n}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},n})\Vert _{l^{2},M}^{2}.\end{eqnarray}$$
The second and third terms on the left-hand side of (4.43) can be transformed into the following.
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{x,n},E_{\tilde{\boldsymbol{u}}}^{x,n})_{l^{2},T,M}+\unicode[STIX]{x1D70C}([{\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}]^{-1}d_{t}E_{\tilde{\boldsymbol{u}}}^{y,n},E_{\tilde{\boldsymbol{u}}}^{y,n})_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{\unicode[STIX]{x1D70C}}{2}d_{t}\Vert \sqrt{({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}}E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x0394}t}{2}\Vert \sqrt{({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}}d_{t}E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,\frac{\unicode[STIX]{x1D70C}}{2}(d_{t}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}E_{\tilde{\boldsymbol{u}}}^{x,n-1},E_{\tilde{\boldsymbol{u}}}^{x,n-1})_{l^{2},T,M}-\frac{\unicode[STIX]{x1D70C}}{2}(d_{t}({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})^{-1}E_{\tilde{\boldsymbol{u}}}^{y,n-1},E_{\tilde{\boldsymbol{u}}}^{y,n-1})_{l^{2},M,T}.\end{eqnarray}$$
The fourth and fifth terms on the left-hand side of (4.43) can be estimated by
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D70C}\left((F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{x}\boldsymbol{U}^{n}])U^{x,n}-\left(\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{x}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{x,n}\right),E_{\tilde{\boldsymbol{u}}}^{x,n}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left((F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}[S_{y}\boldsymbol{U}^{n}])U^{y,n}-\left(\!\frac{\unicode[STIX]{x1D70C}F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}[S_{y}\tilde{\boldsymbol{u}}^{n}]\tilde{u} ^{y,n}\right),E_{\tilde{\boldsymbol{u}}}^{y,n}\right)_{l^{2},M,T}\nonumber\\ \displaystyle & & \displaystyle \quad =\unicode[STIX]{x1D70C}(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}\left([\boldsymbol{S}\boldsymbol{U}^{n}]\boldsymbol{U}^{n}-[\boldsymbol{S}\tilde{\boldsymbol{u}}^{n}]\tilde{\boldsymbol{u}}^{x,n}\right),E_{\tilde{\boldsymbol{u}}}^{n})_{l^{2}}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}\left(\left(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}-\frac{F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}}\right)[\boldsymbol{S}\tilde{\boldsymbol{u}}^{n}]\tilde{\boldsymbol{u}}^{n},E_{\tilde{\boldsymbol{u}}}^{n}\right)_{l^{2}}\!.\end{eqnarray}$$
Similar to the estimate of Lemma 4.5 in Rui et al. (Reference Rui, Zhao and Pan2015), we have
$$\begin{eqnarray}\unicode[STIX]{x1D70C}(F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}\left([\boldsymbol{S}\boldsymbol{U}^{n}]\boldsymbol{U}^{n}-[\boldsymbol{S}\tilde{\boldsymbol{u}}^{n}]\tilde{\boldsymbol{u}}^{n}\right)\!,E_{\tilde{\boldsymbol{u}}}^{n})_{l^{2}}\geqslant 0.\end{eqnarray}$$
Using the uniform Lipschitz continuity of
$F(\cdot )/\sqrt{K(\cdot )}$
and Cauchy–Schwartz inequality, the last term on the right-hand side of (4.46) can be estimated by
$$\begin{eqnarray}\displaystyle & & \displaystyle \left|\unicode[STIX]{x1D70C}\left((F({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})(\sqrt{K({\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n})})^{-1}-\frac{F(\unicode[STIX]{x1D719}^{n})}{\sqrt{K(\unicode[STIX]{x1D719}^{n})}})[\boldsymbol{S}\tilde{\boldsymbol{u}}^{n}]\tilde{\boldsymbol{u}}^{n},E_{\tilde{\boldsymbol{u}}}^{n}\right)_{l^{2}}\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+C\Vert {\mathcal{P}}_{h}\unicode[STIX]{x1D713}^{n}-{\mathcal{P}}_{h}\unicode[STIX]{x1D719}^{n}\Vert _{l^{2}}^{2}+C\Vert {\mathcal{P}}_{h}\unicode[STIX]{x1D719}^{n}-\unicode[STIX]{x1D719}^{n}\Vert _{l^{2}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C(h^{4}+k^{4}).\end{eqnarray}$$
Combining (4.43) with (4.44)–(4.48), multiplying (4.43) by
$\unicode[STIX]{x0394}t$
, summing for
$n$
from 1 to
$m$
,
$m\leqslant N$
and applying Gronwall’s lemma give
$$\begin{eqnarray}\displaystyle & & \displaystyle \Vert E_{p}^{m}+\unicode[STIX]{x1D6FF}_{p}^{\tilde{\boldsymbol{u}},m}\Vert _{l^{2},M}^{2}+\Vert E_{\tilde{\boldsymbol{u}}}^{m}\Vert _{l^{2}}^{2}+\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert \unicode[STIX]{x1D63F}E_{\tilde{\boldsymbol{u}}}^{n}\Vert ^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2})+O(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
This completes the proof.
Lemma 11. The approximate errors of discrete concentration and flux variable satisfy
$$\begin{eqnarray}\displaystyle & & \displaystyle \Vert E_{c_{f}}^{m}+\unicode[STIX]{x1D6FF}_{c_{f}}^{m}\Vert _{l^{2},M}^{2}+C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert E_{\boldsymbol{w}}^{n}\Vert _{l^{2}}^{2}\leqslant C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2})\nonumber\\ \displaystyle & & \displaystyle \quad +\,C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+O(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}),\quad m\leqslant N,\end{eqnarray}$$
where the positive constant
$C$
is independent of
$h,~k$
and
$\unicode[STIX]{x0394}t$
.
Proof. From (3.1), we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n}[d_{t}(c_{f}-\unicode[STIX]{x1D6FF}_{c_{f}})]_{i+1/2,j+1/2}^{n}+c_{f,i+1/2,j+1/2}^{n}[d_{t}\unicode[STIX]{x1D719}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,[d_{x}w^{x}]_{i+1/2,j+1/2}^{n}+[d_{y}w^{y}]_{i+1/2,j+1/2}^{n}+\frac{k_{c}k_{s}}{k_{c}+k_{s}}a_{v}(\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n})(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{i+1/2,j+1/2}\nonumber\\ \displaystyle & & \displaystyle \qquad -\,f_{P,i+1/2,j+1/2}^{n}(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{i+1/2,j+1/2}\nonumber\\ \displaystyle & & \displaystyle \quad =[f_{I}c_{I}]_{i+1/2,j+1/2}^{n}+\mathop{\sum }_{l=1}^{4}R_{l,i+1/2,j+1/2}^{n},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle R_{1,i+1/2,j+1/2}^{n} & = & \displaystyle \left(\unicode[STIX]{x1D719}^{n}d_{t}c_{f}^{n}-\unicode[STIX]{x1D719}^{n}\frac{\unicode[STIX]{x2202}c_{f}^{n}}{\unicode[STIX]{x2202}t}+c_{f}^{n}d_{t}\unicode[STIX]{x1D719}^{n}-c_{f}^{n}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{n}}{\unicode[STIX]{x2202}t}\right)_{i+1/2,j+1/2}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle C\left(\left\Vert \frac{\unicode[STIX]{x2202}^{2}c_{f}}{\unicode[STIX]{x2202}t^{2}}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}+\left\Vert \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}t^{2}}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}\right)\unicode[STIX]{x0394}t,\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle R_{2,i+1/2,j+1/2}^{n} & = & \displaystyle \left(d_{x}w^{x,n}-\frac{\unicode[STIX]{x2202}w^{x,n}}{\unicode[STIX]{x2202}x}+d_{y}w^{y,n}-\frac{\unicode[STIX]{x2202}w^{y,n}}{\unicode[STIX]{x2202}y}\right)_{i+1/2,j+1/2}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle C\left(\left\Vert \frac{\unicode[STIX]{x2202}^{3}w^{x}}{\unicode[STIX]{x2202}x^{3}}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}h^{2}+\left\Vert \frac{\unicode[STIX]{x2202}^{3}w^{y}}{\unicode[STIX]{x2202}y^{3}}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}k^{2}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle R_{3,i+1/2,j+1/2}^{n} & = & \displaystyle -\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n}[d_{t}\unicode[STIX]{x1D6FF}_{c_{f}}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle C\left(\left\Vert \frac{\unicode[STIX]{x2202}^{3}c_{f}}{\unicode[STIX]{x2202}x^{2}\unicode[STIX]{x2202}t}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}h^{2}+\left\Vert \frac{\unicode[STIX]{x2202}^{3}c_{f}}{\unicode[STIX]{x2202}y^{2}\unicode[STIX]{x2202}t}\right\Vert _{L^{\infty }(J;L^{\infty }(\unicode[STIX]{x1D6FA}))}k^{2}\right),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle R_{4,i+1/2,j+1/2}^{n} & = & \displaystyle -\frac{k_{c}k_{s}}{k_{c}+k_{s}}a_{v}(\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n})\unicode[STIX]{x1D6FF}_{c_{f},i+1/2,j+1/2}^{n}+[f_{P}\unicode[STIX]{x1D6FF}_{c_{f}}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle C(h^{2}+k^{2}).\end{eqnarray}$$
Then, subtracting (4.51) from (3.12), we can obtain that
$$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n}[d_{t}(E_{c_{f}}+\unicode[STIX]{x1D6FF}_{c_{f}})]_{i+1/2,j+1/2}^{n}+Z_{i+1/2,j+1/2}^{n}[d_{t}\unicode[STIX]{x1D713}]_{i+1/2,j+1/2}^{n}-c_{f,i+1/2,j+1/2}^{n}[d_{t}\unicode[STIX]{x1D719}]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,[d_{x}E_{\boldsymbol{w}}^{x}]_{i+1/2,j+1/2}^{n}+[d_{y}E_{\boldsymbol{w}}^{y}]_{i+1/2,j+1/2}^{n}-[f_{P}(E_{c_{f}}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})]_{i+1/2,j+1/2}^{n}\nonumber\\ \displaystyle & & \displaystyle \quad =-\frac{k_{c}k_{s}}{k_{c}+k_{s}}(a_{v}(\unicode[STIX]{x1D713}_{i+1/2,j+1/2}^{n})Z_{i+1/2,j+1/2}^{n}-a_{v}(\unicode[STIX]{x1D719}_{i+1/2,j+1/2}^{n})(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{i+1/2,j+1/2})\nonumber\\ \displaystyle & & \displaystyle \qquad -\,E_{\unicode[STIX]{x1D719},i+1/2,j+1/2}^{n}[d_{t}(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n})]_{i+1/2,j+1/2}^{n}-\mathop{\sum }_{l=1}^{4}R_{l,i+1/2,j+1/2}^{n}.\end{eqnarray}$$
Multiplying (4.56) by
$[E_{c_{f}}+\unicode[STIX]{x1D6FF}_{c_{f}}]_{i+1/2,j+1/2}^{n}h_{i+1/2}k_{j+1/2}$
and making summation on
$i,~j$
for
$0\leqslant i\leqslant N_{x}-1,~0\leqslant j\leqslant N_{y}-1$
, we have that
$$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D713}^{n}d_{t}(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}+(Z^{n}d_{t}\unicode[STIX]{x1D713}^{n}-c_{f}^{n}d_{t}\unicode[STIX]{x1D719}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,(d_{x}E_{\boldsymbol{w}}^{x,n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}+(d_{y}E_{\boldsymbol{w}}^{y,n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle =-\frac{k_{c}k_{s}}{k_{c}+k_{s}}(a_{v}(\unicode[STIX]{x1D713}^{n})Z^{n}-a_{v}(\unicode[STIX]{x1D719}^{n})(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,f_{P}\Vert (E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})\Vert _{l^{2},M}^{2}-(E_{\unicode[STIX]{x1D719}}^{n}d_{t}(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \qquad -\mathop{\left(\mathop{\sum }_{l=1}^{4}R_{l}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\right)}\nolimits_{l^{2},M}\!.\end{eqnarray}$$
The first term on the left-hand side of (4.57) can be transformed into
$$\begin{eqnarray}\displaystyle & & \displaystyle (\unicode[STIX]{x1D719}^{n}[d_{t}(E_{c_{f}}+\unicode[STIX]{x1D6FF}_{c_{f}})]^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \quad =\frac{1}{2}d_{t}\Vert \sqrt{\unicode[STIX]{x1D719}^{n}}(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})\Vert _{l^{2},M}^{2}+\frac{\unicode[STIX]{x0394}t}{2}\Vert \sqrt{\unicode[STIX]{x1D719}^{n}}d_{t}(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})\Vert _{l^{2},M}^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad -\frac{1}{2}\Vert \sqrt{d_{t}\unicode[STIX]{x1D719}^{n}}(E_{c_{f}}^{n-1}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n-1})\Vert _{l^{2},M}^{2}.\end{eqnarray}$$
Recalling Lemma 8, the second term on the left-hand side of (4.57) can be estimated as
$$\begin{eqnarray}\displaystyle & & \displaystyle |(Z^{n}d_{t}\unicode[STIX]{x1D713}^{n}-c_{f}^{n}d_{t}\unicode[STIX]{x1D719}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant |(Z^{n}d_{t}\unicode[STIX]{x1D713}^{n}-c_{f}^{n}d_{t}\unicode[STIX]{x1D713}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}|+|(c_{f}^{n}d_{t}\unicode[STIX]{x1D713}^{n}-c_{f}^{n}d_{t}\unicode[STIX]{x1D719}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\Vert E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\Vert _{l^{2},M}^{2}+C\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C(h^{4}+k^{4}).\end{eqnarray}$$
Taking notice of Lemmas 1, 2 and (3.13), the third term on the left-hand side of (4.57) can be transformed into
$$\begin{eqnarray}\displaystyle (d_{x}E_{\boldsymbol{w}}^{x,n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M} & = & \displaystyle -(E_{\boldsymbol{w}}^{x,n},D_{x}(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}))_{l^{2},T,M}\nonumber\\ \displaystyle & = & \displaystyle \left(E_{\boldsymbol{w}}^{x,n},\frac{W^{x,n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{w^{x,n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}+(E_{\boldsymbol{w}}^{x,n},\unicode[STIX]{x1D716}^{x}(c_{f}))_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle +\left(E_{\boldsymbol{w}}^{x,n},\frac{\unicode[STIX]{x1D712}(U^{x,n})({\mathcal{P}}_{h}Z)^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{u^{x,n}c_{f}^{n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}.\end{eqnarray}$$
The first term on the right-hand side of (4.60) can be estimated as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(E_{\boldsymbol{w}}^{x,n},\frac{W^{x,n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{w^{x,n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \left(E_{\boldsymbol{w}}^{x,n},\frac{E_{\boldsymbol{w}}^{x,n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}\right)_{l^{2},T,M}+\left(E_{\boldsymbol{w}}^{x,n},\frac{w^{x,n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{w^{x,n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}\!.\end{eqnarray}$$
The last term on the right-hand side of (4.61) can be estimated as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left|\left(E_{\boldsymbol{w}}^{x,n},\frac{w^{x,n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{w^{x,n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \frac{1}{6D^{\ast }}\Vert E_{\boldsymbol{w}}^{x,n}\Vert _{l^{2},T,M}^{2}+C\Vert ({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}-({\mathcal{P}}_{h}\unicode[STIX]{x1D719})^{n}\Vert _{l^{2},T,M}^{2}+C\Vert ({\mathcal{P}}_{h}\unicode[STIX]{x1D719})^{n}-\unicode[STIX]{x1D719}^{n}\Vert _{l^{2},T,M}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \frac{1}{6D^{\ast }}\Vert E_{\boldsymbol{w}}^{x,n}\Vert _{l^{2},T,M}^{2}+C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C(h^{4}+k^{4}).\end{eqnarray}$$
The second term on the right-hand side of (4.60) can be estimated as
$$\begin{eqnarray}|(E_{\boldsymbol{w}}^{x,n},\unicode[STIX]{x1D716}^{x,n}(c_{f}))_{l^{2},T,M}|\leqslant \frac{1}{6D^{\ast }}\Vert E_{\boldsymbol{w}}^{x,n}\Vert _{l^{2},T,M}^{2}+O(h^{4}+k^{4}).\end{eqnarray}$$
The last term in the right-hand side of (4.60) can be estimated as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left|\left(E_{\boldsymbol{w}}^{x,n},\frac{\unicode[STIX]{x1D712}(U^{x,n})({\mathcal{P}}_{h}Z)^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{u^{x,n}c_{f}^{n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}\right|\leqslant \left|\left(E_{\boldsymbol{w}}^{x,n},\frac{U^{x,n}({\mathcal{P}}_{h}Z)^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{U^{x,n}c_{f}^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}\right)_{l^{2},T,M}\right|\nonumber\\ \displaystyle & & \displaystyle \qquad +\left|\left(E_{\boldsymbol{w}}^{x,n},\frac{(\unicode[STIX]{x1D712}(U^{x,n})-u^{x,n})c_{f}^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}\right)_{l^{2},T,M}\right|+\left|\left(E_{\boldsymbol{w}}^{x,n},\frac{u^{x,n}c_{f}^{n}}{({\mathcal{P}}_{h}\unicode[STIX]{x1D713})^{n}D_{11}^{n}}-\frac{u^{x,n}c_{f}^{n}}{\unicode[STIX]{x1D719}^{n}D_{11}^{n}}\right)_{l^{2},T,M}\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \frac{1}{6D^{\ast }}\Vert E_{\boldsymbol{w}}^{x,n}\Vert _{l^{2},T,M}^{2}+C\Vert E_{c_{f}}^{n}\Vert _{l^{2},M}^{2}+C\Vert E_{\tilde{\boldsymbol{u}}}^{x,n}\Vert _{l^{2},T,M}^{2}+C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+O(h^{4}+k^{4}),\end{eqnarray}$$
where the last inequality holds due to the fact that
$|\unicode[STIX]{x1D712}(U^{x,n})-u^{x,n}|\leqslant |U^{x,n}-u^{x,n}|$
pointwise. The last term on the left-hand side of (4.57) can be estimated similarly to (4.60)–(4.64).
Noticing that
$0\leqslant a_{v}(\unicode[STIX]{x1D713}^{n})\leqslant a_{0}$
and
$|a_{v}(\unicode[STIX]{x1D713}^{n})-a_{v}(\unicode[STIX]{x1D719}^{n})|\leqslant a_{0}/(1-\unicode[STIX]{x1D719}_{0})|E_{\unicode[STIX]{x1D719}}^{n}|$
, the first term on the right-hand side of (4.57) can be estimated as
$$\begin{eqnarray}\displaystyle & & \displaystyle \left|-\frac{k_{c}k_{s}}{k_{c}+k_{s}}(a_{v}(\unicode[STIX]{x1D713}^{n})Z^{n}-a_{v}(\unicode[STIX]{x1D719}^{n})(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \left|-\frac{k_{c}k_{s}}{k_{c}+k_{s}}(a_{v}(\unicode[STIX]{x1D713}^{n})(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\right|\nonumber\\ \displaystyle & & \displaystyle \qquad +\left|-\frac{k_{c}k_{s}}{k_{c}+k_{s}}((a_{v}(\unicode[STIX]{x1D713}^{n})-a_{v}(\unicode[STIX]{x1D719}^{n}))(c_{f}^{n}-\unicode[STIX]{x1D6FF}_{c_{f}}^{n}),E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})_{l^{2},M}\right|\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C\Vert E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\Vert _{l^{2},M}^{2}+C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}.\end{eqnarray}$$
The last term on the right-hand side of (4.57) can be estimated by
$$\begin{eqnarray}\left|-\mathop{\left(\mathop{\sum }_{l=1}^{4}R_{l}^{n},E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\right)}\nolimits_{l^{2},M}\right|\leqslant C\Vert E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\Vert _{l^{2},M}^{2}+O(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
Combining (4.57) with the above estimates, we have
$$\begin{eqnarray}\displaystyle & & \displaystyle \frac{1}{2}d_{t}\Vert \sqrt{\unicode[STIX]{x1D719}^{n}}(E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n})\Vert _{l^{2},M}^{2}+\frac{1}{2D^{\ast }}\Vert E_{\boldsymbol{w}}^{n}\Vert _{l^{2}}^{2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant \frac{1}{2}\Vert \sqrt{d_{t}\unicode[STIX]{x1D719}^{n}}(E_{c_{f}}^{n-1}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n-1})\Vert _{l^{2},M}^{2}+C\Vert E_{c_{f}}^{n}+\unicode[STIX]{x1D6FF}_{c_{f}}^{n}\Vert _{l^{2},M}^{2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,C\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{m}^{2}+C\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+C\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+C(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
Multiplying (4.67) by
$2\unicode[STIX]{x0394}t$
, summing for
$n$
from 1 to
$m$
,
$m\leqslant N$
and applying Gronwall’s lemma give that
$$\begin{eqnarray}\displaystyle \Vert E_{c_{f}}^{m}+\unicode[STIX]{x1D6FF}_{c_{f}}^{m}\Vert _{l^{2},M}^{2}+C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert E_{\boldsymbol{w}}^{n}\Vert _{l^{2}}^{2} & {\leqslant} & \displaystyle C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t(\Vert E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2}+\Vert d_{t}E_{\unicode[STIX]{x1D719}}^{n}\Vert _{l^{2},M}^{2})\nonumber\\ \displaystyle & & \displaystyle +\,C\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert E_{\tilde{\boldsymbol{u}}}^{n}\Vert _{l^{2}}^{2}+O(\unicode[STIX]{x0394}t^{2}+h^{4}+k^{4}).\end{eqnarray}$$
This completes the proof.
Based on the above lemmas, we can obtain
$$\begin{eqnarray}\displaystyle & & \displaystyle \Vert (P-p)^{m}\Vert _{l^{2},M}+\Vert (\boldsymbol{U}-\boldsymbol{u})^{m}\Vert _{l^{2}}+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert \unicode[STIX]{x1D63F}(\boldsymbol{U}-\tilde{\boldsymbol{u}})^{n}\Vert _{l^{2}}^{2}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\Vert (Z-c_{f})^{m}\Vert _{l^{2},M}+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert (\boldsymbol{W}-\boldsymbol{w})^{n}\Vert _{l^{2}}^{2}\right)^{1/2}+\Vert (\unicode[STIX]{x1D713}-\unicode[STIX]{x1D719})^{m}\Vert _{l^{2},M}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{2}),\quad m\leqslant N.\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert (\boldsymbol{U}-\tilde{\boldsymbol{u}})^{n}\Vert _{l^{2}}^{2}\right)^{1/2} & {\leqslant} & \displaystyle \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert \unicode[STIX]{x1D63F}(\boldsymbol{U}-\tilde{\boldsymbol{u}})^{n}\Vert _{l^{2}}^{2}\right)^{1/2}\nonumber\\ \displaystyle & {\leqslant} & \displaystyle C(\unicode[STIX]{x0394}t+h^{2}+k^{2}),\quad m\leqslant N.\end{eqnarray}$$
Recalling the definition of
$\tilde{\boldsymbol{u}}$
, we have the main error estimates below.
Theorem 12. Suppose (H1)–(H5) hold, then there exists a positive constant
$C$
independent of
$h,~k$
and
$\unicode[STIX]{x0394}t$
such that
$$\begin{eqnarray}\displaystyle & & \displaystyle \Vert (P-p)^{m}\Vert _{l^{2},M}+\Vert (\boldsymbol{U}-\boldsymbol{u})^{m}\Vert _{l^{2}}+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert (\boldsymbol{U}-\boldsymbol{u})^{n}\Vert _{l^{2}}^{2}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\Vert (Z-c_{f})^{m}\Vert _{l^{2},M}+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert (\boldsymbol{W}-\boldsymbol{w})^{n}\Vert _{l^{2}}^{2}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\Vert (\unicode[STIX]{x1D713}-\unicode[STIX]{x1D719})^{m}\Vert _{l^{2},M}\leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{2}),\quad m\leqslant N.\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert d_{x}(U^{x,n}-u^{x,n})\Vert _{l^{2},M}^{2}\right)^{1/2}\!+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert d_{y}(U^{y,n}-u^{y,n})\Vert _{l^{2},M}^{2}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{2}),\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & & \displaystyle \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{y}(U^{x,n}-u^{x,n})\Vert _{l^{2},T_{y}}^{2}\right)^{1/2}\!+\left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{x}(U^{y,n}-u^{y,n})\Vert _{l^{2},T_{x}}^{2}\right)^{1/2}\nonumber\\ \displaystyle & & \displaystyle \quad \leqslant C(\unicode[STIX]{x0394}t+h+k).\end{eqnarray}$$
If all
$k_{j+1/2}=k$
, i.e. uniform meshes are used in the
$y$
direction, we further have that
$$\begin{eqnarray}\displaystyle & \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{y}(U^{x,n}-u^{x,n})\Vert _{l^{2},T_{y}}^{2}\right)^{1/2}\leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{3/2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{y}(U^{x,n}-u^{x,n})\Vert _{l^{2},T}^{2}\right)^{1/2}\leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{2}). & \displaystyle\end{eqnarray}$$
If all
$h_{i+1/2}=h$
, i.e. uniform meshes are used in the
$x$
direction, we further have that
$$\begin{eqnarray}\displaystyle & \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{x}(U^{y,n}-u^{y,n})\Vert _{l^{2},T_{x}}^{2}\right)^{1/2}\leqslant C(\unicode[STIX]{x0394}t+h^{3/2}+k^{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \left(\mathop{\sum }_{n=1}^{m}\unicode[STIX]{x0394}t\Vert D_{x}(U^{y,n}-u^{y,n})\Vert _{l^{2},T}^{2}\right)^{1/2}\leqslant C(\unicode[STIX]{x0394}t+h^{2}+k^{2}). & \displaystyle\end{eqnarray}$$
Remark.
We can derive rigorous stability analysis for the finite difference scheme for wormhole propagation with the Darcy–Brinkman–Forchheimer framework by combining the procedure in the above proof of error estimates and that in Li & Rui (Reference Li and Rui2018a ). However, the process is very tedious, so we leave the details to the interested readers.
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
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D63F}=10^{-2}~\unicode[STIX]{x1D644},~K_{0}=1,~T=0.5,\\ \displaystyle \unicode[STIX]{x1D6FC}=k_{c}=k_{s}=\unicode[STIX]{x1D707}=\unicode[STIX]{x1D700}=\unicode[STIX]{x1D70C}=1,\\ \displaystyle a_{0}=0.5,~\unicode[STIX]{x1D70C}_{s}=10,~f_{P}=-1,\end{array}\right\} & \displaystyle\end{eqnarray}$$
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
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}l@{}}\Vert F-f\Vert _{l^{\infty }(i)}=\max _{1\leqslant n\leqslant N}\left\{\Vert (F-f)^{n}\Vert _{i}\right\},\quad i=l^{2},M~\text{or}~l^{2},\\ \Vert F-f\Vert _{l^{2}(i)}=\left(\unicode[STIX]{x0394}t\mathop{\sum }_{n=1}^{N}\left\Vert F^{n}-f^{n}\right\Vert _{i}^{2}\right)^{1/2},\quad i=l^{2}~\text{or}~l^{2},\unicode[STIX]{x1D709},\end{array}\right\} & \displaystyle\end{eqnarray}$$
where
$\unicode[STIX]{x1D709}=M,T,T_{x},T_{y}$
.
The non-uniform mesh generated from the
$20\times 20$
uniform mesh for Example 2.
Error and convergence rates of Example 1 on uniform grids.
Error and convergence rates of Example 1 on uniform grids.
Error and convergence rates of Example 1 on uniform grids.
Error and convergence rates of Example 1 on non-uniform grids.
Error and convergence rates of Example 1 on non-uniform grids.
Error and convergence rates of Example 1 on non-uniform grids.
Error and convergence rates of Example 1 on non-uniform grids.
Error and convergence rates of Example 1 on non-uniform grids.
Error and convergence rates of Example 2 on uniform grids.
Error and convergence rates of Example 2 on uniform grids.
Error and convergence rates of Example 2 on uniform grids.
Numerical porosity of Example 3.
Numerical porosity of Example 4.
Numerical concentration of Example 4.
Numerical pressure of Example 4.
Numerical velocity of Example 4.
Example 1. Here the initial condition and the right-hand side of the equation are computed according to the analytic solution given below.
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle p(\boldsymbol{x},t)=t(x^{2}-1/2)(y^{2}-1/2),\\ \displaystyle \boldsymbol{u}(\boldsymbol{x},t)=t(x^{2}(x-1)^{2}y(y-1)(2y-1),-y^{2}(y-1)^{2}x(x-1)(2x-1)),\\ \displaystyle c_{f}(\boldsymbol{x},t)=tx^{2}(1-x)^{2}y^{2}(1-y)^{2},\\ \displaystyle \unicode[STIX]{x1D719}(\boldsymbol{x},t)=1-\displaystyle \text{e}^{-(1/80)t^{2}x^{2}(1-x)^{2}y^{2}(1-y)^{2}\text{e}^{x+y+1}-(x+y+1)}.\end{array}\right\}\end{eqnarray}$$
Example 2. Here the initial condition and the right-hand side of the equation are computed according to the analytic solution given below.
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\displaystyle p(\boldsymbol{x},t)=t\sin (\unicode[STIX]{x03C0}y),\\ \displaystyle \boldsymbol{u}(\boldsymbol{x},t)=\unicode[STIX]{x03C0}t(\sin ^{2}(\unicode[STIX]{x03C0}x)\sin (2\unicode[STIX]{x03C0}y),-\sin (2\unicode[STIX]{x03C0}x)\sin ^{2}(\unicode[STIX]{x03C0}y)),\\ \displaystyle c_{f}(\boldsymbol{x},t)=t(x+1)^{2}(x-1)^{2}(y+1)^{2}(y-1)^{2},\\ \displaystyle \unicode[STIX]{x1D719}(\boldsymbol{x},t)=1-\displaystyle \text{e}^{-(1/80)t^{2}(x+1)^{2}(x-1)^{2}(y+1)^{2}(y-1)^{2}\text{e}^{x+y+3}-(x+y+3)}.\end{array}\right\}\end{eqnarray}$$
In this example
$\Vert D_{y}E_{\boldsymbol{u}}^{x}\Vert _{l^{\infty }(l^{2},T_{y})}$
also has second-order convergence on uniform grids in the spatial direction because
$\unicode[STIX]{x2202}^{2}u^{x}/\unicode[STIX]{x2202}y^{2}=0$
for
$y=-1$
and
$y=1$
.
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:
$$\begin{eqnarray}\displaystyle & \left.\begin{array}{@{}l@{}}\displaystyle \unicode[STIX]{x1D707}=10^{-2}~\text{Pa s},\quad \unicode[STIX]{x1D70C}_{s}=2500~\text{kg m}^{-2},\quad \unicode[STIX]{x1D70C}=10^{3}~\text{kg m}^{-2},\quad a_{0}=0.5~\text{m}^{-1},\\ \displaystyle \unicode[STIX]{x1D6FC}=0.1~\text{kg mol}^{-1},\quad k_{c}=0.1~\text{m s}^{-1},\quad c_{I}=10~\text{mol m}^{-2},\quad k_{s}=1~\text{m s}^{-1}.\end{array}\right\} & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{I}=\left\{\begin{array}{@{}l@{}}0.1~\text{m s}^{-1},\quad x=1.25\times 10^{-3},\\ 0,\quad \text{otherwise.}\end{array}\right. & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & f_{P}=\left\{\begin{array}{@{}l@{}}-0.1~\text{m s}^{-1},\quad x=1.9875\times 10^{-1},\\ 0,\quad \text{otherwise.}\end{array}\right. & \displaystyle\end{eqnarray}$$
Example 3. In this example,
$\unicode[STIX]{x1D63F}=10^{-7}~\unicode[STIX]{x1D644}$
. The distributions of initial porosity and permeability are listed as follows:
$$\begin{eqnarray}\left.\begin{array}{@{}l@{}}\unicode[STIX]{x1D719}_{0}=0.6,\quad K_{0}=10^{-6},\quad (x,y)=(1.25\times 10^{-3},5.125\times 10^{-2}),\\ \unicode[STIX]{x1D719}_{0}=0.5,\quad K_{0}=10^{-7},\quad (x,y)=(1.25\times 10^{-3},1.0125\times 10^{-1}),\\ \unicode[STIX]{x1D719}_{0}=0.2,\quad K_{0}=10^{-8},\quad \text{otherwise.}\end{array}\right\}\end{eqnarray}$$
Example 4. In this example,
$\unicode[STIX]{x1D63F}=10^{-7}~\unicode[STIX]{x1D644}$
. The initial porosity in the porous medium also obeys the uniform distribution and the range is from 0.05 to 0.35. The initial permeability of the porous media complies with the uniform distribution and the range is from
$10^{-8}$
to
$10^{-7}$
.
The distributions of porosity for examples 3 and 4 are calculated at the end of 25th time step, 50th time step, 75th time step and 100th time step, respectively. We show the distributions of porosity in figures 2 and 3. These results are computed on the grid of
$80\times 80$
cells. In addition, we demonstrate the configurations of concentration, pressure and flow streamlines for example 4 at the end of 25th time step, 50th time step, 75th time step and 100th time step, respectively, in figures 4–6.
From tables 1–16, we can see that the finite difference approximations on non-uniform staggered grids for the pressure, velocity, porosity, concentration and its flux have
$(\unicode[STIX]{x0394}t+h^{2}+k^{2})$
accuracy in discrete norms. These results are consistent with the error estimates in Theorem 12. From figures 2 and 3, it can be observed that the average porosities in all frameworks are increasing, which demonstrates that the matrix is dissolved by acid. Furthermore, we can see that at the beginning of the development, all frameworks display the formation of some small fingers. As time elapses, some major fingers appear. We can observe that the chemical concentration front advances in certain directions more than the other directions in figure 4. Also, as shown in figures 5 and 6, it is obvious that the heterogeneity of porosity and permeability in wormhole formation also has distinct effect on the velocity field, which in turn increases non-uniformity of the chemical reaction front. Those examples clearly show that the fully conservative finite difference scheme on non-uniform staggered grids can simulate wormhole propagation effectively.
Error and convergence rates of Example 2 on non-uniform grids.
Error and convergence rates of Example 2 on non-uniform grids.
Error and convergence rates of Example 2 on non-uniform grids.
Error and convergence rates of Example 2 on non-uniform grids.
Error and convergence rates of Example 2 on non-uniform grids.
6 Conclusion
In this paper, we have developed a finite difference scheme on non-uniform staggered grids for wormhole propagation with the Darcy–Brinkman–Forchheimer framework in porous media, which can be effectively applied to enhance oil and gas production rates. The DBF framework can simulate the convective flow comparatively well under conditions of variable porosity and deal with multiple complicating factors in wormhole propagation, including boundary-layer development, macroscopic and microscopic shear stress, and the inertial force. The coupled analysis approach to handle the fully coupling relation of multivariables is employed. To guarantee full mass conservation, we introduce an auxiliary flux variable. Based on the cutoff operator of velocity and the coupled analysis approach, also inspired by the analysis technique for Darcy–Forchheimer and Stokes equations, we obtain second-order superconvergence of the pressure, velocity, porosity, concentration and its flux in different discrete norms on non-uniform grids. We also obtain second-order superconvergence for some terms of the
$H^{1}$
norm of the velocity on non-uniform grids. Finally, some numerical experiments are presented to verify the theoretical analysis and effectiveness of the given scheme.
Acknowledgements
The authors would like to thank the editor and referees for their valuable comments and suggestions which helped us to improve the results of this paper. The first author is supported by the Postdoctoral Science Foundation of China grant no. BX20190187. The second author is supported by the National Natural Science Foundation of China grant no. 11671233.
