Multiblock Mortar Mixed Approach for Second Order Parabolic Problems

Abstract

In this paper, the multiblock mortar mixed approximation of second order parabolic partial differential equations is considered. In this method, the simulation domain is decomposed into the non-overlapping subdomains (blocks), and a physically-meaningful boundary condition is set on the mortar interface between the blocks. The governing equations hold locally on each subdomain region. The local problems on blocks are coupled by introducing a special approximation space on the interfaces of neighboring subdomains. Each block is locally covered by an independent grid and the standard mixed finite element method is applied to solve the local problem. The unique solvability of the discrete problem is shown, and optimal order convergence rates are established for the approximate velocity and pressure on the subdomain. Furthermore, an error estimate for the interface pressure in mortar space is presented. The numerical experiments are presented to validate the efficiency of the method.

1. Introduction

The numerical approximation of partial differential equations has received considerable attention due to its applications in diverse fields of science and technology. In modern scientific computation, efficiency and accuracy are the main objectives of numerical methods. Generally, numerical methods reduce the solution of partial differential equations into the solution of an algebraic system of equations. In practical applications, using direct methods, it is very complicated to attain both efficiency and accuracy at the same time. For example, fluid flow in an underground reservoir involves a highly heterogeneous medium. Furthermore, the physical properties of the media change rapidly on a large domain; for example, in fluid flow in porous media, the permeability of fluid fluctuates rapidly. In such cases, either accuracy is sacrificed or a very fine mesh is used to resolve heterogeneity in permeability. The latter gives rise to a huge system of coupled equations, which in many cases becomes very challenging to solve computationally. These difficulties motivated discovering alternate techniques that keep the computational burden manageable along with the maximum possible accuracy.

Several domain decomposition methods have been developed to overcome these difficulties. The purpose of all domain decomposition methods is to provide efficient approximation by alleviating the computation burden. This is achieved by splitting the physical problem into smaller subproblems by dividing the domain into a series of subdomains. This technique allows considering different physical models in different regions. It is also capable of employing the most suitable approximation method in different blocks of the computational domain. We consider the second order linear parabolic partial differential equation modeling single-phase Darcy flow in a porous medium written as a system of first order equations:

$$\begin{array}{}\mathrm{(1)}& \hfill \mathit{u}& =-a\left(x\right)\nabla p,\hfill & \mathrm{in}\Omega \times J,\hfill \mathrm{(2)}& \hfill \frac{\partial p}{\partial t}+\nabla ·\mathit{u}& =f\left(x\right)\hfill & \mathrm{in}\Omega \times J,\hfill \mathrm{(3)}& \hfill p& =g\hfill & \mathrm{on}\partial \Omega \times J,\hfill \mathrm{(4)}& \hfill p& =p_0,\hfill & \mathrm{in}\Omega \times \left\{0\right\},\hfill \end{array}$$

where $\Omega \subset {\mathit{R}}^d$, $d=2$ or 3, and $J=[0,T]$. Equation (1) is the Darcy law, and (2) is the mass conservation. The unknown $\mathit{u}$ represents Darcy velocity, and p is the pressure. Assume that a has the continuous bounded derivative up to order two and there exist constants $\alpha$ and $\beta$ such that:

$$0<\alpha \le a\le \beta <\infty .$$

(5)

Moreover, assume that there exist constants $B1$ and $B2$ such that:

$$\begin{array}{c}\left|f\right|,\frac{\partial f}{\partial x},\frac{\partial f}{\partial y}\le B1,\hfill \\ \left|g\right|,\frac{\partial g}{\partial x},\frac{\partial g}{\partial y}\le B2.\hfill \end{array}$$

Mixed finite element methods have emerged with considerable popularity in several areas of science and engineering. The mixed methods enjoy local mass conservation and provide accurate approximation of two variables of physical interest simultaneously, for example velocity and pressure in fluid flow. They are also capable of approximating both variables with the same accuracy. The mixed finite element methods have been extensively used to solve the elliptic [1,2,3,4,5,6,7,8,9,10,11,12,13,14] and parabolic [15,16,17,18,19] partial differential equations. The mixed finite element methods are also used to approximate the partial differential equations governing fluid flow [20,21,22,23,24,25,26,27,28].

The mortar finite element method was first introduced in [29] for the elliptic partial differential equation. It is a nonconforming domain decomposition method in which the weak continuity condition is imposed across the subdomain interfaces. The mixed finite element method on locally-nested refined grids was developed in [30,31]. In this technique, the continuity of flux across the subdomain block interface is enforced by using the concept of slave or worker nodes. Since in this approach, the grids must be nested, it cannot be extended to the non-matching ones.

The mixed version of the mortar method is proposed and analyzed in [20]. The idea was to divide the original domain into smaller blocks and impose the governing partial differential equation on each block. An unknown pressure is introduced on the interfaces between the subdomains known as the Lagrange multiplier. This new variable provides the Dirichlet boundary condition for the local problem on each subdomain block. This procedure splits the large problem into a series of small subproblems. As a result, the local small-sized problems can be solved more comfortably as compared to a single large problem on the entire domain. This method gives better approximation on non-matching grids. Since different grids are used on the adjacent subdomains, the normal trace of the velocity space cannot be used as a Lagrange multiplier space; therefore, a different space, namely the mortar finite element space, is defined on the inter-block boundaries. The method gives the optimal order convergence if the mortar space approximation is one order higher than the normal trace of the subdomain velocity space. The multiblock approach provides the independent partition of each subdomain block so that an appropriate approximation method can be applied to solve the local problem. This approach offers great flexibility to formulate the different physical and mathematical models on different subdomain blocks. Moreover, the non-overlapping domain decomposition algorithms can be applied to implement the scheme in parallel. The multiblock technique has been further explored to solve several problems, for example [32,33,34].

To the best of the authors’ knowledge, there is no article in the literature that presents a complete analysis of the multiblock mortar mixed approach for the time-dependent problem under consideration. The present study is an extension of the results established in [20] to the parabolic equations. The continuous time discretization is considered, and the difficulties arising in the multiblock mortar mixed formulation of time-dependent problems are addressed. First, the existence and uniqueness of the discrete solution is accomplished by using the established theory for ordinary differential equations. Second, the convergence analysis of parabolic problems requires the construction of elliptic projection [35], which maps the continuous solution into the discrete space. The a priori error estimate is derived via mixed elliptic projection for the approximate velocity and pressure on the subdomain. The error estimate for the unknown mortar pressure through the interface bilinear form as described in [20] is not applicable in our case. Therefore, an alternate way is used, and an error bound for the Lagrange multiplier serving as the Dirichlet boundary condition on mortar the interface is proven. The error estimate for the mortar interface variable predicts the accuracy of data flow across the subdomain interface.

Consider the domain $\mathcal{D}\subset {\mathit{R}}^2$, and let us denote by $\partial \mathcal{D}$ the boundary of $\mathcal{D}$. Let $\mathit{n}$ denote the unit outward normal to $\partial \mathcal{D}$. Throughout the article, we shall use ${(·,·)}_{\mathcal{D}}$ to denote the $L^2\left(\mathcal{D}\right)$ or $L^2{\left(\mathcal{D}\right)}^2$ the inner product on $\mathcal{D}$. Furthermore the inner product on $\partial \mathcal{D}$ is denoted by ${\langle ·,·\rangle }_{\partial \mathcal{D}}$. The $L^2\left(\mathcal{D}\right)$ or $L^2{\left(\mathcal{D}\right)}^2$ norm will be denoted by $∥·∥$. For any non-negative integer ℓ and $1\le q\le \infty$, let $W^{\ell ,q}\left(\mathcal{D}\right)=\{\psi \in L^q\left(\mathcal{D}\right)\mid D^{\alpha }\left(\psi \right)\in L^q\left(\mathcal{D}\right)\}$ denote the standard Sobolev space. The norm and semi-norm on $W^{\ell ,q}\left(\mathcal{D}\right)$ are denoted by ${∥·∥}_{\ell ,q,\mathcal{D}}$ and ${|·|}_{\ell ,q,\mathcal{D}}$, respectively. In particular, for $q=2,W^{\ell ,2}\left(\mathcal{D}\right)=H^{\ell }\left(\mathcal{D}\right)$ denotes the space of ℓ times differentiable $L^2\left(\mathcal{D}\right)$ functions. The norm on $H^{\ell }\left(\mathcal{D}\right)$ is denoted by ${∥·∥}_{\ell ,\mathcal{D}}.$ In the notation below, we shall omit the subscript $\mathcal{D}$ when $\mathcal{D}=\Omega$. Throughout the article, the boldface letters are used to denote the vectors. Furthermore, for the weak formulation of mixed finite element methods, the usual function space is defined as:

$$\mathit{H}(div;\mathcal{D})=\{\psi \in L^2{\left(\mathcal{D}\right)}^2\mid \nabla ·\psi \in L^2\left(\mathcal{D}\right)\}.$$

Finally, for a given norm space $\chi =\chi \left(\mathcal{D}\right)$ equipped with the norm ${∥·∥}_{\chi }$, the function space $W^{\ell ,q}(J;\chi \left(\mathcal{D}\right))$ denotes the space of functions $\psi :[0,T]⟶\chi \left(\mathcal{D}\right)$ equipped with the norm, for $1\le q<\infty$:

