Weiwei Wang,Chuanju Xu

School of Mathematical Sciences,Xiamen University,361005 Xiamen,China

Received 17 January 2014;Accepted(in revised version)21 February 2014

A Posteriori Error Estimation of Spectral and Spectral Element Methods for the Stokes/Darcy Coupled Problem

Weiwei Wang,Chuanju Xu?

School of Mathematical Sciences,Xiamen University,361005 Xiamen,China

Received 17 January 2014;Accepted(in revised version)21 February 2014

.In this paper,we carry out an a posteriori error analysis of Legendre spectral approximations to the Stokes/Darcy coupled equations.The spectral approximations are based on a weak formulation ofthe coupled equations by using the Beavers-Joseph-Saffman interface condition.The main contribution of the paper consists of deriving a number of posteriori error indicators and their upper and lower bounds for the single domain case.An extension of the upper bounds to the multi-domain case in the spectral element framework is also given.

AMS subject classifications:35Q35,76D03

Aposteriori error,Stokes/Darcy coupled equations,Spectralmethod,Spectralelement method.

1 Introduction

The model of the Stokes equations coupled with the Darcy equations has been a subject of interest in a large variety of different fields,see,e.g.[9-12,14,16].Recently,we have introduced a new formulation for the Stokes/Darcy coupled equations,subject respectively to the Beavers-Joseph-Saffman interface condition and an alternative matching interface condition[22].Some spectral approximations are proposed and a priori error estimates are derived therein.In this paper we consider a posteriori error analysis for the above mentioned spectral approximations.The motivation of this consideration is that a posteriori error estimators are computable quantities in terms of the discrete solution,and can be used to measure the actual approximation errors without the knowledge of exact solutions.They are essential for designing algorithms with adaptive mesh refinement with minimal computational cost.On the other side,there are few work on a posteriori error analysis of the spectral method,and it is not clear if the adaptive strategy in the spectral method can be as efficient as in the finite element framework.Therefore this paper can be regarded as a step towards a better understanding about the adaptive spectral method,with particular attention to the Stokes/Darcy coupled equations.

Some a posteriori erroranalysis ofthe finite elementapproximation to the Stokes Darcy coupled equations have been carried in[1,8].The work[1]used the Lagrangian multiplier in their variational formulation,while[8]replaced the Darcy equations with a Poisson-like equation.In[1]the finite element subspaces consist of Bernardi-Raugel and Raviart-Thomas elements for the velocities,piecewise constants for the pressures,and continuous piecewise linear elements for the Lagrange multiplier defined on the interface.They have derived a residual-based a posteriori error estimate for the Stokes/Darcy coupled problem.The finite element spaces adopted in[8]are the Hood-Taylor element for the velocity and the pressure in the Stokes equation and conforming piecewise quadratic element for the Darcy pressure.The a posteriori error analysis was based on a suitable evaluation of the residual of the finite element solution.

In contrast to the lower order methods,a posteriori error estimation for high order methods such as spectral method is much less developed,although there exist a few papers on this topic for the elliptic problems(see,e.g.,[3,6,13,19]).The purpose of this work is to carry out an a posteriori error analysis for the spectral approximation of the Stokes/Darcy coupled equations.The analysis will be based on the formulation introduced in our previous work,which allows to extend the idea from[1,6,8,19]to derive the residual-based a posteriori error estimator in the framework of spectral element method.

The rest of the paper is organized as follows.In section 2 we briefly recall the formulation proposed in[22]for the Stokes/Darcy coupled problem.The core of the work is given by section 3 and section 4,where we develop the a posteriori error analysis.In section 3 we derive a residual-based a posteriori error estimate.The efficiency of this estimate is given in section 4.In section 5,we extend the results to the case of multi-domain in the framework of the spectral element method.

2 The Stokes/Darcy coupled problem

We are interested in the following Stokes/Darcy coupled equations in two dimensions:

