Aposteriori error estimation of Subgrid multiscale stabilized finite element method for transient Stokes model
aa r X i v : . [ m a t h . A P ] J a n Aposteriori error estimation of Subgrid multiscalestabilized finite element method for transientStokes model
Manisha Chowdhury ∗ Indian Institute of Technology KanpurKanpur, Uttar Pradesh, India
Abstract
In this study, we present a novel stabilized finite element analysis fortransient Stokes model. The algebraic subgrid multiscale approach hasbeen employed to arrive at the stabilized coupled variational formulation.Derivation of the stabilized form as well as stability analysis of it’s fullydiscrete formulation are presented elaborately. Discrete inf - sup conditionfor pressure stabilization has been proven. For the time discretization thefully implicit schemes have been used. A detailed derivation of the apos-teriori error estimate for the stabilized subgrid multiscale finite elementscheme has been presented. Numerical experiment has been carried outto verify theoretically established order of convergence. For a long period of time different numerical methods such as finite difference, fi-nite volume and finite element methods have been employed to study the Stokesequations, widely used fluid flow mathematical model. This mixed problem doesnot satisfy inf - sup condition for equal-order velocity pressure interpolation andsometimes shows pressure instability. In this regard several stabilization tech-niques such as the streamline-upwind Petrov–Galerkin (SUPG) [5] -[7], the pres-sure stabilized Petrov–Galerkin (PSPG) [8]-[9], the Discontinuous Galerkin [10],the symmetric pressure stabilized Galerkin method [11], have been introduced todeal with this instability problem. In this paper we present Subgrid multicsalestabilized finite element analysis for transient Stokes model. Hughes in [3] hasintroduced the concept of stabilized multiscale subgrid method for Helmholtzequation and further developments are going on afterwards. Generally two ap-proaches of SGS stabilized formulation, namely algebraic approach, abbreviatedas
ASGS and orthogonal projection approach, known as
OSGS method, have ∗ Email addresses: [email protected](M.Chowdhury) STATEMENT OF THE PROBLEM been studied.
Badia and
Codina in [2] have studied both the approaches forunified Stokes-Darcy fluid flow problem and experimentally established equallywell performances of both the stabilized formulations. In [4] Codina presentsa study on comparison of stabilized finite element methods viz.
SU P G , GLS , SGS , T aylor − Galerkin etc. for solving diffusion-convection-reaction equationand experimentally shows that
SGS performs well in compared to other stabi-lized method. Here we have derived algebraic subgrid scale (
ASGS ) stabilizedfinite element method for the model using continuous velocities and pressurespaces across the inter-element boundaries.Therefore it is possible to eliminatethe jump terms in the variational formulation.In particular the
ASGS approach consists of algebraic approximation of thesubscales that arise from the decomposition of the exact solution field into re-solvable scale and unresolvable scale, have been used for finite element schemestabilization. Stabilization parameters are derived following the approach in [2]for
ASGS method. For time discretization fully implicit schemes have been cho-sen. In this study to ensure pressure stability discrete inf − sup condition hasbeen established for the choice of finite dimensional spaces and the stabilizedbilinear form is shown coercive. Furthermore aposteriori error estimate for thecurrent stabilized ASGS finite element method for the unsteady Stokes modelhas been derived. Generally error estimation provides important informationfor finding out the convergence rate of the numerical method and it’s depen-dency on the parameters present in the model. We have derived residual basedaposteriori error estimate which is computable and very beneficial for studyingadaptivity and control of solution error. This estimation provides convergencerates with respect to space and time both and the norm defined here to findthe estimation includes the standard norms involving both velocity and pressurevariable. Consequently it gives a wholesome information about the convergenceof the method. Numerical studies have shown the realization of theoretical or-der of convergence and the robustness of current stabilized ASGS finite elementmethod for Stokes system.Organisation of the paper is as follows: Section 2 introduces the model andit’s weak form. Next section presents space and time discretized variationalformulation. Section 4 describes the derivation of stabilized multiscale sub-grid formulation and stability analysis of the fully-discrete form. Next sectionhas elaborately described the derivation of aposteriori error estimation for thissubgrid formulation. At last section 6 contains numerical results to verify thenumerical performance of the method.
Let Ω ⊂ R d , d=2,3 be an open bounded domain with piecewise smooth bound-ary ∂ Ω. For the sake of simplicity in further calculations, we have consideredtwo dimensional model, but it can be easily extended for three dimensionalmodel. Let us now introduce here the system of fluid flow in Ω: Find u : Ω × STATEMENT OF THE PROBLEM [0,T] → R , p : Ω × [0,T] → R such that, ∂ u ∂t − µ ∆ u + ▽ p = f in Ω × [0 , T ] ▽ · u = 0 in Ω × [0 , T ] u = on ∂ Ω × [0 , T ] u = u at t = 0 (1)where ( u , p ) is a pair of Stokes velocities and pressure respectively and µ is theStokes dynamic viscosity and f is the force term.In operator form the above system of equations can be written as: find U =( u ,p) M ∂ t U + L ( U ) = F (2)where M, a matrix = diag(1,1,0), ∂ t U = [ ∂ u ∂t , ∂p∂t ] T L ( U ) = (cid:20) − µ ∆ u + ▽ p ▽ · u (cid:21) and F = (cid:20) f (cid:21) Let us introduce the adjoint L ∗ of L as follows, L ∗ ( U ) = (cid:20) − µ ∆ u − ▽ p − ▽ · u (cid:21) Now we impose suitable assumptions, that are necessary to conclude the resultsfurther, on the coefficients mentioned above. (i) µ is positive constant. (ii) The spaces of continuous solutions u = ( u , u ) and p are taken as thefollowing: (a) The velocities u , u ∈ L ∞ (0 , T ; H (Ω)) T C (0 , T ; H (Ω)) (b) and the pressures p ∈ L ∞ (0 , T ; H (Ω)) T C (0 , T ; L (Ω)). (iii) Additional assumptions imposed on continuous velocity solution are: u tt and u ttt are taken to be bounded functions on Ω. Weak formulation : To present variational formulation let us first introducethe spaces V u and Q for velocity and pressure respectively in the following: V u = { v ∈ ( H (Ω)) : v = on ∂ Ω } Q = { q ∈ L (Ω) : Z Ω q d Ω = 0 } (3)The variational formulation of (1) is to find U = ( u , p ) ∈ V u × Q (= V ) suchthat ∀ V = ( v , q ) ∈ V ( M ∂ t U , V ) + a S ( u , v ) − b ( v , p ) + b ( u , q ) = L ( V ) (4)where ( M ∂ t U , V ) = ( ∂ t u , v ) ; a S ( u , v ) = R Ω S µ ▽ u : ▽ v ; b ( u , q ) = R Ω ( ▽ · u ) q and L ( V ) = ( f , v ) S The above variational formulation can be written in the following compact way:find U ∈ V such that( M ∂ t U , V ) + B ( U , V ) = L ( V ) ∀ V ∈ V (5)3 DISCRETE FORMULATIONS let the domain Ω be discretized into finite numbers of subdomains Ω k fork=1,2,..., n el , where n el is the total number element subdomains. Let h k bethe diameter of each sub-domain Ω k .Let h = max k =1 , ,...n el h k and ˜Ω = S n el k =1 Ω k be the union of interior elements.Let V h u and Q h be finite dimensional subspaces of V u and Q respectively suchthat V h u = { v ∈ V u : v (Ω k ) = P l (Ω k ) } , and Q h = { q ∈ Q : q (Ω k ) = P l − (Ω k ) } where P l (Ω k ) denotes complete polynomial of order l respectively over each Ω k for k=1,2,..., n el .The discrete formulation is to find U h ∈ V h u × Q h (= V h ) such that( M ∂ t U h , V h ) + B ( U h , V h ) = L ( V h ) ∀ V h ∈ V h (6) For time discretization let us introduce the following uniform partition of thetime interval [0,T]: for time step size dt = TN , where N is a positive integer, n -thtime step t n = ndt and for given 0 ≤ θ ≤ f n := f ( · , t n ) f or ≤ n ≤ Nf n,θ := 12 (1 + θ ) f ( n +1) + 12 (1 − θ ) f n f or ≤ n ≤ N − dt := t n +1 − t n ∂ t f n,θ := f n,θ − f n θ dt (7)We see for θ = 0 the discretization follows Crank-Nicolson formula and for θ = 1it is backward Euler discretization rule.For sufficiently smooth function f ( t ), using the Taylor series expansion aboutt= t n,θ , we will have f n +1 = f ( t n,θ ) + (1 − θ ) dt ∂f∂t ( t n,θ ) + (1 − θ ) dt ∂ f∂t ( t n,θ ) + O ( dt ) f n = f ( t n,θ ) − (1 + θ ) dt ∂f∂t ( t n,θ ) + (1 + θ ) dt ∂ f∂t ( t n,θ ) + O ( dt ) (8)We have considered here t n,θ − t n = (1+ θ ) dt Multiplying the above first and second sub-equations in (14) by θ and − θ respectively and then adding them we will have the following f n,θ = f ( t n,θ ) + 18 (1 + θ )(1 − θ ) dt ∂ f∂t ( t n,θ ) + O ( dt ) (9)4 STABILIZED MULTISCALE FORMULATION
Let u n,θ , p n,θ be approximations of u ( x , t n,θ ) , p ( x , t n,θ ) respectively. Now byTaylor series expansion [15],we have u n +1 − u n dt = u t ( x , t n,θ ) + T E | t = t n,θ ∀ x ∈ Ω (10)where the truncation error
T E | t = t n,θ ≃ T E n,θ depends upon time-derivativesof the respective variables and dt [15]. Now applying assumption (iii) on u tt and u ttt we will have another property as follows: k T E n,θ k ≤ ( ˜ C dt if θ = 1˜ C dt if θ = 0 (11)According to this rule we need to solve for U n,θh , ∀ V h ∈ V h ( M ∂ t U n,θh , V h ) + B ( U n,θh , V h ) = L ( V h ) + ( TE n,θ , V h ) (12)and the exact solution U n,θ satisfies the above equation in the following way: ∀ V h ∈ V h ( M ∂ t U n,θ , V h ) + B ( U n,θ , V h ) = L ( V h ) + ( TE n,θ , V h ) (13) Now we start deriving stabilized formulation with decomposing additively theexact solution into the resolvable scale U h , which is a finite element solutionand unresolvable scale term U ′ , known as subgrid scale. U = U h + U ′ (14)where U ′ belongs to a space V ′ which completes U h in U . The test function V can too be decomposed likewise into the components V h and V ′ . Followingthe classical approach we substitute (13) in (4) as follows:( M ∂ t U h , V h ) + B ( U h , V h ) + ( M ∂ t U ′ , V h ) + B ( U ′ , V h ) = L ( V h ) ∀ V h ∈ V h ( M ∂ t U h , V ′ ) + B ( U h , V ′ ) + ( M ∂ t U ′ , V ′ ) + B ( U ′ , V ′ ) = L ( V ′ ) ∀ V ′ ∈ V ′ (15)Now integrating the second sub-equation of (14) and applying suitable boundaryconditions we have Z ˜Ω V ′ · [ M ∂ t U ′ + L U ′ ] = Z ˜Ω V ′ · [ F − M ∂ t U h + L U h ] (16)Consideration of continuous velocities, pressure at the inter-element boundariesmakes the jump term vanishes in the above equation and we obtain over eachelement Ω k M ∂ t U ′ + L U ′ = F − M ∂ t U h + L U h = T he Residual ( R ( U h )) (17)5 STABILIZED MULTISCALE FORMULATION along with boundary condition on U ′ which is not known. Now we solve for U ′ the above equation. Let us assume an approximation of the differential operatorin this way: L ≈ τ − k , where τ k is a matrix whose components are known asstabilization parameters [2]. Now substituting this in (16) and applying timediscretization rule we have1 θ dt M ( U ′ − U ′ n ) + τ − k U ′ = R ( U h ) (18)this implies an expression for subgrid scale term in the following U ′ = ( 1 θ dt M + τ − k ) − { R ( U h ) + 1 θ dt M U ′ n } = τ ′ k ( R ( U h ) + d ) (19)where the matrices τ ′ k = ( θ dt M + τ − k ) − and d = θ dt M U ′ n .Now substituting the result (16) in the first sub-problem of (14) followed byintegrating the fourth term once and substituting the expression (18) obtainedfor U ′ , we have( M ∂ t U h , V h ) + B ( U h , V h ) + n el X k =1 ( R ( U h ) − τ − k τ ′ k ( R ( U h ) + d ) , V h ) k + n el X k =1 ( τ ′ k ( R ( U h ) + d ) , L ∗ V h ) k = L ( V h ) ∀ V h ∈ V h (20)Now expanding the residual term we have the final form of stabilized formulationin the following( M ∂ t U h , V h ) + B ASGS ( U h , V h ) = L ASGS ( V h ) ∀ V h ∈ V h (21)where B ASGS ( U h , V h ) = B ( U h , V h )+ P n el k =1 ( τ ′ k ( M ∂ t U h + L U h − d ) , −L ∗ V h ) Ω k − P n el k =1 (( I − τ − k τ ′ k )( M ∂ t U h + L U h ) , V h ) Ω k − P n el k =1 ( τ − k τ ′ k d , V h ) Ω k L ASGS ( V h ) = L ( V h ) + P n el k =1 ( τ ′ k F , −L ∗ V h ) Ω k − P n el k =1 (( I − τ − k τ ′ k ) F , V h ) Ω k where the stabilization parameter τ k is in matrix form as τ k = diag ( τ k , τ k , τ k ) = (cid:20) τ k I × τ k (cid:21) and τ ′ k = ( 1 dt M + τ − k ) − = (cid:20) τ k dtdt + ρτ k I × τ k (cid:21) = diag ( τ ′ k , τ ′ k , τ ′ k )and the matrix d = P n +1 i =1 ( dt M τ ′ k ) i ( F − M ∂ t U h − L ( U h )) =[ d , d ] T It can be easily observed that d is always 0 due to the matrix M. Now theexpressions of stabilization parameters [2] are τ k = h c µ and τ k = c τ k , where c , c are constants.According to the time discretization rule we need to solve for U n,θh , ∀ V h ∈ V h ( M ∂ t U n,θh , V h ) + B ASGS ( U n,θh , V h ) = L ASGS ( V h ) + ( TE n,θ , V h ) (22)6 .1 Stability analysis of fully-discrete stabilized form4 STABILIZED MULTISCALE FORMULATION Let us first mention discrete inf - sup condition for discrete formulation. Lemma 1.
Discrete inf-sup condition : There exists a constant β > ,independent of h , such that inf q h ∈ Q h sup v h ∈ V dh | b ( q h , v h ) |k v h k k q h k ≥ β Proof.
Let π h be L projection on V h satisfying the following relations:(i)Stability estimates [13] for u ∈ H (Ω), k π h u k ≤ C k u k and k π h u k ≤ C k u k (ii) Interpolation estimate [14] for v ∈ ( H (Ω)) d k ▽ · ( v − π h v ) k ≤ ¯ Ch r k ▽ · v k r for 0 ≤ r ≤ . Let q h ∈ Q h . From [12] there exists v ∈ ( H (Ω)) d such that ▽ · v = q h and k v k ≤ C k q h k k q h k = ( q h , ▽ · v ) = ( q h , ▽ · v − ▽ · π h v ) + ( q h , ▽ · π h v ) ≤ k q h k k ▽ · ( v − v h ) k + | b ( q h , π h v ) |≤ k q h k k v − v h k + | b ( q h , π h v ) |≤ k q h k ¯ C k v k + | b ( q h , π h v ) |≤ ¯ CC k q h k + | b ( q h , π h v ) | (23)which implies | b ( q h , π h v ) | ≥ (1 − ¯ CC ) k q h k ≥ (1 − ¯ CC ) C k q h k k v k ≥ (1 − ¯ CC ) C C k q h k k π h v k (24)Choose β = (1 − ¯ CC ) C C provided ¯ CC < U nh ∈ V h find U n +1 h ∈ V h such that ∀ V h ∈ V h ( Mdt U n +1 h , V h )+ B ( U n +1 h , V h ) − n el X k =1 (( I − τ − k τ ′ k )( Mdt U n +1 h + L ( U n +1 h )) , V h ) Ω k + n el X k =1 ( τ ′ k ( Mdt U n +1 h + L ( U n +1 h )) , −L ∗ ( V h )) Ω k = L ASGS ( V h ) + ( Mdt U nh , V h ) − n el X k =1 (( I − τ − k τ ′ k )( Mdt U n +1 h , V h ) Ω k + n el X k =1 ( τ − k τ ′ k d , V h ) Ω k − n el X k =1 ( τ ′ k d , L ∗ ( V h )) Ω k + n el X k =1 ( τ ′ k ( Mdt U nh , −L ∗ ( V h )) Ω k (25)7 .1 Stability analysis of fully-discrete stabilized form4 STABILIZED MULTISCALE FORMULATION Now dropping the superscripts let us denote the unknown part in the following¯ B ASGS ( U h , V h ; dt ) := B ( U h , V h ) − n el X k =1 (( I − τ − k τ ′ k )( Mdt U h + L ( U h )) , V h ) Ω k + ( Mdt U h , V h ) + n el X k =1 ( τ ′ k ( Mdt U h + L ( U h )) , −L ∗ ( V h )) Ω k (26) Theorem 1.
For regular partitions satisfying inverse inequalities and assuminga condition dt > ¯ Ch there exists positive parameters C i s (for i=1,2,3) depend-ing upon h such that ¯ B ASGS ( U h , U h ; dt ) ≥ C k u h k + C k ▽ u h k + C k ▽ p h k Proof.
Expanding the terms of (20) we have B ASGS ( U h , U h ; dt ) = 1 dt k u h k + a S ( u h , u h ) − n el X k =1 ( τ k dt + τ k I × ( u h dt − µ ∆ u h + ▽ p h ) , u h ) Ω k + n el X k =1 ( τ k dtdt + τ k I × ( u h dt − µ ∆ u h + ▽ p h ) , µ ∆ u h + ▽ p h ) Ω k (27)let us look at the terms separately. a S ( u h , u h ) = µ k ▽ u h k Now applying Cauchy-Schwarz inequality on the next terms we have( τ k dt + τ k I × u h dt , u h ) Ω k ≤ dt ( P n el k =1 | τ k || dt + τ k | ) k u h k k and applying inverse inequalities in the following intermediate steps we have( τ k dtdt + τ k I × u h dt , µ ∆ u h + ▽ p h ) Ω k ≤ τ k dt + τ k ( u h , µ ∆ u h + ▽ p h ) Ω k ≤ τ k dt + τ k { µ k u h k k k ∆ u h k k + k u h k k k ▽ p h k k }≤ τ k dt + τ k { ǫ k u h k k + µ C I ǫ h k ▽ u h k k + 12 ǫ k ▽ p h k k } (28)Now the next term similarly as above( τ k dt + τ k I × ( − µ ∆ u h + ▽ p h ) , u h ) Ω k ≤ τ k dt + τ k { ǫ k u h k k + µ C I ǫ h k ▽ u h k k + 12 ǫ k ▽ p h k k } (29)and( τ k dtdt + ρτ k I × ( − µ ∆ u h + ▽ p h ) , µ ∆ u h + ▽ p h ) Ω k ≥ τ k dtdt + ρτ k {k▽ p h k k − µ k ∆ u h k k } ≥ τ k dtdt + ρτ k {k▽ p h k k − µ C I h k▽ u h k k } (30)8 ERROR ESTIMATION
Now combining all the results¯ B ASGS ( U h , U h ; dt ) ≥ C k u h k + C k ▽ u h k + C k ▽ p h k (31)where C = (1 − ǫ τ ) dt + τ = τ − − ǫ dtτ − +1 = c µh − ǫ dtτ − +1 = c µ − ǫ h h ( dtτ − +1) C = τ [ µ h ( c − C I ǫ ( dt + τ ) − C I dt ( dt + τ ) )] ≥ τ ¯ Cτ − = ¯ C C = τ dt + τ ( dt − ǫ ) = ( dt − ǫ − ) dtτ − +1 Choose arbitrary ǫ in such a way that dt > ǫ − > ¯ Ch holds and therefore C and C can be made positive. ¯ C is positive constants by choice of stabilizationparameters. We start this section with introducing the projection operator corresponding toeach unknown variable followed by notation of error and it’s component wisesplitting. Later we derive aposteriori error estimates.
Let us introduce the projection operator for each of these error components.(I)For any u ∈ H (Ω) × H (Ω) we assume that there exists an interpolation I h u : H (Ω) × H (Ω) −→ V h u satisfying b ( u − I h u u , q h ) = 0 ∀ q h ∈ Q h (II) Let I hp : H (Ω) −→ Q hs be the L orthogonal projection given by R Ω ( p − I hp p ) q h = 0 ∀ q h ∈ Q h and for any p ∈ H (Ω)Let e = ( e u , e p ) denote the error where the components are e u = ( e u , e u ) =( u − u h , u − u h ) and e p = ( p − p h ). Now each component of the error can besplit into two parts interpolation part, E I and auxiliary part, E A as follows: e u = ( u − u h ) = ( u − I hu u ) + ( I hu u − u h ) = E Iu + E Au Similarly, e u = E Iu + E Au and e p = E Ip + E Ap .At this point let us mention the standard interpolation estimation result [1]in the following: for any exact solution with regularity upto (m+1) k v − I hv v k l = k E Iv k l ≤ C ( p, Ω) h m +1 − l k v k m +1 (32)where l ( ≤ m + 1) is a positive integer and C is a constant depending on mand the domain. For l=0 and 1 it implies standard L (Ω) and H (Ω) normsrespectively. For simplicity we will use k · k instead of k · k to denote L (Ω)norm. Now we put some results using the properties of projection operators andthese results will be used in error estimations. Result 1. ( ∂∂t E I,n u , v h ) = 0 v h ∈ V h (33)9 .2 Aposteriori error estimation 5 ERROR ESTIMATION Result 2.
For any given auxiliary error E A,n and unknown E A,n +1 ( ∂∂t E A,n , E
A,n,θ ) ≥ dt ( k E A,n +1 k − k E A,n k ) (34) In this section we are going to derive residual based aposteriori error estimation.Before deriving error estimations let us define required norms the estimations.Let us consider the space ˜ V := L (0 , T ; V s ) T L ∞ (0 , T ; Q s ) and it’s associatednorm is denoted by ˜ V -norm. For the functions g , g , g belonging to the spaces L (0 , T ; L (Ω)), L (0 , T ; H (Ω)), ˜ V respectively norms over these spaces, ab-breviated as L ( L ), L ( H ), ˜ V are defined in the following k g k L ( L ) = N − X n =0 Z t n +1 t n Z Ω | g n,θ | dt k g k L ( H ) = N − X n =0 Z t n +1 t n ( Z Ω | g n,θ | + Z Ω | ∂g ∂x n,θ | + Z Ω | ∂g ∂y n,θ | ) dt k g k V = max ≤ n ≤ N k g n k + k g k L ( H ) (35) Theorem 2.
For computed velocity u h and pressure p h belonging to V h u and Q h satisfying (32)-(33), assume dt is sufficiently small and positive, and sufficientregularity of exact solution in equation (1). Then there exists a constant ¯ C ,independent of u , p and depending on the residual such that k u − u h k V + k p − p h k L ( L ) ≤ ¯ C ( R )( h + dt r ) (36) where R is the residual vector and r = ( , if θ = 12 , if θ = 0 (37) Proof.
We estimate aposteriori error by dividing the procedure into two parts.In the first part we find error bound corresponding to velocity and concentration followed by the second part estimating error associated with the pressure term.Let us first introduce the residual vector corresponding to each equations R = (cid:20) f − { ∂ u h ∂t − µ ∆ u h + ▽ p h }− ▽ · u h (cid:21) = (cid:20) R R (cid:21) First part : We have ∀ V ∈ Vµ | v | ≤ B ( V , V ) = a S ( v , v ) (38)Since e ∈ V we substitute the errors e u into the above relation and adding fewterms in both sides we have( ∂e u ∂t , e u )+ µ l k e u k ≤ ( ∂e u ∂t , e u )+ a S ( e u , e u )+ b ( e u , e p ) − b ( e u , e p )+ µ l k e u k (39)10 .2 Aposteriori error estimation 5 ERROR ESTIMATION Now first we will find a lower bound of
LHS and then upper bound for
RHS and finally combining them we will get aposteriori error estimate. Applying(33) on the first term of
LHS we have the following relations( e n +1 u − e n u dt , e n,θ u ) ≥ dt ( k e n +1 u k − k e n u k ) (40)Hence 12 dt ( k e n +1 u k − k e n u k ) + µ k e n,θ u k ≤ LHS ≤ RHS (41)Now our job is to find upper bound for
RHS and to reach at the desired es-timates let us divide it into two broad parts by splitting errors in each of theterms in the following way:
RHS = [( e n +1 u − e n u dt , E I,n,θ u ) + a S ( e n,θ u , E I,n,θ u ) + b ( e n,θ u , E I,n,θp ) − b ( E I,n,θ u , e n,θp )]+ [( e n +1 u − e n u dt , E A,n,θ u ) + a S ( e n,θ u , E A,n,θ u ) + b ( e n,θ u , E A,n,θp ) − b ( E A,n,θ u , e n,θp )] + µ k e n,θ u k = RHS I + RHS A + µ k e n,θ u k (42)Our aim is to bring residual into context and for this purpose RHS I involvinginterpolation error terms can be estimated as follows: we have ∀ V ∈ V ( e n +1 u − e n u dt , v )+ a S ( e n,θ u , v ) − b ( v , e n,θp )+( ▽· e n,θ u , q ) = Z Ω R n,θ · v + Z Ω R n,θ q (43)Now substituting v , q in the above expressions by E I,n,θ u , E I,n,θp respectively, wehave the
RHS I as, RHS I = Z Ω ( R n,θ · E I,n,θ u + R n,θ E I,n,θp ) ≤ h {k R n,θ k ( 1 + θ k u n +1 k + 1 − θ k u n k ) + C k R n,θ k ( 1 + θ k p n +1 k + 1 − θ k p n k ) }≤ h ( ¯ C k R n,θ k + ¯ C k R n,θ k ) (44)The parameters ¯ C i , for i=1,2,3,4, are coming from imposing assumption (iv) .Now we are going to estimate of RHS A . For that we employ subgrid formulation(20). Subtracting (20) from the variational finite element formulation satisfied11 .2 Aposteriori error estimation 5 ERROR ESTIMATION by the exact solution we have ∀ V h ∈ V h ( e n +1 u − e n u dt , v h ) + a S ( e n,θ u , v h ) − b ( v h , e n,θp ) + b ( e n,θ u , q h )= n el X k =1 { ( τ ′ k ( R n,θ + d ) , −L ∗ V h ) Ω k − (( I − τ − k τ k ) R n,θ , V h ) Ω k + ( τ − k τ k d , V h ) Ω k } + ( TE n,θ , v h ) (45)where the column vector d = [ d , d ] T = P n +1 i =1 ( dt M τ ′ k ) i ( F − M ∂ t U h − L U h ) = P n +1 i =1 ( dt M τ ′ k ) i R . Hence we have the components d = ( P n +1 i =1 ( dt τ ′ ) i ) I d × d R n,θ and d = 0.Now expanding the terms in (44) further and substituting V h by ( E A,n,θ u , E A,n,θp )in the above equation we have
RHS A as follows RHS A = n el X k =1 [( τ ′ I × { R n,θ + d } , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k + τ ′ ( R n,θ , ▽ · E A,n,θ u ) Ω k + ((1 − τ − τ ′ ) I × R n,θ , E A,n,θ u ) Ω k + ( τ − τ ′ I × d , E A,n,θ u ) Ω k ]+( TE n,θ , E A,n,θ u ) (46)Now we estimate each term separately. Before going to further calculationslet us mention an important observation: By the virtue of the choices of thefinite element spaces V h u and Q h , we can clearly say that Ω k every functionbelonging to that spaces and their first and second order derivatives all arebounded over each element sub-domain. We can always find positive finite realnumbers to bound each of the functions over element sub-domain. ApplyingCauchy-Schwarz inequality followed by this observation on the following termsof (45) we have n el X k =1 ( τ ′ I × R n,θ , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k ≤ | τ | TT − C τ ( n el X k =1 D B k ) k R n,θ k n el X k =1 ((1 − τ − τ ′ ) I × R n,θ , E A,n,θ u ) Ω k ≤ | τ | T − C τ ( X i =1 n el X k =1 B i k ) k R n,θ k (47)where the constants D B k and B i k appear due to imposing bounds on ∆ E A,n,θ u , ▽ E A,n,θp and E A,n,θ u over each sub-domain Ω k and C τ is upper bound on τ .Now the terms containing the components of d can be estimated in the following12 .2 Aposteriori error estimation 5 ERROR ESTIMATION way: n el X k =1 ( τ ′ I × d , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k = n el X k =1 ( τ ′ { n +1 X i =1 ( 1 dt τ ′ ) i } I × R n,θ , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k ≤ n el X k =1 ( τ ′ { ∞ X i =1 ( 1 dt τ ′ ) i } I × R n,θ , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k = n el X k =1 ( τ dt + τ I × R n,θ , µ ∆ E A,n,θ u + ▽ E A,n,θp ) Ω k ≤ | τ | C τ T − C τ ( n el X k =1 D B k ) k R n,θ k n el X k =1 ( τ − τ ′ I × d , E A,n,θ u ) Ω k ≤ | τ | T − C τ ( X i =1 n el X k =1 B i k ) k R n,θ k (48)Now we estimate the remaining terms as follows: n el X k =1 τ ′ ( R n,θ , ▽ · E A,n,θ u ) Ω k = n el X k =1 τ ′ ( R n,θ , ▽ · e n,θ u ) Ω k − n el X k =1 τ ′ ( R n,θ , ▽ · E I,n,θ u ) Ω k = n el X k =1 τ ′ ( ▽ · e n,θ u , ▽ · e n,θ u ) Ω k − n el X k =1 τ ′ ( ▽ · e n,θ u , ▽ · E I,n,θ u ) Ω k (49)Applying Cauchy-Schwarz inequality followed by Young’s inequality in the fol-lowing steps we have: ≤ C τ ( X i =1 k ∂e n,θui ∂x i k ) + C τ ( X i =1 k ∂e n,θui ∂x i k )( X i =1 k ∂E I,n,θui ∂x i k ) ≤ C τ X i =1 k ∂e n,θui ∂x i k + ǫ ′ C τ X i =1 k ∂e n,θui ∂x i k + C τ ǫ ′ X i =1 k ∂E I,n,θui ∂x i k ≤ C τ [(2 + ǫ ′ ) X i =1 k e n,θui k + h ǫ ′ X i =1 ( 1 + θ k u n +1 i k + 1 − θ k u ni k ) ] ≤ (2 + ǫ ′ ) C τ k e n,θ u k + h C τ ǫ ′ ¯ C (50)where the parameter ¯ C comes for applying assumption (iv) . Now the estima-13 .2 Aposteriori error estimation 5 ERROR ESTIMATION tion of terms involving trancation is in the following:( TE n,θ , E A,n,θ u ) = ( TE n,θ , e n,θ u ) − ( TE n,θ , E I,n,θ u ) ≤ ǫ ′ k TE n,θ k + ǫ ′ k e n,θ u k + k E I,n,θ u k ) ≤ ǫ ′ k TE n,θ k + ǫ ′ {k e n,θ u k + h ( 1 + θ k u n +1 k + 1 − θ k u n k ) }≤ ǫ ′ k TE n,θ k + ǫ ′ k e n,θ u k + h ǫ ′ C (51)and this completes estimating all the terms of RHS A . Now applying P oincare inequality on last term in
RHS in (41) we have µ k e n,θ u k ≤ µC P | e n,θ u | ≤ µC P k e n,θ u k (52)Now this completes finding bounds for each term in the RHS of (41). Puttingcommon terms all together in the left hand side and multiplying them by 2 andthen integrating both sides over ( t n , t n +1 ) for n = 0 , ..., ( N −
1) , we have N − X n =0 ( k e n +1 u k − k e n u k ) + { µ l − ǫ ′ ) C τ − µC P − ǫ ′ } N − X n =0 Z t n +1 t n k e n,θ u k dt ≤ h N − X n =0 Z t n +1 t n { ¯ C k R n,θ k + ¯ C k R n,θ k + 2 C τ ǫ ′ ¯ C + h ǫ ′ ¯ C } dt +2 | τ | T − C τ [( T + C τ )( n el X k =1 D B k ) + 2 X i =1 n el X k =1 B i k ] N − X n =0 Z t n +1 t n k R n,θ k dt + 2 ǫ ′ N − X n =0 Z t n +1 t n k TE n,θ k dt (53)Choose the arbitrary parameters including C τ and the P oincare constant C P in such a way that all the coefficients in the left hand side can be made positive.Then taking minimum over the coefficients in the left hand side let us divideboth sides by them. Using property (10) associated with both implicit timediscretisation scheme and the fact that τ is of order h , we have arrived at thefollowing relation: k e u k V ≤ C ′ ( R )( h + dt r ) (54)where r = ( , if θ = 1 f or backward Euler rule , if θ = 0 f or Crank − N icolson scheme (55)This only completes one part of aposteriori estimation and in the next part wecombine the corresponding pressure part.14
NUMERICAL EXPERIMENT
Second part : By Galerkin orthogonality followed by property (I) of projectionoperator we have b ( v h , I h p − p h ) = ( ∂e u ∂t , v h ) + a S ( e u , v h ) (56)Integrating both sides with respect time and later applying Cauchy − Schwarz ’sinequality,
Y oung ’s inequality and the above result (53) we have N − X n =0 Z t n +1 t n b ( v h , E A,n,θp ) dt = N − X n =0 Z t n +1 t n { ( e n +1 u − e n u dt , v h ) + a S ( e n,θ u , v h ) } dt ≤ ¯ C ′ ( R )( h + dt r ) k v h k (57)Now applying this result on the following relation k I h p − p h k L ( L ) = k E Ap k L ( L ) = N − X n =0 Z t n +1 t n k E A,n,θp k dt ≤ N − X n =0 Z t n +1 t n sup v h b ( v h , E A,n,θp ) k v h k dt ≤ ¯ C ′ ( R )( h + dt r ) (58)Now combining the results obtained in the first and second part and applyinginterpolation estimate on pressure interpolation term E Ip , we finally arrive atthe following aposteriori error estimate: k u − u h k V + k p − p h k L ( L ) ≤ ¯ C ( R )( h + dt r ) (59) Remark 1.
These estimations clearly imply that the scheme is f irst orderconvergent in space with respect to total norm, whereas in time it is f irst or-der convergent for backward Euler time discretization scheme and second orderconvergent for Crank-Nicolson method.
In this section we have numerically verified the convergence rate establishedtheoretically under stabilized method in the previous section. For simplicity wehave considered bounded square domain Ω= (0,1) × (0,1). We have taken con-tinuous piecewise linear finite element(P1) space into account for approximatingboth the velocity and pressure variables and applied backward Euler time dis-cretization rule. Let us mention here the exact solutions : u = ( e − t x ( x − CONCLUSION
Time Grid ASGS methodstep size Total error RoC0.1 10 ×
10 0.06516110.05 20 ×
20 0.0349207 0.8999280.025 40 ×
40 0.0182159 0.938830.0125 80 ×
80 0.00931727 0.967220.00625 160 ×
160 0.0047026 0.986447Table 1: Total error and Rate of convergence(RoC) under
ASGS method fortransient Stokes model at T = 11) y ( y − y − , − e − t y ( y − x ( x − x − p = e − t (2 x − y − µ = 0 . τ = h µ and τ = τ Remark 2.
Table 1 is showing total error under
ASGS method at each meshsize and at different time steps.Clearly the order of convergence under
ASGS method is 1, which justifies theoretically established result.
The paper presents
ASGS stabilized finite element analysis of transient Stokesfluid flow model. Whereas this paper in one hand has elaborately derived thestabilized formulation, on other hand it has presented stability analysis of thestabilized formulation where the stabilized bilinear form is shown coercive anddiscrete inf - sup condition has been established to ensure pressure stability. Aswell as aposteriori error estimation has been carried out elaborately. It is es-sential to mention that the norm employed for error estimation consists of thefull norms corresponding to each variable belonging to their respective spaces.Therefore it provides a wholesome information about convergence of the method.Theoretically the rate of convergence for aposteriori error estimation turns outto be O ( h ) in space for different time discretization rules. In numerical sectionthe theoretically established result is well verified through a test case problem.This piece of work will definitely be useful for studying adaptivity and controlof solution error. 16 EFERENCES REFERENCES
References [1] Y. Amanbek, M.F. Wheeler, A priori error analysis for transient problemsusing Enhanced Velocity approach in the discrete-time setting, Journal ofComputational and Applied Mathematics 361, 459-471(2019).[2] S. Badia, R. Codina, Unified stabilized finite element formulations for theStokes and Darcy problems, SIAM J. Numer. Anal., 47(3), 1971–2000(2009).[3] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet toNeumann formulation, subgrid scale models, bubbles and the origins of sta-bilized methods, Comput. Methods Appl. Mech. Engrg. 127,387-401(1995).[4] R. Codina, Comparison of some finite element methods for solving thediffusion-convection-reaction equation, Comput. Methods Appl. Mech. En-grg. 156, 185-210 (1998).[5] J. Douglas Jr., J.P. Wang, An absolutely stabilized finite element methodfor the Stokes problem, Math. Comp. 52 (186), 495–508(1989)[6] L. Tobiska, R. Verfurth, Analysis of a streamline diffusion finite elementmethod for the Stokes and Navier–Stokes equations, SIAM J. Numer. Anal.33 (1), 107–127(1996).[7] C. Johnson, J. Saranen, Streamline diffusion methods for the incompressibleEuler and Navier–Stokes equations, Math. Comp. 47 (175) , 1–18(1986).[8] T.J.R. Hughes, L.P. Franca, M. Balestra, A new finite element formulationfor computational fluid dynamics. V. Circumventing the Babuska–Brezzi con-dition: a stable Petrov–Galerkin formulation of the Stokes problem accom-modating equal-order interpolations, Comput. Methods Appl. Mech. Engrg.59 (1), 85–99. (1986).[9] E. Burman, M. A. Fernández, Analysis of the PSPG method for the transientStokes’ problem,Comput. Methods Appl. Mech. Engrg. 200, 2882-2890 (2011).[10] E. Burman, M. A. Fernández, Galerkin finite element methods with sym-metric pressure stabilization for the transient Stokes equations: stability andconvergence analysis, SIAM J. Numer. Anal. 47, 409–439 (2008) .[11] N. Ahmed, S. Becher, G. Matthies, Higher-order discontinuous Galerkintime stepping and local projection stabilization techniques for the transientStokes problem, Comput. Methods Appl. Mech. Engrg. 313, 28-52 (2017).[12] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equa-tions: Theory and Algorithms, vol. 5 of Springer series in computationalmathematics. Springer, Berlin (1986).17
EFERENCES REFERENCES [13] E. Burman, M. A. Fernández, Continuous interior penalty finite elementmethod for the time-dependent Navier–Stokes equations: space discretizationand convergence, Numer. Math., 107,39–77(2007).[14] J. A. Wheeler, M. F. Wheeler, I. Yotovc, Enhanced velocity mixed finiteelement methods for flow in multiblock domains, Computational Geosciences6: 315–332(2002).[15] B. Rivi` ee