$${∥\psi ∥}_{W^{\ell ,q}(J;\chi \left(\mathcal{D}\right))}={\left\{{\displaystyle \sum _{i=0}^{\ell }}{\int }_J{∥\frac{{\partial }^i}{{\partial t}^i}\psi (·,t)∥}_{\chi (\mathcal{D})}^{\ell }dt\right\}}^{1/\ell },$$

and if $q=\infty$, the integral is replaced by an essential supremum.

The remainder of the article is outlined as follows. In the next section, we discuss the domain decomposition and the method formulation and present the unique solvability of the discrete problem. The mixed method elliptic projection and related error estimates are given in Section 3. In Section 4, we provide optimal order convergence results. Section 5 deals with the error estimates for unknown mortar. The numerical experiments are provided in Section 6. Finally, the paper is concluded in Section 7.

2. Domain Decomposition and Method Formulation

Let us consider the non-overlapping decomposition of the domain $\Omega$ into the series of subdomains ${\left\{{\Omega }_i\right\}}_{i=1}^{nb}$ such that:

$$\overline{\Omega }={\cup }_{i=1}^{nb}{\Omega }_i\mathrm{and}{\Omega }_i\cap {\Omega }_j\ne \varnothing .$$

Let us define the interfaces:

$$\begin{array}{c}\hfill {\Gamma }_{ij}:=\partial {\Omega }_i\cap \partial {\Omega }_j,\Gamma :={\cup }_{1\le i<j\le nb}{\Gamma }_{ij},{\Gamma }_i:=\partial {\Omega }_i\cap \Gamma =\partial {\Omega }_i\backslash \partial \Omega .\end{array}$$

Here, ${\Gamma }_{ij}$ denotes the interface between two subdomains, $\Gamma$ represents the interior skeleton of decomposition, and ${\Gamma }_i$ is the interface related to a single subdomain in the decomposition.

We now define the function spaces to formulate the weak and discrete problems. Let us define the local weak spaces as ${\mathit{S}}_i=\mathit{H}(div;{\Omega }_i),U_i=L^2\left({\Omega }_i\right),{\Sigma }_{ij}=H^{1/2}\left({\Gamma }_{ij}\right).$ Then, their global counterparts are given by:

$$\begin{array}{c}\hfill \mathit{S}={\oplus }_{i=1}^{nb}{\mathit{S}}_i,U={\oplus }_{i=1}^{nb}{U}_i=L^2(\Omega ),\Sigma ={\oplus }_{1\le i<j\le nb}{\Sigma }_{ij}.\end{array}$$

In order to define the discrete finite element spaces, we first partition each subdomain into the local conforming and quasi-uniform mesh denoted by ${\mathcal{T}}_{h,i}$, that is $\overline{\Omega }={\cup }_{E\in {\mathcal{T}}_{h,i}}E$, where E denotes an element in the mesh with maximum diameter $h_i$. For the purpose of analysis, we define $h={max}_{1\le i\le nb}{h}_i$. We allow the possibility that the local grids ${\mathcal{T}}_i$ and ${\mathcal{T}}_j$ on adjacent subdomains ${\Omega }_i$ and ${\Omega }_j$ do not necessarily align on the interface. Let us denote by ${\mathcal{T}}_h$ the union of partitions ${\mathcal{T}}_{h,i}$ on subdomains ${\Omega }_i$, $1\le i\le nb$, that is ${\mathcal{T}}_h:={\cup }_{i=1}^{nb}{\mathcal{T}}_{h,i}.$ Let us consider the mixed finite element spaces ${\mathit{S}}_{h,i}\subset {\mathit{S}}_i,U_{h,i}\subset U_i$ of order $k+1$. These mixed finite element spaces can be either Raviart–Thomas ($R{T}_k$) [12], or Brezzi–Douglas–Fortin–Marini ($BDF{M}_k$) [2], or Brezzi–Douglas–Fortin ($BD{F}_k$) [1]. Moreover, we assume that the order of approximation is the same on each subdomain. It is well known that the property $\nabla ·{\mathit{S}}_{h,i}=U_{h,i}$ holds in standard mixed finite element spaces. Then, the corresponding global mixed spaces on ${\mathcal{T}}_h$ are defined by:

$${\mathit{S}}_h={\oplus }_{i=1}^{nb}{\mathit{S}}_{h,i},U_h={\oplus }_{i=1}^{nb}{U}_{h,i}.$$

Note that the normal components of vectors in the space ${\mathit{S}}_h$ are continuous between the elements of local grid on each subdomain block ${\Omega }_i$, but not across the interface $\Gamma$.

We consider the quasi-uniform mesh on mortar interface ${\Gamma }_{ij}$, denoted by ${\mathcal{T}}_{h,ij}$. This mortar interface mesh will be used to define a finite element space that provides the approximation for the Dirichlet boundary condition assigned to the mortar interface. Let us denote by ${\Sigma }_{h,ij}\subset L^2\left({\Gamma }_{ij}\right)$, the mortar space consisting either of continuous or discontinuous piecewise polynomials of degree $k+1$ on the interface mesh, where k is associated with the degree of polynomial in ${\mathit{S}}_h.\mathit{n}$. Note that for $d=2$ and the lowest order Raviart–Thomas space $\left({RT}_0\right)$, the mortar space ${\Sigma }_{h,ij}$ contains piecewise linear polynomials. In the case of $d=3$ and $R{T}_0$, the mortar space consists of a bilinear polynomial on the interface. Then, the global mortar space on $\Gamma$ is given by:

$${\Sigma }_h={\oplus }_{1\le i<j\le nb}{\Sigma }_{h,ij}.$$

In the following, we treat any function $\mu \in \Sigma$ as extended by zero on $\partial \Omega$. Also in the case of discontinuous space, ${\mathcal{T}}_{h,ij}$ need not to be conforming.

A unique solution for the discrete problem exists, if the following technical condition holds.

Assumption 1.

Let ${\mathcal{Q}}_h^i:L^2\left({\Gamma }_i\right)⟶{\mathit{S}}_{h,i}·{\mathit{n}}_i{|}_{{\Gamma }_i}$ be the $L^2$- orthogonal projection on each subdomain ${\Omega }_i$, $1\le i\le nb$, satisfying:

$$\begin{array}{c}\hfill {\langle \psi -{\mathcal{Q}}_h^i\psi ,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}=0,\mathit{v}\in {\mathit{S}}_{h,i},\psi \in L^2\left({\Gamma }_i\right),\end{array}$$

(6)

then assume that there exists a constant C, independent of h, such that for all $\mu \in {\Sigma }_h$:

$$\begin{array}{c}\hfill {\parallel \mu \parallel }_{{\Gamma }_{ij}}\le C(\parallel {\mathcal{Q}}_h^i{\mu \parallel }_{{\Gamma }_{ij}}+\parallel {\mathcal{Q}}_h^j\mu {\parallel }_{{\Gamma }_{ij}}),1\le i<j\le n.\end{array}$$

(7)

The condition (7) enforces a limit on the mortar degrees of freedom, and it can easily be satisfied in practice [36,37].

In the analysis below, several projection operators will be needed. We now introduce these projection operators and state their approximation properties. Recall that in the mixed finite element spaces, there exists the projection operator ${\Pi }_i:{\left(H^{1/2+ϵ}\left({\Omega }_i\right)\right)}^d\cap {\mathit{S}}_i⟶{\mathit{S}}_{h,i}$ satisfying for any $\mathit{\psi }\in {\left(H^{1/2+ϵ}\left({\Omega }_i\right)\right)}^d\cap {\mathit{S}}_i$ and $ϵ>0$:

$$\begin{array}{}\mathrm{(8)}& {(\nabla ·({\Pi }_i\mathit{\psi }-\mathit{\psi }),w)}_{{\Omega }_i}=0,\hfill & w\in U_h,\hfill \mathrm{(9)}& {\langle (\mathit{\psi }-{\Pi }_i\mathit{\psi })·{\mathit{n}}_i,\mathit{v}·{\mathit{n}}_i\rangle }_{\partial {\Omega }_i}=0,\hfill & \mathit{v}\in {\mathit{S}}_{h,i}.\hfill \end{array}$$

The projection operators ${\Pi }_i$ and ${\mathcal{Q}}_h^i$ defined above have the following approximation properties:

$$\begin{array}{}\mathrm{(10)}& {∥\mathit{\psi }-{\Pi }_i\mathit{\psi }∥}_{{\Omega }_i}\le C{∥\mathit{\psi }∥}_{r,{\Omega }_i}{h}^r,\hfill & 1\le r\le k+1,\hfill \mathrm{(11)}& {∥\nabla ·(\mathit{\psi }-{\Pi }_i\mathit{\psi }∥}_{{\Omega }_i}\le C{∥\nabla ·\mathit{\psi }∥}_{r,{\Omega }_i}{h}^r,\hfill & 0\le r\le l+1,\hfill \mathrm{(12)}& {∥(\mathit{\psi }-{\Pi }_i\mathit{\psi })·{\mathit{n}}_i∥}_{-s,{\Gamma }_{ij}}\le C{∥\mathit{\psi }∥}_{r,{\Gamma }_{ij}}{h}^{r+s},\hfill & 0\le r\le k+1,0\le s\le k+1,\hfill \mathrm{(13)}& {∥\psi -{\mathcal{Q}}_h^i\psi ∥}_{-s,{\Gamma }_{ij}}\le C{∥\psi ∥}_{r,{\Gamma }_{ij}}{h}^{r+s},\hfill & 0\le r\le k+1,0\le s\le k+1.\hfill \end{array}$$

Furthermore, if $\mathit{\psi }\in {\left(H^{ϵ}\left({\Omega }_i\right)\right)}^d,0<ϵ<1$, [38,39]:

$$\begin{array}{}\mathrm{(14)}& {∥{\Pi }_i\mathit{\psi }∥}_{{\Omega }_i}\le C{∥\mathit{\psi }∥}_{ϵ,{\Omega }_i},\hfill \mathrm{(15)}& {∥{\Pi }_i\mathit{\psi }-\mathit{\psi }∥}_{{\Omega }_i}\le C{h}^{ϵ}{∥\mathit{\psi }∥}_{ϵ,{\Omega }_i},\hfill \mathrm{(16)}& {∥{\Pi }_i\mathit{\psi }∥}_{{\Omega }_i}\le C({∥\mathit{\psi }∥}_{ϵ,{\Omega }_i}+{∥\nabla ·\mathit{\psi }∥}_{{\Omega }_i}).\hfill \end{array}$$

Here, ${\parallel ·\parallel }_r$ is the $H^r$-norm and ${\parallel ·\parallel }_{-s}$ is the norm of $H^{-s}$, the dual of $H^s$.

In addition, there exists the $L^2$-orthogonal projection operator $P_h:U⟶U_h$ such that for any $w\in U$, $P_{h}w\in U_h$ is defined by:

$$\begin{array}{cc}\hfill & (w-P_{h}w,\xi )=0,\xi \in U_h,\hfill \end{array}$$

(17)

and satisfies:

$$\begin{array}{}\mathrm{(18)}& (\nabla ·\mathit{v},w-P_{h}w)=0,\hfill & \mathit{v}\in {\mathit{S}}_h,w\in U,\hfill \mathrm{(19)}& ∥w-P_{h}w∥\le C{∥w∥}_r{h}^r,\hfill & 0\le r\le l+1,\hfill \end{array}$$

where l is the degree of the polynomial in $U_h$. Moreover, the above-defined projection operator ${\Pi }_i$ and $P_h$ satisfies:

$$\begin{array}{}\mathrm{(20)}& \nabla ·{\Pi }_i\mathit{\psi }=P_h\nabla ·\mathit{\psi },\hfill \mathrm{(21)}& \left({\Pi }_i\mathit{\psi }\right)·{\mathit{n}}_i={\mathcal{Q}}_h^i(\mathit{\psi }·{\mathit{n}}_i),\hfill \end{array}$$

where ${\mathcal{Q}}_h^i$ is defined in (6).

Finally, let ${\mathcal{P}}_h$ be the $L^2(\Gamma )$ projection onto ${\Sigma }_h$ satisfying:

$$\begin{array}{}\mathrm{(22)}& {\langle \phi -{\mathcal{P}}_h\phi ,\mu \rangle }_{\Gamma }=0,\hfill & \phi \in L^2(\Gamma ),\mu \in {\Sigma }_h,\hfill \mathrm{(23)}& {∥\phi -{\mathcal{P}}_h\phi ∥}_{-s,{\Gamma }_{ij}}\le C{∥\phi ∥}_{r,{\Gamma }_{ij}}{h}^{r+s},\hfill & 0\le r\le k+2,0\le s\le k+2.\hfill \end{array}$$

Note that the bounds (11)–(13), (19), and (23), are $L^2$-projection approximations [40]. The bound (10) can be found in [3,13].

In the subsequent analysis, we shall make use of the following inequalities:

$$\begin{array}{}\mathrm{(24)}& {∥\psi ∥}_{r,{\Gamma }_{ij}}\le C{∥\psi ∥}_{r+1/2,{\Omega }_i},\hfill \mathrm{(25)}& {∥\mathit{\psi }·\mathit{n}∥}_{\partial {\Omega }_i}\le C{h}^{-1/2}{∥\mathit{\psi }∥}_{{\Omega }_i}.\hfill \end{array}$$

Equation (24) is the nonstandard trace theorem (see [38], Theorem 1.5.2.1), and (25) is the local inverse inequality (see [20], Lemma 4.1).

2.1. Multiblock Variational Formulation

In the above-defined framework, the weak formulation of (1)–(4) is stated as:

Seeking the map $\{\mathit{u},p,\lambda \}:J⟶\mathit{S}\times U\times \Sigma$, such that:

$$\begin{array}{}\mathrm{(26)}& {(a^{-1}\mathit{u},v)}_{{\Omega }_i}={(p,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle \lambda ,\mathit{v}·n_i\rangle }_{{\Gamma }_i}-{\langle g,\mathit{v}·n_i\rangle }_{\partial {\Omega }_i\backslash \Gamma },\hfill & v\in {\mathit{S}}_i,\hfill \mathrm{(27)}& {(\frac{\partial p}{\partial t},w)}_{{\Omega }_i}+{(\nabla ·\mathit{u},w)}_{{\Omega }_i}={(f,w)}_{{\Omega }_i},\hfill & w\in U_i,\hfill \mathrm{(28)}& {\displaystyle \sum _{i=1}^{nb}}{\langle \mathit{u}·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0,\hfill & \mu \in \Sigma ,\hfill \mathrm{(29)}& p=p_0,\hfill & t=0,\hfill \end{array}$$

where ${\mathit{n}}_i$ denotes the outer unit normal to $\partial {\Omega }_i$ and $\lambda$ is the pressure on interfaces ${\Gamma }_i$. Equation (26) is obtained using integration by parts, and Equation (28) imposes the weak continuity of normal flux across the interface. Furthermore, we assume that the local problem over each subdomain ${\Omega }_i$ is at least $H^{3/2+ϵ}$ regular.

2.2. Multiblock Mixed Approximation

In the mixed finite element approximation of (26)–(29), seeking for $\{{\mathit{u}}_h,p_h,{\lambda }_h\}:J⟶{\mathit{S}}_h\times U_h\times {\Sigma }_h$ such that for each ${\Omega }_i$, $1\le i\le nb$,

$$\begin{array}{}\mathrm{(30)}& {(a^{-1}{\mathit{u}}_h,\mathit{v})}_{{\Omega }_i}={(p_h,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle {\lambda }_h,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}-{\langle g,\mathit{v}·{\mathit{n}}_i\rangle }_{\partial {\Omega }_i\backslash \Gamma },\hfill & \mathit{v}\in {\mathit{S}}_{h,i},\hfill \mathrm{(31)}& {(\frac{\partial p_h}{\partial t},w)}_{{\Omega }_i}+{(\nabla ·{\mathit{u}}_h,w)}_{{\Omega }_i}={(f,w)}_{{\Omega }_i},\hfill & w\in U_{h,i},\hfill \mathrm{(32)}& {\displaystyle \sum _{i=1}^{nb}}{\langle {\mathit{u}}_h·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0,\hfill & \mu \in {\Sigma }_h,\hfill \end{array}$$

with the initial condition $p_h\left(0\right)=P_h{p}_0$, the $L^2(\Omega )$-projection of initial data function $p_0$ onto $U_h$. We applied the standard mixed finite element over each subdomain ${\Omega }_i$. The local mass conservation over each subdomain grid cell is enforced by Equation (31). Note that the variable ${\lambda }_h$ represents the approximation to pressure on the mortar interface, and the continuity of normal flux across the interface is imposed by Equation (32). For the sake of analysis, let us define the weak continuous velocities [20]:

$${\mathit{S}}_h^0=\left\{\mathit{v}\in {\mathit{S}}_h:{\displaystyle \sum _{i=1}^{nb}}{\langle \mathit{v}{\mid }_{{\Omega }_i}·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0\forall \mu \in {\Sigma }_h\right\}.$$

(33)

The following lemma can be found in [20].

Lemma 1.

Under the assumption (7), there exists a projection operator ${\Pi }_0:\left(H^{1/2+ϵ}(\Omega )\right)\cap \mathit{S}⟶{\mathit{S}}_h^0$ such that:

$${(\nabla ·({\Pi }_0\mathit{\psi }-\mathit{\psi }),\tau )}_{\Omega }=0,\tau \in U_h,$$

(34)

satisfying:

$$\begin{array}{}\mathrm{(35)}& ∥{\Pi }_0\mathit{\psi }-\Pi \mathit{\psi }∥\le C{\displaystyle \sum _{i=1}^{nb}}{∥\mathit{\psi }∥}_{r+1/2,{\Omega }_i}{h}^{r+1/2},\hfill & 0\le r\le k+1,\hfill \mathrm{(36)}& ∥{\Pi }_0\mathit{\psi }-\mathit{\psi }∥\le C{\displaystyle \sum _{i=1}^{nb}}{∥\mathit{\psi }∥}_{r,{\Omega }_i}{h}^r,\hfill & 1\le r\le k+1.\hfill \end{array}$$

By using the weak velocities (33), unknown mortar pressure ${\lambda }_h$ can be eliminated, then our method (30)–(32) takes the following equivalent form:

$$\begin{array}{}\mathrm{(37)}& (a^{-1}{\mathit{u}}_h,\mathit{v})={\displaystyle \sum _{i=1}^{nb}}\left\{{(p_h,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle g,\mathit{v}·{\mathit{n}}_i\rangle }_{\partial {\Omega }_i\backslash \Gamma }\right\},\hfill & \mathit{v}\in {\mathit{S}}_h^0,\hfill \mathrm{(38)}& {\displaystyle \sum _{i=1}^{nb}}\left\{{(\frac{\partial p_h}{\partial t},w)}_{{\Omega }_i}+{(\nabla ·{\mathit{u}}_h,w)}_{{\Omega }_i}\right\}=(f,w),\hfill & w\in U_h.\hfill \end{array}$$

Next, we show the existence and uniqueness of the solution for the semidiscrete scheme (37) and (38). By using the basis function in the finite element spaces, the system (37) and (38) can be written as follows. Let N denote the number of edges and M be the number of elements in finite element discretization. Let $\phi$ and $\psi$ denote the basis functions in ${\mathit{S}}_h^0$ and $U_h$, respectively. Then, the finite element solutions ${\mathit{u}}_h$ and $p_h$ are written as:

$$\begin{array}{}\mathrm{(39)}& {\mathit{u}}_h={\displaystyle \sum _{i=1}^N}{\tilde{\mathit{u}}}_{h,i}\left(t\right){\phi }_i\left(x\right),\hfill \mathrm{(40)}& p_h={\displaystyle \sum _{i=1}^M}{\tilde{p}}_{h,i}{\psi }_i\left(x\right).\hfill \end{array}$$

Using (39) and (40), our method (37) and (38) can be written in the matrix form:

$$\begin{array}{}\mathrm{(41)}& A{\tilde{\mathit{u}}}_h-B{\tilde{p}}_h=-G,\hfill \mathrm{(42)}& C{\tilde{p}}_{ht}+B^T{\tilde{\mathit{u}}}_h=F,\hfill \end{array}$$

with ${\tilde{p}}_h\left(0\right)=0$ and ${\tilde{p}}_{ht}=\frac{\partial \widetilde{p_h}}{\partial t}.$ Observe that the matrices A and C are invertible, then (41) gives the relation:

$${\tilde{\mathit{u}}}_h=A^{-1}B{\tilde{p}}_h-A^{-1}G.$$

(43)

Combining (43) and (42), we arrive at:

$$C{\tilde{p}}_{ht}+B^T{A}^{-1}B{\tilde{p}}_h=F+B^T{A}^{-1}G.$$

(44)

Note that (44) is an ordinary partial differential equation and has a unique solution with smoothness conditions. Therefore, the discrete system (37) and (38) has the unique solution.

3. Elliptic Projection

For the purpose of analysis, we introduce the elliptic projection $\{{\mathcal{R}}_h,{\mathcal{S}}_h,{\mathcal{L}}_h\}$ of solution $\{\mathit{u},p,\lambda \}$ into the space ${\mathit{S}}_h\times U_h\times {\Sigma }_h$ by mapping $\{{\mathcal{R}}_h,{\mathcal{S}}_h,{\mathcal{L}}_h\}:J⟶{\mathit{S}}_h\times U_h\times {\Sigma }_h$ given by:

$$\begin{array}{}\mathrm{(45)}& \begin{array}{c}{(a^{-1}{\mathcal{R}}_h\mathit{u},\mathit{v})}_{{\Omega }_i}={({\mathcal{S}}_{h}p,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle {\mathcal{L}}_h\lambda ,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\hfill \\ -{\langle g,\mathit{v}·{\mathit{n}}_i\rangle }_{\partial {\Omega }_i\backslash \Gamma },\hfill \end{array}\hfill & \mathit{v}\in {\mathit{S}}_{h,i},\hfill \mathrm{(46)}& {(\nabla ·{\mathcal{R}}_h\mathit{u},w)}_{{\Omega }_i}={(f-\frac{\partial p}{\partial t},w)}_{{\Omega }_i},\hfill & w\in U_{h,i},\hfill \mathrm{(47)}& {\displaystyle \sum _{i=1}^{nb}}{\langle {\mathcal{R}}_h\mathit{u}·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0,\hfill & \mu \in {\Sigma }_h.\hfill \end{array}$$

Subtracting (45)–(47) from (26)–(28), the elliptic projection satisfies the system:

$$\begin{array}{}\mathrm{(48)}& {(a^{-1}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),\mathit{v})}_{{\Omega }_i}={(p-{\mathcal{S}}_{h}p,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle \lambda -{\mathcal{L}}_h\lambda ,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i},\hfill & \mathit{v}\in {\mathit{S}}_{h,i},\hfill \mathrm{(49)}& {(\nabla ·(\mathit{u}-{\mathcal{R}}_h\mathit{u}),w)}_{{\Omega }_i}=0,\hfill & w\in U_{h,i},\hfill \mathrm{(50)}& {\displaystyle \sum _{i=1}^{nb}}{\langle (\mathit{u}-{\mathcal{R}}_h\mathit{u})·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0,\hfill & \mu \in {\Sigma }_h.\hfill \end{array}$$

By using the weak velocities defined in (33), the system (48)–(50) has the following equivalent form:

$$\begin{array}{}\mathrm{(51)}& (a^{-1}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),\mathit{v})={\displaystyle \sum _{i=1}^{nb}}\left\{{(p-{\mathcal{S}}_{h}p,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle p,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\right\},\hfill & \mathit{v}\in {\mathit{S}}_h^0,\hfill \mathrm{(52)}& {\displaystyle \sum _{i=1}^{nb}}{(\nabla ·(\mathit{u}-{\mathcal{R}}_h\mathit{u}),w)}_{{\Omega }_i}=0,\hfill & w\in U_h.\hfill \end{array}$$

Then, the error estimates stated in Theorem 1 are derived in [20].

Theorem 1.

For the velocity ${\mathcal{R}}_h\mathit{u}$ and pressure ${\mathcal{S}}_{h}p$ of mixed elliptic projection (48)–(50), if Condition (7) holds, then there exists a positive constant independent of h such that:

$$\begin{array}{}\mathrm{(53)}& ∥\nabla ·(\mathit{u}-{\mathcal{R}}_h\mathit{u})∥\le C{\displaystyle \sum _{i=1}^{nb}}{∥\nabla ·\mathit{u}∥}_{r,{\Omega }_i}{h}^r,\hfill & 1\le r\le l+1,\hfill \mathrm{(54)}& ∥\mathit{u}-{\mathcal{R}}_h\mathit{u}∥\le C{\displaystyle \sum _{i=1}^{nb}}[{∥p∥}_{r+1,{\Omega }_i}+{∥\mathit{u}∥}_{r,{\Omega }_i}]h^r,\hfill & 1\le r\le k+1,\hfill \mathrm{(55)}& \parallel p-{\mathcal{S}}_{h}p\parallel \le C{\displaystyle \sum _{i=1}^{nb}}{[\parallel p\parallel }_{r+1,{\Omega }_i}+{\parallel \mathit{u}\parallel }_{r,{\Omega }_i}{+\parallel \nabla ·\mathit{u}\parallel }_{r,{\Omega }_i}]h^r,\hfill & 1\le r\le min(k+1,l+1).\hfill \end{array}$$

Sequentially, we shall need the error estimates for $\frac{\partial }{\partial t}(\mathit{u}-{\mathcal{R}}_h\mathit{u})$ and $\frac{\partial }{\partial t}(p-{\mathcal{S}}_{h}p)$. In the next theorem, we derive these results.

Theorem 2.

There exists a positive constant C independent of h such that:

$$\begin{array}{}\mathrm{(56)}& \parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u})\parallel \le C{\displaystyle \sum _{i=1}^{nb}}[\parallel {\displaystyle \frac{\partial \mathit{u}}{\partial t}}{\parallel }_{r,{\Omega }_i}+\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}]h^r,\hfill \mathrm{(57)}& \parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p)\parallel \le C{\displaystyle \sum _{i=1}^{nb}}[\parallel {\displaystyle \frac{\partial \mathit{u}}{\partial t}}{\parallel }_{r,{\Omega }_i}+\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}]h^r.\hfill \end{array}$$

Proof. 

Using the projection operators $P_h,{\mathcal{P}}_h$ and ${\Pi }_0$, the system (51) and (52) can be modified as follows:

$$\begin{array}{ll}(a^{-1}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),\mathit{v})=& {\sum }_{i=1}^{nb}\left\{{(P_{h}p-{\mathcal{S}}_{h}p,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle p-{\mathcal{P}}_{h}p,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\right\}\\ & +(a^{-1}({\Pi }_0\mathit{u}-\mathit{u}),\mathit{v}),\mathit{v}\in {\mathit{S}}_h^0,\hfill \end{array}$$

(58)

$$\begin{array}{c}\hfill \sum _{i=1}^{nb}{(\nabla ·({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),w)}_{{\Omega }_i}=0,w\in U_h.\end{array}$$

(59)

Differentiating (58) and (59) with respect to t, we have the relation:

$$\begin{array}{cc}(a^{-1}\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),\mathit{v})=\hfill & {\sum }_{i=1}^{nb}\left\{{(\frac{\partial }{\partial t}(P_{h}p-{\mathcal{S}}_{h}p),\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle \frac{\partial }{\partial t}(p-{\mathcal{P}}_{h}p),\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\right\}\\ & +(a^{-1}\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-\mathit{u}),\mathit{v}),\mathit{v}\in {\mathit{S}}_h^0,\hfill \end{array}$$

(60)

$$\begin{array}{c}\sum _{i=1}^{nb}{(\nabla ·\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),w)}_{{\Omega }_i}=0,w\in U_h.\hfill \end{array}$$

(61)

Taking $\mathit{v}={\displaystyle \frac{\partial }{\partial t}}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),w={\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p)$ into (60) and (61) and the summing the resultant equations give:

$$\begin{array}{cc}(a^{-1}\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}))\hfill & =(a^{-1}\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-\mathit{u}),\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}))\hfill \\ & +{\sum }_{i=1}^{nb}{\langle \frac{\partial }{\partial t}({\mathcal{P}}_{h}p-p),\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u})·{\mathit{n}}_i\rangle }_{{\Gamma }_i},\hfill \end{array}$$