where ui=u|Ωiand pi=p|Ωi,i=s,d,with u and p denoting the velocity and pressure respectively,

ν＞0 is the kinematic viscosity of the fluid,f is a given volumetric force,κis defined by

withκ1andκ2parameters associated to the kinematic viscosity of fluid,the permeability and porosity of the porous media in x and y direction respectively.

In(2.1),the computational domain Ω is assumed to be an open bounded subset of R2,with Lipschitz boundary?Ω.Ωsand Ωdare respectively the fluid and porous media subdomains of Ω,such that Ωs∩Ωd=?,ˉΩs∪ˉΩd=ˉΩ,Γs=?Ωs∩?Ω,Γd=?Ωd∩?Ω;see Fig.1 for an example of such a domain.The unit normal vector nson the boundary Γsis chosen pointing outwards from Ωs(similarly for the notation nd).

Figure 1:The model computational domain of the coupled problem.

Mathematically,it is known that some suitable conditions on the interface Γ:=?Ωs∩?Ωdare needed to close Eq.(2.1).Here we consider the following matching conditions on the interface:

The first condition guarantees that the exchange of fluid between the two domains is conservative.The second one guarantees the balance of two driving forces.The third condition is usually called Beavers-Joseph-Saffman interface condition,which has been derived following the work of[16,17,21,22].

2.1 Weak formulation

To construct the weak formulation of(2.1)-(2.2),we need some basic notations.We use the notations L2,H1,H10,and so on,to mean the usual Sobolev spaces.Let(u,v)Ω=RΩuv,〈u,v〉Γ=RΓuv,L2(Ω)=(L2(Ω))2,H1(Ω)=(H1(Ω))2,and L20(Ω)={q:q∈L2(Ω),(q,1)Ω=0}.We introduce the following functional spaces:

The spaces X and M are respectively equipped with the norms

Proceeding in the usual way(see,e.g.[22]),we find that the variational formulation of(2.1)-(2.2)reads:Find(u,p)∈X×M such that

where a(·,·)and b(·,·)are two bilinear forms,defined respectively by

F:X→R is the linear functional:

Theorem 2.1.The weak problem(2.3)admits a unique solution(u,p)∈X×M.Moreover this solution satisfies:

Proof.The well-posedness of(2.3)follows from a straightforward application of the saddle point theory by verifying that a(·,·)and b(·,·)are continuous,a(·,·)is coercive,and b(·,·)satisfies the following LBB compatibility condition[23]:there existsβ＞0 such that

then we have

Lemma 2.1.There exists a constant c＞0,such that

Proof.The proof is similar to Lemma 2.1 of[8].

Lemma 2.2.Let(u,p)∈X×M be the unique solution of(2.3),then

Proof.Takingvd=0 in the first equation of(2.3),we get

Integrating by parts and rearranging the terms,we obtain that

from which we conclude(2.5).

2.2 Spectral discretizations

We consider the spectral method to approximate coupled problems(2.3).For ease of presentation,let’s first assume that both Ωsand Ωdare rectangular domains.Define two discrete spaces:

where QNis the space of algebraic polynomials of degree less than or equal to N with respect to each single variable x or y.Then the spectral approximation to(2.3)reads:find(uN,pN)∈XN×MN,such that:

where aN(·,·),bN(·,·)are two bilinear forms,defined by:

with(·,·)GL,(·,·)Gbe evaluations of the continuous inner product(·,·)by the Gauss-Lobatto and Gauss quadratures respectively.

Theorem 2.2.The spectral discrete problem(2.6)is well posed.

Proof.Similarly to the proof of Theorem 2.1,this theorem can be proved by directly verifying the continuity and coercivity ofthe form aN(·,·)in XN×XN,the continuity of bN(·,·)in XN×MN,and the following LBB condition(Lemma 3.1 of[22]):there exists a positiveβN,which may be dependent on N,such that