(62)

which implies:

$$\begin{array}{cc}(a^{-1}\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}),\frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u}))\le \hfill & C[\parallel \frac{\partial }{\partial t}({\Pi }_0\mathit{u}-\mathit{u})\parallel \parallel \frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u})\parallel \hfill \\ & +{\sum }_{i=1}^{nb}\parallel \frac{\partial }{\partial t}({\mathcal{P}}_{h}p-p){\parallel }_{{\Gamma }_i}{\parallel \frac{\partial }{\partial t}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u})·{\mathit{n}}_i\parallel }_{{\Gamma }_i}].\hfill \end{array}$$

(63)

Using (23)–(25), (36), and the fact that ${\displaystyle \frac{\partial }{\partial t}}\left({\mathcal{P}}_{h}p\right)={\mathcal{P}}_h({\displaystyle \frac{\partial p}{\partial t}}),{\displaystyle \frac{\partial }{\partial t}}\left({\Pi }_0\mathit{u}\right)={\Pi }_0({\displaystyle \frac{\partial \mathit{u}}{\partial t}})$, we have:

$$\begin{array}{c}\hfill \parallel {\displaystyle \frac{\partial }{\partial t}}({\Pi }_0\mathit{u}-{\mathcal{R}}_h\mathit{u})\parallel \le C{\displaystyle \sum _{i=1}^{nb}}[\parallel {\displaystyle \frac{\partial \mathit{u}}{\partial t}}{\parallel }_{r,{\Omega }_i}+\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}]h^r\end{array}$$

(64)

Using the triangle inequality, (36) and (64), we have:

$$\begin{array}{c}\hfill \parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u})\parallel \le C\sum _{i=1}^{nb}[\parallel {\displaystyle \frac{\partial \mathit{u}}{\partial t}}{\parallel }_{r,{\Omega }_i}+\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}]h^r.\end{array}$$

(65)

Next, we derive the bound for ${\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p)$. For this, we shall use the duality argument. Let φ be a solution of:

$$\begin{array}{}\mathrm{(66)}& \hfill -\nabla ·(a\nabla \phi )& =-{\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p),\hfill & \mathrm{in}\Omega \hfill \mathrm{(67)}& \hfill \phi & =0,\hfill & \mathrm{on}\partial \Omega ,\hfill \end{array}$$

then by elliptic regularity, we have:

$${∥\phi ∥}_2\le C\parallel {\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p)\parallel .$$

(68)

Differentiating (51) with respect to t, we have:

$$\begin{array}{cc}\hfill (a^{-1}{\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),\mathit{v})& =\sum _{i=1}^{nb}\left\{{({\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p),\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\right\}.\hfill \end{array}$$

(69)

Let $\mathit{v}={\Pi }_{0}a\nabla \phi$, then:

$$\begin{array}{cc}\parallel {\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p){\parallel }^2\hfill & ={\sum }_{i=1}^{nb}{({\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p),\nabla ·{\Pi }_{0}a\nabla \phi )}_{{\Omega }_i}\hfill \\ & ={\sum }_{i=1}^{nb}\{{(a^{-1}{\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),{\Pi }_{0}a\nabla \phi )}_{{\Omega }_i}\hfill \\ & +{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},{\Pi }_{0}a\nabla \phi ·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\}.\hfill \end{array}$$

(70)

Using (36), the first term on the right-hand side of (70) is estimated as:

$$\begin{array}{cc}{\sum }_{i=1}^{nb}{(a^{-1}{\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),{\Pi }_{0}a\nabla \phi )}_{{\Omega }_i}=\hfill & {\sum }_{i=1}^{nb}\{{(a^{-1}{\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),{\Pi }_{0}a\nabla \phi -a\nabla \phi )}_{{\Omega }_i}\hfill \\ & +{(a^{-1}{\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}),a\nabla \phi )}_{{\Omega }_i}\}\hfill \\ & \le C{\sum }_{i=1}^{nb}[\parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h{\mathit{u}\parallel }_{{\Omega }_i}{\parallel {\Pi }_{0}a\nabla \phi -a\nabla \phi \parallel }_{{\Omega }_i}\hfill \\ & +\parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}){\parallel }_{{\Omega }_i}{\parallel a\nabla \phi \parallel }_{{\Omega }_i}]\hfill \\ & \le C{\sum }_{i=1}^{nb}[\parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}){\parallel }_{{\Omega }_i}{\parallel \phi \parallel }_{2,{\Omega }_i}].\hfill \end{array}$$

(71)

For the second term on the right-hand side of (70), we use projection operator ${\Pi }_i$ to obtain:

$$\begin{array}{cc}{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},{\Pi }_{0}a\nabla \phi ·{\mathit{n}}_i\rangle }_{{\Gamma }_i}=\hfill & {\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},({\Pi }_{0}a\nabla \phi -{\Pi }_{i}a\nabla \phi )·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\hfill \\ & +{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},({\Pi }_{i}a\nabla \phi -a\nabla \phi )·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\hfill \\ & +{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},a\nabla \phi ·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\hfill \\ & \le {\sum }_j\parallel {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{{\Gamma }_{ij}}[\parallel ({\Pi }_{0}a\nabla \phi -{\Pi }_{i}a\nabla \phi ).{\mathit{n}}_i{\parallel }_{{\Gamma }_{ij}}\hfill \\ & +\parallel ({\Pi }_{i}a\nabla \phi -a\nabla \phi )·{\mathit{n}}_i{\parallel }_{{\Gamma }_{ij}}]\hfill \\ & +{\sum }_j\parallel {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{-1/2,{\Gamma }_{ij}}{\parallel a\nabla \phi ·{\mathit{n}}_i\parallel }_{1/2,{\Gamma }_{ij}}\hfill \end{array}$$

(72)

Next, by using (10), (23), (24), (25), and (35), terms on the right-hand side of (72) are bounded as follows:

$$\begin{array}{}\mathrm{(73)}& \parallel {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{{\Gamma }_{ij}}\le {\parallel {\displaystyle \frac{\partial p}{\partial t}}\parallel }_{r+1,{\Omega }_i}{h}^{r+1/2},\hfill \mathrm{(74)}& \parallel ({\Pi }_{0}a\nabla \phi -{\Pi }_{i}a\nabla \phi )·{\mathit{n}}_i{\parallel }_{{\Gamma }_{ij}}\le C{h}^{1/2}{\parallel \phi \parallel }_{2,{\Omega }_i},\hfill \mathrm{(75)}& \parallel ({\Pi }_i\nabla \phi -a\nabla \phi )·{\mathit{n}}_i{\parallel }_{{\Gamma }_{ij}}\le C{h}^{1/2}{\parallel \phi \parallel }_{2,{\Omega }_i},\hfill \mathrm{(76)}& \parallel {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{-1/2,{\Gamma }_{ij}}\le C{\parallel {\displaystyle \frac{\partial p}{\partial t}}\parallel }_{r+1,{\Omega }_i}{h}^{r+1},\hfill \mathrm{(77)}& \parallel a\nabla \phi ·{\mathit{n}}_i{\parallel }_{1/2,{\Gamma }_{ij}}\le {\parallel \phi \parallel }_{2,{\Omega }_i}.\hfill \end{array}$$

Combining (73)–(77) with (72), we arrive at:

$$\begin{array}{c}{\langle {\displaystyle \frac{\partial p}{\partial t}}-{\mathcal{P}}_h{\displaystyle \frac{\partial p}{\partial t}},{\Pi }_{0}a\nabla \phi ·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\le \parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}{h}^{r+1}{\parallel \phi \parallel }_{2,{\Omega }_i}.\hfill \end{array}$$

(78)

Substituting (71) and (78) in (70), we obtain:

$$\parallel {\displaystyle \frac{\partial }{\partial t}}(P_{h}p-{\mathcal{S}}_{h}p){\parallel }_{\Omega }\le C{\displaystyle \sum _{i=1}^{nb}}[\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}{h}^{r+1}+\parallel {\displaystyle \frac{\partial }{\partial t}}(\mathit{u}-{\mathcal{R}}_h\mathit{u}){\parallel }_{{\Omega }_i}].$$

(79)

Using the triangle inequality (19), (65), and (79), we have:

$$\begin{array}{c}\hfill \parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p)\parallel \le C{\displaystyle \sum _{i=1}^{nb}}[\parallel {\displaystyle \frac{\partial \mathit{u}}{\partial t}}{\parallel }_{r,{\Omega }_i}+\parallel {\displaystyle \frac{\partial p}{\partial t}}{\parallel }_{r+1,{\Omega }_i}]h^r.\end{array}$$

(80)

which is required. □

4. A Priori Error Estimates for the Semidiscrete Scheme

In this section, we shall derive the error estimate for velocity $\mathit{u}$ and pressure p for our semidiscrete scheme.

Theorem 3.

There exists a positive constant C independent of h such that:

$$\begin{array}{}\mathrm{(81)}& \parallel \mathit{u}-{\mathit{u}}_h{\parallel }_{L^2(J:L^2(\Omega ))}\le C[\parallel \mathit{u}-{\mathcal{R}}_h{\mathit{u}\parallel }_{L^2(J:L^2(\Omega ))}+\parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p){\parallel }_{L^2(J:L^2(\Omega ))}].\hfill \mathrm{(82)}& \parallel p-p_h{\parallel }_{L^{\infty }(J:L^2(\Omega ))}\le C[\parallel \mathit{u}-{\mathcal{R}}_h{\mathit{u}\parallel }_{L^2(J:L^2(\Omega ))}+\parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p){\parallel }_{L^2(J:L^2(\Omega ))}].\hfill \end{array}$$