Moreover,if we define the bilinear form AN(·,·)by

Then,as a consequence of the ellipticity property of the bilinear form aN(·,·)and the inf-sup condition on bN(·,·),we have the following result.

Lemma 2.3.It holds

where cNis a positive constant depending onβN.

3 A posteriori error estimation

This section is devoted to deriving a posteriori error estimation for the spectral approximation of the Stokes/Darcy coupled equations.

3.1 Upper bound estimation

Our analysis starts with the upper bound estimation for the error in terms of the error indicator.The analysis makes use of some known results on the polynomial approximation theory.

·Polynomial approximation theory

We first recall some well-known spectral projection operators,which will play an important role in the analysis of the error upper bound.Detailed proofs of the results presented in the following can be found in[2,5,20].

It is known(see[5])that

where Lkm(x)=Lk(x)Lm(y),Lkis the Legendre polynomial of degree k.

where INis the Lagrange interpolation operator.

·Upper bound

In the following theorem,we give the main result of the a posteriori error estimation for the Stokes/Darcy coupled equations.

Theorem 3.1.Assume that the datafbelongs toHr(Ωs),r＞1.Let(u,p)be the solution of(2.3),(uN,pN)be the solution of(2.6).Then it holds

with

and fN=INf.

Proof.Let eu=u-uN,ep=p-pN,it follows from Lemma 2.1 that there exists(v,q)∈X×M such that

Next,we estimate the terms on the right-hand side.For the last term,we have

Using(3.6),(3.7)and the exactitude of the Gauss-Lobatto quadratures,we obtain

Integrating by parts leads to

By using Lemma 3.1,we get

It remains to evaluate the first two terms on the right-hand side of(3.5).We deduce from(2.3)and(2.6)that,for all vN∈XN,

By using the definition(2.4),we get,for all vN∈XN,

Applying the integration by parts and(2.5)leads to

Moreover,it is known[20]that

If we take

then vN∈XN,and

Combining(3.2),(3.5),(3.8),(3.9),(3.10),and(3.11),we deduce

Using Lemma 3.1 and some rearrangement,we obtain

The desired estimate(3.3)is then obtained by dividing both sides by kvkX+k q kM.

3.2 Lower bound estimation

Now itis interesting to see whetherηis also a lower bound of ku-uNkX+k p-pNkM.Unfortunately,it is not true due to the poor inverse inequalities for polynomials.Nevertheless,inspired by the existing results on the a posteriori error estimation for hp-FEMofelliptic equations[19],we are able to derive a lower bound for the error ku-uNkX+k p-pNkMin term of the modified indicatorηθ,θ∈[0,1],to be defined below.

The estimation of the lower bound makes use of some polynomial inverse inequalities in weighted Sobolev spaces,which we recall below.

·Polynomial inverse estimates

Let Λ =(-1,1),we define the weight function ΦΛ=(1-x2).Then there holds the following inverse estimates.

Lemma 3.3.[3,19]Let-1＜α＜β,δ∈[0,1].Then the following estimates hold for all polynomials qN∈QN(Λ),

The generalization of the above inequalities to 2-dimensional case can be done by introducing the distance function as follows:

Then we have some similar inverse inequalities in 2D,as stated in the following lemma.

Lemma 3.4.[19]Letα,β∈R satisfy-1＜α＜βandδ∈[0,1].Then it holds for all polynomials qN∈QN(?Ω),

We will also need a known polynomial lifting result for the extension from an edge to the domain.

where cαis a constant depending only onα.

·Lower bound

Let Fi(i=s,d)be the mappings from?Ω to the element Ωi,FΓand FΓdbe the mappings from Λ to Γ and Γdrespectively.We define the following weight functions:

and the following error indicatorηθ,θ∈[0,1]:

The following lemmas are the main results about the lower bound estimates.

Lemma 3.6.Letθ∈[0,1],ζ＞0,then we have