Proof. 

Observe that from the triangle inequality, we have:

$$\begin{array}{}\mathrm{(83)}& \parallel \mathit{u}-{\mathit{u}}_h\parallel \le \parallel \mathit{u}-{\mathcal{R}}_h\mathit{u}\parallel +\parallel {\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h\parallel .\hfill \mathrm{(84)}& \parallel p-p_h\parallel \le \parallel p-{\mathcal{S}}_{h}p\parallel +\parallel {\mathcal{S}}_{h}p-p_h\parallel .\hfill \end{array}$$

The terms $\parallel \mathit{u}-{\mathcal{R}}_h\mathit{u}\parallel$ and $\parallel p-{\mathcal{S}}_{h}p\parallel$ in Inequalities (83) and (84) are estimated by (54) and (55), respectively. Therefore, it remains to bound the terms $\parallel {\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h\parallel$ and $\parallel {\mathcal{S}}_{h}p-p_h\parallel .$ Subtracting (30)–(32) from (45)–(47), we obtain:

$$\begin{array}{}\mathrm{(85)}& {(a^{-1}({\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h),\mathit{v})}_{{\Omega }_i}={({\mathcal{S}}_{h}p-p_h,\nabla ·\mathit{v})}_{{\Omega }_i}-{\langle {\mathcal{L}}_h\lambda -{\lambda }_h,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i},\hfill & \hfill \mathit{v}\in {\mathit{S}}_{h,i},\mathrm{(86)}& {({\displaystyle \frac{\partial {\mathcal{S}}_{h}p}{\partial t}}-{\displaystyle \frac{\partial p_h}{\partial t}},w)}_{{\Omega }_i}+{(\nabla ·({\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h),w)}_{{\Omega }_i}=-{({\displaystyle \frac{\partial p}{\partial t}}-{\displaystyle \frac{\partial {\mathcal{S}}_{h}p}{\partial t}},w)}_{{\Omega }_i},\hfill & \hfill w\in U_{h,i},\mathrm{(87)}& {\displaystyle \sum _{i=1}^{nb}}{\langle ({\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h)·{\mathit{n}}_i,\mu \rangle }_{{\Gamma }_i}=0,\hfill & \hfill \mu \in {\Sigma }_h.\end{array}$$

Taking $\mathit{v}={\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h,w={\mathcal{S}}_{h}p-p_h$, $\mu ={\mathcal{L}}_h\lambda -{\lambda }_h$ in (85)–(87) and adding the resultant equations, we obtain:

$$\begin{array}{cc}{(a^{-1}({\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h),({\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h))}_{{\Omega }_i}+\hfill & {({\displaystyle \frac{\partial {\mathcal{S}}_{h}p}{\partial t}}-{\displaystyle \frac{\partial p_h}{\partial t}},{\mathcal{S}}_{h}p-p_h)}_{{\Omega }_i}=\hfill \\ & -{({\displaystyle \frac{\partial p}{\partial t}}-{\displaystyle \frac{\partial {\mathcal{S}}_{h}p}{\partial t}},{\mathcal{S}}_{h}p-p_h)}_{{\Omega }_i}.\hfill \end{array}$$

(88)

Clearly:

$$\begin{array}{cc}{({\displaystyle \frac{\partial {\mathcal{S}}_{h}p}{\partial t}}-{\displaystyle \frac{\partial p_h}{\partial t}},{\mathcal{S}}_{h}p-p_h)}_{{\Omega }_i}\hfill & ={({\displaystyle \frac{\partial }{\partial t}}({\mathcal{S}}_{h}p-p_h),{\mathcal{S}}_{h}p-p_h)}_{{\Omega }_i}\hfill \\ & ={\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\parallel {\mathcal{S}}_{h}p-p_h\parallel }_{{\Omega }_i}^2.\hfill \end{array}$$

(89)

Using (5), (89), the Cauchy Schwarz inequality, and $ϵ$ inequality, we have from (88):

$$\begin{array}{cc}\hfill {\displaystyle \frac{1}{c_0}}\parallel {\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h{\parallel }_{{\Omega }_i}^2+{\displaystyle \frac{1}{2}}{\displaystyle \frac{d}{dt}}{\parallel {\mathcal{S}}_{h}p-p_h\parallel }_{{\Omega }_i}^2& \le C[\parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p){\parallel }_{{\Omega }_i}\parallel {\mathcal{S}}_{h}p-p_h{\parallel }_{{\Omega }_i}]\hfill \\ & \le C[\parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p){\parallel }_{{\Omega }_i}^2+\parallel {\mathcal{S}}_{h}p-p_h{\parallel }_{{\Omega }_i}^2].\hfill \end{array}$$

Summing up over subdomains ${\Omega }_i$, integrating with respect to t, and using Gronwall’s inequality, we obtain:

$$\begin{array}{c}\parallel {\mathcal{R}}_h\mathit{u}-{\mathit{u}}_h{\parallel }_{L^2(J:L^2(\Omega ))}+{\parallel {\mathcal{S}}_{h}p-p_h\parallel }_{L^{\infty }(J:L^2(\Omega ))}\hfill \\ \le C[\parallel \mathit{u}-{\mathcal{R}}_h\mathit{u}{\parallel }_{L^2(J:L^2(\Omega ))}+\parallel {\displaystyle \frac{\partial }{\partial t}}(p-{\mathcal{S}}_{h}p){\parallel }_{L^2(J:L^2(\Omega ))}].\hfill \end{array}$$

(90)

Combining (83) and (84) and (88) and (90), we conclude (81) and (82). □

5. An Interface Error Estimate

This section is devoted to deriving an error estimate for the Lagrange multiplier introduced on the interface. The interface error estimation reflects the accuracy of data transmitted across the subdomain interfaces. For this, subtract (30) from (26) to get:

$$\begin{array}{cc}\hfill {\langle \lambda -{\lambda }_h,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}={(p-p_h,\nabla ·\mathit{v})}_{{\Omega }_i}-{(a^{-1}(\mathit{u}-{\mathit{u}}_h),\mathit{v})}_{{\Omega }_i},& \mathit{v}\in {\mathit{S}}_{h,i}.\hfill \end{array}$$

(91)

For $\mu \in L^2(\Gamma )$, let us define the norm:

$${\left|\left|\left|\mu \right|\right|\right|}_{\Gamma }:=\underset{\begin{array}{c}\mathit{v}\in {\mathit{S}}_h\\ \mathit{v}\ne 0\end{array}}{sup}\frac{{\sum }_{i=1}^{nb}{\langle \mu ,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}}{{\parallel \mathit{v}\parallel }_{\mathit{S}}}.$$

(92)

and:

$${\left|\left|\left|\mu \right|\right|\right|}_{J\times \Gamma }:={\left({\int }_J{\left|\left|\left|\mu (.,t)\right|\right|\right|}_{\Gamma }^2\right)}^{1/2}.$$

(93)

The relation (91) implies:

$$\begin{array}{cc}{\langle \lambda -{\lambda }_h,\mathit{v}·{\mathit{n}}_i\rangle }_{{\Gamma }_i}\hfill & \le C[\parallel p-p_h{\parallel }_{{\Omega }_i}\parallel \nabla ·\mathit{v}{\parallel }_{{\Omega }_i}+\parallel \mathit{u}-{\mathit{u}}_h{\parallel }_{{\Omega }_i}\parallel \mathit{v}{\parallel }_{{\Omega }_i}],\hfill \\ & \le C[\parallel p-p_h{\parallel }_{{\Omega }_i}\parallel \mathit{v}{\parallel }_{\mathit{H}(div;{\Omega }_i)}+\parallel \mathit{u}-{\mathit{u}}_h{\parallel }_{{\Omega }_i}\parallel \mathit{v}{\parallel }_{{\Omega }_i}],\hfill \end{array}$$

(94)

then by using (92), we see that:

$$\begin{array}{cc}\hfill {\left|\left|\left|\lambda -{\lambda }_h\right|\right|\right|}_{\Gamma }& \le C{\left(\sum _{i=1}^{nb}(\parallel p-p_h{\parallel }_{{\Omega }_i}+\parallel \mathit{u}-{\mathit{u}}_h{{\parallel }_{{\Omega }_i})}^2\right)}^{1/2}\hfill \\ & \le C\left[\parallel p-p_h\parallel +\parallel \mathit{u}-{\mathit{u}}_h\parallel \right].\hfill \end{array}$$

Integrating with respect to time, we obtain:

$${\left|\left|\left|\lambda -{\lambda }_h\right|\right|\right|}_{J\times \Gamma }\le C\left[\parallel p-p_h{\parallel }_{L^2(J;L^2(\Omega ))}+{\parallel \mathit{u}-{\mathit{u}}_h\parallel }_{L^2(J;L^2(\Omega ))}\right].$$

We have proven the following theorem.

Theorem 4.

There exists a positive constant C independent of h such that:

$${\left|\left|\left|\lambda -{\lambda }_h\right|\right|\right|}_{J\times \Gamma }\le C\left[\parallel p-p_h{\parallel }_{L^2(J;L^2(\Omega ))}+{\parallel \mathit{u}-{\mathit{u}}_h\parallel }_{L^2(J;L^2(\Omega ))}\right].$$

Note that $\parallel p-p_h{\parallel }_{L^2(J;L^2(\Omega ))}\parallel \mathit{u}-{\mathit{u}}_h{\parallel }_{L^2(J;L^2(\Omega ))}$ are estimated in Theorem 3. Thus, Theorem 4 gives an error bound for the interface pressure variable.

6. Numerical Results

This section presents the results of numerical examples to illustrate the theoretical findings. For the sake of simplicity, we chose the unit square $\Omega ={[0,1]}^2$ as the spatial domain and the time interval $J=[0,1]$. We constructed two examples to verify the convergence order of the scheme and considered two cases for each example. In the first case, we considered two subdomains ($nb=2$). The second case consisted of four subdomains $(nb=4)$. On four subdomains, the interface was along $x=1/2,y=1/2$, while for the two subdomains, it was at $y=1/2$. We used the lowest order Raviart–Thomas $\left(R{T}_0\right)$ space on each subdomain, which approximates the unknowns $({\mathit{u}}_h,p_h)$ in $(P_0,R{T}_0)$. The piecewise continuous linear polynomials were used to approximate the pressure variable on the mortar interface between the subdomains. The error in mortar interface space was calculated by using the discrete $L^2$-norm. The errors were reported at time $t=1$ for each case. We also present the discrete errors for each time step from $t=0$ to $t=1$.

Example 1.

Consider $\Omega ={[0,1]}^2$, $a=\frac{1}{4}(x^2+y^2+1)$, the source term f, and the velocity $\mathit{u}$ are calculated by using the exact solution:

$$p(x,y)=txy(1-x)(1-y).$$

(95)

The results of numerical experiment for fixed h and varying time steps for Example 1 are displayed in

Table 1. Discrete errors with time steps for Example 1, $nb=2,h=1/64,\delta t=0.1$.

${\mathit{t}}^{\mathit{n}}$	$\parallel {\mathit{u}}^{\mathit{n}}-{\mathit{u}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{p}}^{\mathit{n}}-{\mathit{p}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{\lambda }}^{\mathit{n}}-{\mathit{\lambda }}_{\mathit{h}}^{\mathit{n}}\parallel$
0.1	1.5833 × 10${}^{-4}$	7.7669 × 10${}^{-5}$	2.7748 × 10${}^{-6}$
0.2	3.1668 × 10${}^{-4}$	1.5534 × 10${}^{-4}$	5.6298 × 10${}^{-6}$
0.3	4.7503 × 10${}^{-4}$	2.3302 × 10${}^{-4}$	8.5262 × 10${}^{-6}$
0.4	6.3338 × 10${}^{-4}$	3.1070 × 10${}^{-4}$	1.1446 × 10${}^{-5}$
0.5	7.9173 × 10${}^{-4}$	3.8838 × 10${}^{-4}$	1.4380 × 10${}^{-5}$
0.6	9.5008 × 10${}^{-4}$	4.6606 × 10${}^{-4}$	1.7322 × 10${}^{-5}$
0.7	1.1084 × 10${}^{-3}$	5.4374 × 10${}^{-4}$	2.0268 × 10${}^{-5}$
0.8	1.2668 × 10${}^{-3}$	6.2142 × 10${}^{-4}$	2.3217 × 10${}^{-5}$
0.9	1.4251 × 10${}^{-3}$	6.9910 × 10${}^{-4}$	2.6167 × 10${}^{-5}$
1.0	1.5835 × 10${}^{-3}$	7.7678 × 10${}^{-4}$	2.9119 × 10${}^{-5}$

and

Table 2. Discrete errors with time steps for Example 1, $nb=4,h=1/64,\delta t=0.1$.

${\mathit{t}}^{\mathit{n}}$	$\parallel {\mathit{u}}^{\mathit{n}}-{\mathit{u}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{p}}^{\mathit{n}}-{\mathit{p}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{\lambda }}^{\mathit{n}}-{\mathit{\lambda }}_{\mathit{h}}^{\mathit{n}}\parallel$
0.1	2.1943 × 10${}^{-4}$	1.0941 × 10${}^{-4}$	3.9312 × 10${}^{-6}$
0.2	4.3889 × 10${}^{-4}$	2.1884 × 10${}^{-4}$	7.9759 × 10${}^{-6}$
0.3	6.5835 × 10${}^{-4}$	3.2827 × 10${}^{-4}$	1.2080 × 10${}^{-5}$
0.4	8.7782 × 10${}^{-4}$	4.3769 × 10${}^{-4}$	1.6217 × 10${}^{-5}$
0.5	1.0973 × 10${}^{-3}$	5.4712 × 10${}^{-4}$	2.0374 × 10${}^{-5}$
0.6	1.3167 × 10${}^{-3}$	6.5656 × 10${}^{-4}$	2.4541 × 10${}^{-5}$
0.7	1.5362 × 10${}^{-3}$	7.6599 × 10${}^{-4}$	2.8716 × 10${}^{-5}$
0.8	1.7557 × 10${}^{-3}$	8.7542 × 10${}^{-4}$	3.2893 × 10${}^{-5}$
0.9	1.9751 × 10${}^{-3}$	9.8485 × 10${}^{-4}$	3.7073 × 10${}^{-5}$
1.0	2.1946 × 10${}^{-3}$	1.0943 × 10${}^{-3}$	4.1255 × 10${}^{-5}$

. We choose the fixed $h=1/64$ along with $\delta t=1/10$, the time step. We report the errors on two and four subdomains.

The discrete norm errors and convergence rates on two subdomains for Example 1 are given in

Table 3. Discrete errors and convergence rates for Example 1 at time $T=1$ with $\delta t=0.1$ and $nb=2$.

h	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{p}-{\mathit{p}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{\lambda }-{\mathit{\lambda }}_{\mathit{h}}\parallel$	CR
1/4	2.6484 × 10${}^{-2}$		1.3578 × 10${}^{-2}$		6.3437 × 10${}^{-3}$
1/8	1.2845 × 10${}^{-2}$	1.0439	6.3866 × 10${}^{-3}$	1.0882	1.7877 × 10${}^{-3}$	1.8272
1/16	6.3564 × 10${}^{-3}$	1.0149	3.1288 × 10${}^{-3}$	1.0294	4.6132 × 10${}^{-4}$	1.9543
1/32	3.1693 × 10${}^{-3}$	1.0040	1.5558 × 10${}^{-3}$	1.0080	1.1625 × 10${}^{-4}$	1.9885
1/64	1.5835 × 10${}^{-3}$	1.0010	7.7678 × 10${}^{-4}$	1.0021	2.9119 × 10${}^{-5}$	1.9972
1/128	7.9160 × 10${}^{-4}$	1.0003	3.8825 × 10${}^{-4}$	1.0005	7.2832 × 10${}^{-6}$	1.9993
1/256	3.9578 × 10${}^{-4}$	1.0001	1.9411 × 10${}^{-4}$	1.0001	1.8210 × 10${}^{-6}$	1.9998

. We observe the convergence rate of order $\mathcal{O}\left(h\right)$ for both subdomain variables, that is pressure p and velocity $\mathit{u}$. The convergence order for pressure on the mortar interface was of order $\mathcal{O}\left(h^2\right)$.

The errors and convergence rates on four subdomains $(nb=4)$ are provided in

Table 4. Discrete errors and convergence rates for Example 1 at time $T=1$ with $\delta t=0.1$ and $nb=4$.

h	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{p}-{\mathit{p}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{\lambda }-{\mathit{\lambda }}_{\mathit{h}}\parallel$	CR
1/4	3.6671 × 10${}^{-2}$		1.9202 × 10${}^{-2}$		9.3169 × 10${}^{-3}$
1/8	1.7799 × 10${}^{-2}$	1.0428	9.0122 × 10${}^{-3}$	1.0913	2.5699 × 10${}^{-3}$	1.8581
1/16	8.8092 × 10${}^{-3}$	1.0147	4.4097 × 10${}^{-3}$	1.0312	6.5739 × 10${}^{-4}$	1.9669
1/32	4.3924 × 10${}^{-3}$	1.0040	2.1919 × 10${}^{-3}$	1.0085	1.6500 × 10${}^{-4}$	1.9943
1/64	2.1946 × 10${}^{-3}$	1.0011	1.0943 × 10${}^{-3}$	1.0022	4.1255 × 10${}^{-5}$	1.9998
1/128	1.0971 × 10${}^{-3}$	1.0003	5.4693 × 10${}^{-4}$	1.0006	1.0309 × 10${}^{-5}$	2.0007
1/256	5.4853 × 10${}^{-4}$	1.0001	2.7344 × 10${}^{-4}$	1.0001	2.5765 × 10${}^{-6}$	2.0004

. The optimal order convergence rates of order $\mathcal{O}\left(h\right)$ were obtained for the subdomains pressure and velocity. The convergence order for the pressure $\lambda$ on interface is of order $\mathcal{O}\left(h^2\right)$.

The exact and computed solution for Example 1 are shown in Figure 1.

Example 2.

Consider $\Omega ={[0,1]}^2$, $a=e^{x+y}$, the source term f, and the velocity $\mathit{u}$ are calculated by using the exact solution:

$$p(x,y)=tsin\left(\pi x\right)sin\left(\pi y\right).$$

(96)

The discrete norm errors with respect to the variable temporal steps and fixed h for two subdomains are presented in

Table 5. Discrete errors with time steps for Example 2, $nb=2,h=1/64,\delta t=0.1$.

${\mathit{t}}^{\mathit{n}}$	$\parallel {\mathit{u}}^{\mathit{n}}-{\mathit{u}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{p}}^{\mathit{n}}-{\mathit{p}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{\lambda }}^{\mathit{n}}-{\mathit{\lambda }}_{\mathit{h}}^{\mathit{n}}\parallel$
0.1	1.5684 × 10${}^{-2}$	1.1572 × 10${}^{-3}$	1.8226 × 10${}^{-5}$
0.2	3.1369 × 10${}^{-2}$	2.3143 × 10${}^{-3}$	3.6545 × 10${}^{-5}$
0.3	4.7053 × 10${}^{-2}$	3.4715 × 10${}^{-3}$	5.4869 × 10${}^{-5}$
0.4	6.2738 × 10${}^{-2}$	4.6287 × 10${}^{-3}$	7.3194 × 10${}^{-5}$
0.5	7.8422 × 10${}^{-2}$	5.7859 × 10${}^{-3}$	9.1519 × 10${}^{-5}$
0.6	9.4107 × 10${}^{-2}$	6.9430 × 10${}^{-3}$	1.0984 × 10${}^{-4}$
0.7	1.0979 × 10${}^{-1}$	8.1002 × 10${}^{-3}$	1.2817 × 10${}^{-4}$
0.8	1.2548 × 10${}^{-1}$	9.2574 × 10${}^{-3}$	1.4649 × 10${}^{-4}$
0.9	1.4116 × 10${}^{-1}$	1.0415 × 10${}^{-2}$	1.6482 × 10${}^{-4}$
1.0	1.5684 × 10${}^{-1}$	1.1572 × 10${}^{-2}$	1.8314 × 10${}^{-4}$

and

Table 6. Discrete errors with time steps for Example 2, $nb=4,h=1/64,\delta t=0.1$.

${\mathit{t}}^{\mathit{n}}$	$\parallel {\mathit{u}}^{\mathit{n}}-{\mathit{u}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{p}}^{\mathit{n}}-{\mathit{p}}_{\mathit{h}}^{\mathit{n}}\parallel$	$\parallel {\mathit{\lambda }}^{\mathit{n}}-{\mathit{\lambda }}_{\mathit{h}}^{\mathit{n}}\parallel$
0.1	2.1580 × 10${}^{-2}$	1.6280 × 10${}^{-3}$	2.5642 × 10${}^{-5}$
0.2	4.3160 × 10${}^{-2}$	3.2560 × 10${}^{-3}$	5.1413 × 10${}^{-5}$
0.3	6.4741 × 10${}^{-2}$	4.8840 × 10${}^{-3}$	7.7192 × 10${}^{-5}$
0.4	8.6321 × 10${}^{-2}$	6.5120 × 10${}^{-3}$	1.0297 × 10${}^{-4}$
0.5	1.0790 × 10${}^{-1}$	8.1400 × 10${}^{-3}$	1.2875 × 10${}^{-4}$
0.6	1.2948 × 10${}^{-1}$	9.7680 × 10${}^{-3}$	1.5453 × 10${}^{-4}$
0.7	1.5106 × 10${}^{-1}$	1.1396 × 10${}^{-2}$	1.8031 × 10${}^{-4}$
0.8	1.7264 × 10${}^{-1}$	1.3024 × 10${}^{-2}$	2.0609 × 10${}^{-4}$
0.9	1.9422 × 10${}^{-1}$	1.4652 × 10${}^{-2}$	2.3187 × 10${}^{-4}$
1.0	2.1580 × 10${}^{-1}$	1.6280 × 10${}^{-2}$	2.5765 × 10${}^{-4}$

. The tables shows that the errors increased with time steps.

Next, the discrete errors for subdomain pressure p, the subdomain velocity $\mathit{u}$, and the mortar interface pressure $\lambda$ are presented in

Table 7. Discrete errors and convergence rates for Example 2 at time $T=1$ with $\delta t=0.1$ and $nb=2$.

h	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{p}-{\mathit{p}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{\lambda }-{\mathit{\lambda }}_{\mathit{h}}\parallel$	CR
1/4	$2.6702\times {10}^0$		1.9026 × 10${}^{-1}$		4.0137 × 10${}^{-2}$
1/8	$1.2791\times {10}^0$	1.0618	$9.3369\times {10}^{-2}$	1.0270	1.1145 × 10${}^{-2}$	1.8485
1/16	$6.3045\times {10}^{-1}$	1.0207	4.6387 × 10${}^{-2}$	1.0092	2.8930 × 10${}^{-3}$	1.9458
1/32	3.1400 × 10${}^{-1}$	1.0056	2.3154 × 10${}^{-2}$	1.0025	7.3067 × 10${}^{-4}$	1.9853
1/64	1.5684 × 10${}^{-1}$	1.0015	1.1572 × 10${}^{-2}$	1.0006	1.8314 × 10${}^{-4}$	1.9963
1/128	7.8403 × 10${}^{-2}$	1.0003	5.7852 × 10${}^{-3}$	1.0002	4.5816 × 10${}^{-5}$	1.9990
1/256	3.9199 × 10${}^{-2}$	1.0001	2.8925 × 10${}^{-3}$	1.0000	1.1456 × 10${}^{-5}$	1.9990

and

Table 8. Discrete errors and convergence rates for Example 2 at time $T=1$ with $\delta t=0.1$ and $nb=4$.

h	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{p}-{\mathit{p}}_{\mathit{h}}\parallel$	CR	$\parallel \mathit{\lambda }-{\mathit{\lambda }}_{\mathit{h}}\parallel$	CR
1/4	3.6498 × 10${}^0$		2.6843 × 10${}^{-1}$		4.8133 × 10${}^{-2}$
1/8	1.7572 × 10${}^0$	1.0545	1.3147 × 10${}^{-1}$	1.0298	1.4917 × 10${}^{-2}$	1.6901
1/16	8.6713 × 10${}^{-1}$	1.0190	6.5276 × 10${}^{-2}$	1.0101	3.9978 × 10${}^{-3}$	1.8997
1/32	4.3200 × 10${}^{-1}$	1.0052	3.2576 × 10${}^{-2}$	1.0027	1.0223 × 10${}^{-3}$	1.9674
1/64	2.1580 × 10${}^{-1}$	1.0013	1.6280 × 10${}^{-2}$	1.0007	2.5765 × 10${}^{-4}$	1.9883
1/128	1.0788 × 10${}^{-1}$	1.0003	8.1390 × 10${}^{-3}$	1.0002	6.4627 × 10${}^{-5}$	1.9952
1/256	5.3935 × 10${}^{-2}$	1.0001	4.0694 × 10${}^{-3}$	1.0000	1.6180 × 10${}^{-5}$	1.9979

.

The discrete errors and convergence rates for Example 2 on four subdomains are represented in Table 8.

The results of our numerical experiments matched the theory established. The exact and computed solution for Example 2 are presented in Figure 2.

It is observed from the above tables that the discrete norm errors increased with the time steps. Note that the accuracy of discrete norm error on two subdomain was better than that on four subdomains. The reason for the better accuracy on the two subdomains is stated below in the discussion of the convergence rate.

Comparing the results on two and four subdomains, we note the slightly better accuracy on the two subdomains than the four subdomains. The reason for the better accuracy on the two subdomains may be that in the case of two subdomains, we had a single interface; therefore, transmission of data across the interface was more accurate as compared to that of four subdomains. However, on two subdomains, we had large local problem, which is expensive to solve numerically. Therefore, by scarifying accuracy slightly, we have an efficient solution for the local subdomain problems.

7. Conclusions

In this article, the multiblock mortar mixed technique proposed and analyzed in [20] was extended to the second order linear parabolic problems. This approach is a combination of non-overlapping domain decomposition and the mortar finite element method. This method allowed different grids on different regions of the domain. The continuity of the flux across the subdomain interior boundary was imposed via the mortar finite element on the interface grid.

We considered the semi-discrete formulation of second order parabolic partial differential equation. The existence and uniqueness of solution for the discrete system encountered by the scheme were established. The optimal order convergence was shown for both the velocity and pressure approximations. The theoretical findings were confirmed by performing the numerical experiments for benchmark problems.