Using-?·(-psI+2νD(us))=f,we obtain

Integrating by parts gives

(In the above estimation,we have used the factthatthe weighted function Φsis bounded).

Combining(3.14)and(3.15)completes the proof.

Proof.Note that?·us=0 in Ωs,then

Integrating by parts,and using Cauchy-Schwarz inequality,we deduce

Now we want to bound the H1-semi norm of v.

· Whenθ＞12,using the inverse estimate(3.13a)of Lemma 3.4 and the af fine transformation from the reference element Ω to Ωd,we can get:

This completes the proof.

Lemma 3.9.For the error indicatorηθ,4,we have

Proof.Note thatk-1ud+?pd=0 in Ωd,then

Lemma 3.10.For anyθ∈[0,1],ζ＞0,we have

Proof.Let

where cθdepends onθ.Using(2.5),we obtain

Integrating by parts yields

Plugging the first equation of(2.1)into the above equation gives

Then by using Cauchy-Schwarz inequality and Lemma 3.6 forθ=0,we obtain

whereζ＞0 is from Lemma 3.6.In order to apply Lemma 3.5,we distinguish two cases.

· Caseθ＞12.In this case Lemma 3.5 is directly applicable,which,together with the af fine equivalence and(3.18)withε=N-2,yields

Finally,combining(3.19)and(3.20)leads to the desired estimate.

Similarly,we will be able to derive some estimates forηθ,6andηθ,7.These estimates are stated in the following lemma without a detailed proof.

Lemma 3.11.For anyθ∈[0,1],ζ＞0,we have

By collecting the above results,we are now in a position to give the main result of this section,i.e.,the error lower bound of the numerical solution,in the following theorem.

Theorem 3.2.For anyθ∈[0,1],ζ＞0,it holds

Proof.It is a direct consequence of Lemma 3.6 to Lemma 3.11.

Remark 3.2.Compared to the well known classical inverse inequality(see,e.g.,[5]),the weighted inverse inequalities given in the lemma 3.3 have much weaker powers on the polynomial degree N.These better inequalities have allowed us to derive shaper lower bounds for the error indicatorsηθ,i.Obviously,the bigger isθthe shaper is the estimate.On the other side biggerθmeans heavier weight.In general there does not exist optimal choice ofθandζ.It depends on both the boundary distance function and solution regularity.

Remark 3.3.By using the inequalityη.Nθηθ,we can easily derive the lower bound of the error with the indicatorη.

4 Numerical tests

We now carry out some numerical tests to investigate the behavior of the numerical solution with respect to polynomial degree N.The main purpose is to verify the error indicators provided in the previous section.

The computational domain is Ω =(-1,3)×(-1,1)with Ωs=(-1,1)×(-1,1),Ωd=(1,3)×(-1,1).By using the nodal basis for the discrete spaces XNand MN,the spectral approximation(2.6)results in a discrete saddle point problem.This problem is then split,by applying the Uzawa algorithm,into two positive definite symmetric systems:one for the pressure and anotherfor the velocity.The pressure systemis solved by an inner/outer preconditioned conjugate gradient iteration.We consider the following exact analytical solution:

In Fig.1,we plot the error estimators as functions of polynomial degree N.The result presented in this figure indicates the exponential decay rate of the errors,as in this semilog representation the error curves are all straight lines.Note that the error indicatorη4is vanishing up to the machine precision for all polynomial degree N.Indeed,from the spectral approximation(2.6)it can be easily verified thatκ-1us,N+?pd,Nis identically zero for any N.

In Fig.2,we compare the lower error indicatorηand the sum of the velocity error and pressure error versus the polynomial degree N.We observe that the two error curves are very close to each other and nearly have the same slope.This is an implication that the lower bound of the velocity and pressure error given in Theorem 3.2 is optimal for this test problem.

5 Spectral element method

In this section,we extend the above estimates to the spectral element method.The domain Ω is split into a number of subdomains as follows:

where Ksand Kdstand for the element numbers in the Stokes and Darcy domain,respectively,KΓis the number of the sides that Γ contains.Denote the triangulation by

Figure 1:Error estimators versus the polynomial degree N.

Figure 2:Total error indicator and the error of the velocity and pressure solution as functions of N.

We intend to derive a posteriori error estimate for the spectral element method based on the triangulation T.We denote by?Ω the reference square,and suppose that each element T in T is the image of?Ω under an affine map FT:?Ω→T.

We denote by E(T)the set of the edges of element T,and let

In what follows,hT:=diam T,and hestands for the length of the edge e in E(T).

We define the piecewise polynomial spaces as follows:

Then let

The spectral element approximation to the problem(2.3)reads:find(uδ,pδ)∈ Xδ×Mδ,such that

where the bilinear forms aδ(·,·)and bδ(·,·)are defined respectively by:

We firstpresenta lemma which is simplification ofthe well-known resultsfor Cl′ementtype and Scott-Zhang type quasi interpolation operators[18]in the spectral element context.

where ΩT,Ωeare patches covering T and e with a few layers,respectively.

Next we derive an a posteriori error estimate for the solution of problem(5.1).To this end,we define for each T∈Tsthe a posteriori error indicator

where[φ]denotes the jump of the functionφacross the edge.Similarly,for each T∈Td,we set

Theorem5.1.Assume thatthe dataf∈Hr(Ωs),r＞1.Let(u,p)and(uδ,pδ)be the solutions of(2.3)and(5.1)respectively.Then we have

Proof.Similar to the single domain spectral method in Theorem 3.1,we denote eu=u-uδ,ep=p-pδ,then there existsv∈X,q∈M such that

We estimate the terms on the right sides.Setting

where Iδ-1q:C0(Ω)→PN-1,K(Ω),then?Iv∈Xδ,?Iq∈Mδ.

Now we turn to estimate the terms on the right-hand side of(5.3).

For the first two terms,we have,for allvδ∈Xδ,

By using integration by parts,Lemma 2.2,and the first and fourth equations in(2.1),we obtain

Now combining(5.4)～(5.7),using the Cauchy-Schwarz inequality and Lemma 5.1 gives

For the last term in(5.3)we have b(eu,q)=-b(uδ,q),and bδ(uδ,?Iq)=0.

Thus

Using(5.9)and the exactitude of the Gauss type quadratures,we obtain

Then using lemma 5.1,we obtain

Finally a direct combination of(5.3),(5.8),and(5.10)leads to the desired result(5.2).

Remark 5.1.Based on a similar technique as for the single domain spectral method,a lower bound for the spectral element case can be equally obtained.We omit the details of the proof to avoid a too technical discussion.

6 Concluding remarks

In this paper,we have considered a spectral(element)approximation to the Stokes/Darcy coupled equations.The purpose was to derive some a posteriori error estimates for the solutions of the discrete problems.The main results of this work include:1)We obtained the a posteriori error indicators for the single domain spectral method,and established their lower bounds and upper bounds;2)Some numerical tests are carried out to show sharpness of these estimates;3)A generalization to the spectral element case has been discussed,and a posteriori error estimate of the spectral element solution has been obtained.

Acknowledgments

The authors are grateful to the anonymous referees for their useful comments and suggestion.This work is partially supported by NSF of China(Grant numbers 11071203 and 91130002).

[1]I.Babuska and G.N.Gatica,A residual-based a posteriori error estimator for the Stokes-Darcy coupled problem,SIAM J.Numer.Anal.,48(2)(2010),498-523.

[2]C.Bernardi,Indicateurs d’erreur en h-n version des elements spectraux,RAIROModelisation Math et Analyse Numerique,30(1)(1996),1-38.

[3]C.Bernardi,N.Fi′etier,and R.G.Owens,An error indicator for mortar element solutions to the stokes problem.,IMA J.Numer.Anal.,21(4)(2001),857-886.

[4]C.Bernardi and Y.Maday,Spectral Methods,Handbook of numerical analysis,5(1997),209-485.

[5]C.Canuto,M.Y.Hussaini,A.Quarteroni,and T.A.Zang,Spectral methods:Fundamentals in single domains,Springer Verlag,2006.

[6]Y.Chen,N.Yi,and W.Liu,A Legendre-Galerkin spectral method for optimal control problems governed by elliptic equations,SIAM J.Numer.Anal,46(5)(2008),2254-2275.

[7]A.Chernov,Optimal convergence estimates for the trace of the polynomial-projection operator on a simplex,Math.Comp.,81(278)(2011),765-787.

[8]M.Cui and N.Yan,A posteriori error estimate for the Stokes-Darcy system,Math.Method Appl.Sci.,34(9)(2011),1050-1064.

[9]D.B.Das,V.Nassehi,and R.J.Wakeman,A finite volume model for the hydrodynamics of combined free and porous flow in sub-surface regions,Adv.Environ.Res.,7(1)(2002),35-58.

[10]M.Discacciati and A.Quarteroni,Navier-Stokes/Darcy coupling:modeling,analysis,and numerical Approximation,Rev.Mat.Complut,22(2)(2009),315-426.

[11]V.J.Ervin,E.W.Jenkins,and S.Sun,Coupled generalized non-linear Stokes flow with flow through a porous media,SIAM J.Numer.Anal,47(2)(2009),929-952.

[12]J.K.Guest and J.H.Pr′evost,Topology optimization of creeping fluid flows using a Darcy-Stokes finite element,Internat.J.Numer.Methods Engrg.,66(3)(2006),461-484.

[13]B.Guo.Recent Progress on A-posteriori Error Analysis for the p and hp Finite Element Methods,Contemporary Mathematics,383(2005),47-62.

[14]N.S.Hanspal,A.N.Waghode,V.Nassehi,and R.J.Wakeman,Development of a predictive mathematical model for coupled Stokes/Darcy flows in cross-flow membrane filtration,Chem.Eng.J.,149(1-3)(2009),132-142.

[15]P.Houston,C.Schwab,and E.S¨uli,Discontinuous hp-finite element methods for advectiondiffusion-reaction problems,SIAM J.Numer.Anal.,39(6)(2002),2133-2163.

[16]W.J¨ager and A.Mikeli′c,On the interface boundary condition of Beavers,Joseph,and Saffman,SIAM J.Appl.Math.,60(4)(2000),1111-1127.

[17]M.Le Bars and M.G.Worster,Interfacial conditions between a pure fluid and a porous medium:implications for binary alloy solidification,Int.J.Fluid Mech.Res.,550(2006),149-173.

[18]J.M.Melenk,hp-Interpolation of nonsmooth functions and an application to hp-A posteriori error estimation,SIAM J.Numer.Anal.,(2006),127-155.

[19]J.M.Melenk and B.I.Wohlmuth,On residual-based a posteriori error estimation in hp-FEM,Adv.Comput.Math.,15(1)(2001),311-331.

[20]A.Quarteroni,A.M.Quarteroni,and A.Valli,Numerical approximation of partial differential equations,volume 23.Springer Verlag,2008.

[21]P.G.Saffman,On the boundary condition at the surface of a porous medium,Stud.Appl.Math.,50(2)(1971),93-101.

[22]W.Wang and C.Xu,Spectral Methods Based on New Formulations for Coupled Stokes and Darcy Equations,J.Comput.Phys.,257(2014),126-142.

[23]C.Xu and Y.Maday.A global algorithm in the numerical resolution of the viscous/inviscid coupled problem,Chin.Ann.Math.Ser.B,18(1997),191-200.

?Corresponding author.Email addresses:cjxu@xmu.edu.cn(C.J.Xu)