A superconvergent HDG method for Stokes flow with strongly enforced symmetry of the stress tensor
Matteo Giacomini Alexandros Karkoulias Ruben Sevilla Antonio Huerta

=0.5em

























































































































































































keywordsKeywords:
theoremTheorem
lemmaLemma
definition
definitionDefinition
remarkRemark




figs/









A superconvergent HDG method for Stokes flow with strongly enforced symmetry of the stress tensor
    
			  Matteo Giacomini

             Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain

             Alexandros Karkoulias

			  Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain

             Centre Internacional de Metodes Numerics en Enyinyeria (CIMNE), 

             Campus Nord UPC, Barcelona, Spain

			  Ruben Sevilla

             Zienkiewicz Centre for Computational Engineering, 

             College of Engineering, Swansea University, Wales, UK 

             Antonio Huerta

             Laboratori de Càlcul Numèric (LaCàN), 

             ETS de Ingenieros de Caminos, Canales y Puertos, 

             Universitat Politècnica de Catalunya, Barcelona, Spain


    December 30, 2023
============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================





This work proposes a superconvergent hybridizable discontinuous Galerkin (HDG) method for the approximation of the Cauchy formulation of the Stokes equation using same degree of polynomials for the primal and mixed variables.

The novel formulation relies on the well-known Voigt notation to strongly enforce the symmetry of the stress tensor.

The proposed strategy introduces several advantages with respect to the existing HDG formulations.

First, it remedies the suboptimal behavior experienced by the classical HDG method for formulations involving the symmetric part of the gradient of the primal variable. 

The optimal convergence of the mixed variable is retrieved and an element-by-element post-process procedure leads to a superconvergent velocity field, even for low-order approximations.

Second, no additional enrichment of the discrete spaces is required and a gain in computational efficiency follows from reducing the quantity of stored information and the size of the local problems.

Eventually, the novel formulation naturally imposes physical tractions on the Neumann boundary. 

Numerical validation of the optimality of the method and its superconvergent properties is performed in 2D and 3D using meshes of different element types.




Hybridizable discontinuous Galerkin,
Stokes flow,
Cauchy stress formulation,
Voigt notation,
Superconvergence






§ INTRODUCTION AND MOTIVATIONS




The interest in discontinuous Galerkin (DG) methods <cit.> has increased in the past years owing to their ability to construct high-order discretizations on unstructured meshes and to their flexibility in performing p-adaptivity.

Among the different techniques proposed in the literature to approximate incompressible flow problems, the recent growing interest towards hybridizable discontinuous Galerkin (HDG) methods <cit.> is due to multiple advantages these formulations have with respect to classical DG ones.

Concerning Stokes flow, several HDG formulations have been proposed in the literature <cit.> and the interested reader is referred to <cit.> for an overview on the topic.

The use of hybridization was first introduced with the local discontinuous Galerkin (LDG) method to circumvent the construction of divergence-free approximations of the velocity field <cit.>, see also <cit.>.

Moreover, owing to hybridization <cit.>, the globally coupled unknowns are defined on the boundary of the mesh elements and are connected solely to neighboring elements. Thus, the size of the global problem is greatly reduced.

In addition, it is worth noting that HDG allows equal interpolation for velocity, pressure and strain rate tensor, owing to an appropriate definition of the numerical flux and to the introduction of a stabilization parameter.

Thus, the limitations of using equal-order approximations for velocity and pressure in the incompressible limit, through the fulfillment of the Ladyzhenskaya-Babuška-Brezzi (LBB) condition, are circumvented by HDG.

In particular, Cockburn and co-workers <cit.> proved solvability and stability under the aforementioned assumptions, without the need of an enriched space for the mixed variable, or a reduced space for the hybrid one.

In <cit.>, optimal convergence rates of order k+1 are obtained for all the variables using equal degree of approximation k, whereas classical DG display suboptimal convergence of order k for the pressure and the gradient of the velocity.

A key aspect of HDG is the ability to construct a post-processed velocity field superconverging with order k+2 <cit.>. 

This is crucial when the superconvergent solution is sought to devise automatic procedures to perform p-adaptivity (cf. e.g. <cit.>).

Nevertheless, the classical HDG equal-order approximation for the Cauchy formulation is known to experience suboptimal convergence of the mixed variable and a loss of superconvergence of the post-processed velocity field using low-order approximations <cit.>.

Recently, in a series of publications <cit.>, Cockburn and co-workers devoted a great effort to develop a general framework, namely the M-decomposition, to devise superconvergent HDG discretizations.

This approach relies on enriching the local spaces for the approximation of the mixed variable by adding extra basis functions.

The number of these additional basis functions is not significantly big and in most cases it depends on the type of element under analysis and not on the degree of approximation k.

Despite only the size of the local problems increases and the additional computational effort is limited, it induces a more complex implementation compared to standard HDG methods.

Alternative HDG formulations achieve convergence of order k+2 for the velocity field when polynomials of degree k are chosen to approximate the hybrid variable <cit.>. These methods rely on utilizing smaller spaces for the mixed variable and larger ones for the velocity and exploiting a special stabilization function, the so-called reduced stabilization, to handle them.

Closely related approaches, namely the hybrid high-order (HHO) <cit.> and the hybridized weak Galerkin (HWG) <cit.> methods can also achieve the same orders of convergence. 

The present work devises a superconvergent hybridizable discontinuous Galerkin method for the Cauchy formulation of the Stokes equation using the same degree of approximation for the primal and mixed variables.

The rest of this paper is organized as follows.

In Section <ref>, first, the equations governing the Stokes flow are recalled. 

Then, according to the rationale introduced in <cit.> for the linear elasticity equation, the symmetry of the stress tensor is strongly enforced by means of a technique well known in the computational mechanics community, namely the Voigt notation for symmetric tensors.

The corresponding Cauchy formulation of the Stokes equation with strongly enforced symmetry of the stress tensor is derived.

In Section <ref>, an HDG discretization is introduced.

A local post-process procedure providing a superconvergent velocity field even for low-order approximations is discussed without resorting to the complex framework of the M-decomposition.

Moreover, contrary to other HDG formulations, the proposed method features a reduced number of degrees of freedom for the mixed variable and is computationally more efficient since the resulting local problems are smaller.

The novel HDG formulation is validated in Section <ref>. Extensive analysis of the optimal convergence and superconvergence rates of the primal, mixed and post-processed variables, for two and three dimensional problems is provided by means of numerical simulations. Special emphasis is placed on the influence of the stabilization parameter and on the robustness of the method using meshes of different element types.

Eventually, the capability of the method to accurately compute quantities of interest depending on the solution of the Stokes equation (e.g. the drag force) is discussed and Section <ref> summarizes the results of this paper.





§ STOKES FLOW WITH STRONGLY ENFORCED SYMMETRY OF THE STRESS TENSOR




In this section, the framework to handle symmetric tensors discussed by Fish and Belytschko in <cit.> is recalled and the governing equations of a Stokes flow are formulated using Voigt notation.

First, the Cauchy formulation of the Stokes equation is recalled.





 §.§ Cauchy formulation of the Stokes flow




Consider an open bounded domain Ω⊂ℝ^ with boundary ∂Ω = Γ_D ∪Γ_N, Γ_D ∩Γ_N = ∅ and  being the number of spatial dimensions. 

The strong form of the problem under analysis reads as follows:


    {
    - ·   = s      in Ω,
    ·   = 0        in Ω,
       = -p  + 2 ν      in Ω,
       = _D        on Γ_D,
    ·   = t      on Γ_N,
    .


where the pair (,p) represents the velocity and pressure fields and  is the Cauchy stress tensor.

The terms s, _D and t respectively are the volumetric source term, the Dirichlet boundary datum to impose the value of the velocity on Γ_D and the traction applied on the Neumann boundary Γ_N.

The third equation, known as Stokes law, provides the relationship between the stress tensor and the velocity and pressure variables, through the viscosity coefficient ν > 0, the × identity matrix  and the strain rate tensor , := 1/2(  + ^T ) being the symmetric part of the gradient.

It is well-known that the Cauchy and the velocity-pressure formulations of the Stokes equation are equivalent from the variational point of view.

Nevertheless, a major difference arises when considering the imposition of Neumann boundary conditions. 

On the one hand, natural boundary conditions for the Cauchy formulation enforce the value of the normal stress which represents a physical traction. 

On the other hand, the velocity-pressure formulation only accounts for the gradient of the velocity field instead of its symmetric part, leading to the imposition of the so-called pseudo-tractions.

Hence, the physical interpretation is lost <cit.>.

Within this context, an artificial handling of Neumann boundary conditions is required to impose physically meaningful tractions. 

This represents a drawback when dealing with real-life and industrial applications in which the enforcement of physically relevant quantities is a major constraint to perform reliable numerical simulations and compare them with experimental data.





 §.§ Voigt notation for symmetric tensors




The so-called Voigt notation relies on the idea of storing a second-order tensor in a vectorial format by appropriately rearranging its diagonal and off-diagonal components.

Consequently, the application of differential operators (e.g. symmetric gradient, divergence and ) and the geometrical projections (e.g. in the normal and tagential directions to a surface) may be expressed as matrix equations.

For this purpose, the rationale for the construction of differential operator and geometrical quantities using Voigt notation is recalled.

Consider the previously defined strain rate tensor .

Owing to its symmetry, only = (+1)/2 components (i.e. three in 2D and six in 3D) have to be stored and the following column vector in ^ is obtained:


    :=[e_11 ,  e_22 ,  e_12]^T
       in 2D,
    [e_11 ,  e_22 ,  e_33 ,  e_12 ,  e_13 ,  e_23]^T
       in 3D.


The components of the strain rate in Equation (<ref>) read as


    e_ij := ∂ u_i/∂ x_j + (1-δ_ij) ∂ u_j/∂ x_i,   for  i,j = 1,…, ,


where δ_ij is the classical Kronecker delta.

In order to retrieve the aforementioned strain rate tensor , the off-diagonal terms e_ij,  i ≠ j have to be multiplied by a factor 1/2, namely


    :=[   e_11 e_12/2; e_12/2   e_22 ]   in 2D,
    [   e_11 e_12/2 e_13/2; e_12/2   e_22 e_23/2; e_13/2 e_23/2   e_33 ]   in 3D.


Similarly, the symmetry of the stress tensor  is exploited to store only  components in the column vector


    :=[σ_11 , σ_22 , σ_12]^T
       in 2D,
    [σ_11 , σ_22 , σ_33 , σ_12 , σ_13 , σ_23]^T
       in 3D.






  §.§.§ Differential operators using Voigt notation




Following <cit.>, the strain rate tensor can be written as = by introducing the × matrix 


    :=[ ∂/∂ x_1       0 ∂/∂ x_2;       0 ∂/∂ x_2 ∂/∂ x_1 ]^T
       in 2D,
    [ ∂/∂ x_1       0       0 ∂/∂ x_2 ∂/∂ x_3       0;       0 ∂/∂ x_2       0 ∂/∂ x_1       0 ∂/∂ x_3;       0       0 ∂/∂ x_3       0 ∂/∂ x_1 ∂/∂ x_2 ]^T
       in 3D.


As previously done for the strain rate tensor by introducing the matrix  accounting for the symmetric part of the gradient, the vorticity vector is handled through its skew-symmetric part.

More precisely, the vorticity := × may be expressed in terms of Voigt notation as = through the × matrix , with =(-1)/2 being the number of rigid body rotations in the space (i.e. one in 2D and three in 3D):


    :=[-∂/∂ x_2 , ∂/∂ x_1 ]
       in 2D,
    [        0 -∂/∂ x_3  ∂/∂ x_2;  ∂/∂ x_3        0 -∂/∂ x_1; -∂/∂ x_2  ∂/∂ x_1        0 ]   in 3D.



The  of a vector  in two dimensions is a scalar quantity.

Nevertheless, it can also be computed by embedding  in the three dimensional space ^3 and setting its third component equal to zero. 

Within this contect, × may be interpreted as a vector whose magnitude is given by the aforementioned matrix operation  and pointing along the third direction.






 §.§ Cauchy formulation of the Stokes flow using Voigt notation




Owing to the notation introduced in this section, the Stokes constitutive law may be expressed as = - p +, where the vector ∈^ and the matrix ∈^× read as


    :=[1 ,  1 ,  0 ]^T
       in 2D,
    [1 ,  1 ,  1 ,  0 ,  0 ,  0 ]^T
       in 3D.     :=[      2ν   0_× 1; 0_× 1^T       ν ]   in 2D,
    [ 2ν 0_; 0_  ν ]   in 3D.


Moreover, the Neumann boundary condition applied on Γ_N can be written as ^T  = by introducing the × matrix


    :=[ n_1   0 n_2;   0 n_2 n_1 ]^T
       in 2D,
    [ n_1   0   0 n_2 n_3   0;   0 n_2   0 n_1   0 n_3;   0   0 n_3   0 n_1 n_2 ]^T
       in 3D.


accounting for the normal direction to the boundary.

Similarly, the projection of a vector along the tangential direction , namely a tangent line in 2D and a tangent surface in 3D, reads as · = ^T, being ∈^× the matrix


    :=[n_2 ,  -n_1 ]^T
       in 2D,
    [    0 -n_3  n_2;  n_3    0 -n_1; -n_2  n_1    0 ]   in 3D.


In order to rewrite Equation (<ref>) using Voigt notation, the divergence of a symmetric tensor is expressed in terms of the transpose of the matrix  accounting for the symmetric part of the gradient <cit.>.

In a similar fashion, recall that · = () and observe that the trace operator may be expressed via the vector  introduced in Equation (<ref>).

Combining the matrix forms of the symmetric gradient, the Stokes law and the normal direction presented above, the following formulation of the Stokes equation using Voigt notation is obtained:


    {
    -^T    = s      in Ω,
    
    ^T    = 0              in Ω,
       = - p +       in Ω,
       = _D        on Γ_D,
    ^T    =       on Γ_N.
    .






 §.§ Fundamental theorems using Voigt notation




In <cit.>, a generalized version of the Gauss's and Stokes' theorems using Voigt notation has been introduced. 

In order to construct the variational formulation of the problem under analysis, the following two lemmas are recalled.

 
	Consider a vector ∈^ and a symmetric × tensor  whose counterpart in Voigt notation is . It holds:
	
	
    ∫_∂Ω( ^T ) · dΓ = ∫_Ω·( ) dΩ + ∫_Ω( ^T ) · dΩ .

	


 
	Consider a vector ∈^. It holds:
	
	
    ∫_Ω dΩ = ∫_∂Ω^T  dΓ .

	


The proofs follow straightforwardly by rewriting (<ref>)-(<ref>) using the corresponding continuous differential operators, see <cit.>.





§ A HYBRIDIZABLE DISCONTINUOUS GALERKIN METHOD




HDG is a discontinuous Galerkin method with hybridization based on a mixed formulation.

First, it is defined the so-called broken computational domain by introducing a partition of the domain Ω in  disjoint subdomains Ω_e with boundaries ∂Ω_e.

The internal interface Γ reads as


    Γ := [ ⋃_e=1^∂Ω_e ] ∖∂Ω ,


whereas the mesh skeleton is given by the union of internal and Neumann boundary faces, namely Γ∪Γ_N.

In what follows, the classical  internal products for vector-valued functions in Ω_e ⊂Ω and ∂Ω_e ⊂Γ∪∂Ω are considered:


    (,)_Ω_e := ∫_Ω_e· dΩ  ,     ⟨, ⟩_∂Ω_e := ∑_Γ_i ⊂∂Ω_e∫_Γ_i· dΓ .


Moreover, owing to the piecewise discontinuous nature of the functions involved in the HDG formulation, the jump operator · is defined along each portion of the interface as the sum of the values from the element on the right and the left, Ω_e and Ω_l <cit.>:


    ⊙ = ⊙_e + ⊙_l .


The second-order problem in Equation (<ref>) may thus be written as a system of first-order equations as follows:


    { +    = 0      in Ω_e, and for e=1,… ,,
    
    ^T (  +   p )    = s      in Ω_e, and for e=1,… ,,
    
    ^T    = 0               in Ω_e, and for e=1,… ,,
       = _D           on Γ_D,
    
    ^T ( +   p)    = -      on Γ_N,
    ⊗   = 0      on Γ,
    
    ^T ( +  p)   = 0      on Γ,
    .


where  is the so-called mixed variable and the last two equations are the transmission conditions enforcing the continuity of respectively the velocity and the flux across the interface Γ. 


In the case of purely Dirichlet boundary conditions (i.e. Γ_N = ∅), an additional constraint is required to avoid the indeterminacy of the pressure. 

A common choice relies on imposing zero mean value of the pressure on the boundary (cf. e.g. <cit.>):


    1/|∂Ω|⟨ p, 1⟩_∂Ω = 0 .







 §.§ Strong form of the local and global problems




In a series of papers by Cockburn and co-workers <cit.>, the hybridizable discontinuous Galerkin formulation for Stokes flow has been theoretically and numerically analyzed.

Starting from the mixed formulation on the broken computational domain in Equation (<ref>), HDG features two stages.

First, a set of  local problems are defined element-by-element to compute (_e,_e,p_e) for e=1,… ,:


    {_e + _e    = 0      in Ω_e
    
    ^T _e + ^T   p_e    = s      in Ω_e
    
    ^T _e    = 0                 in Ω_e
    _e    = _D           on ∂Ω_e ∩Γ_D,
    _e    =       on ∂Ω_e ∖Γ_D,
    .


where  is an independent variable representing the trace of the velocity on the mesh skeleton Γ∪Γ_N.

Remark that Equation (<ref>) is a purely Dirichlet boundary value problem. As previously observed, an additional constraint has to be added to remove the indeterminacy of the pressure, namely


    1/|∂Ω_e|⟨ p_e, 1 ⟩_∂Ω_e = ρ_e ,


where ρ_e denotes the mean pressure on the boundary of the element Ω_e.

Hence, for e=1,… , the local problem in Equation (<ref>) provides (_e,_e,p_e) in terms of the global unknowns  and ρ.

The trace of the velocity  and the mean pressure ρ on the element boundaries are determined by solving the global problem accounting for the transmission conditions and the Neumann boundary condition:


    {⊗   = 0      on Γ,
    
    ^T ( +  p)   = 0      on Γ,
    
    ^T ( +  p)    = -      on Γ_N.
    .


The first equation is automatically satisfied due to the Dirichlet boundary condition _e = imposed in the local problems and the unique definition of the hybrid variable  on each face of the mesh skeleton.

Moreover, the divergence-free condition in the local problem induces the following compatibility condition for each element Ω_e,  e=1,… ,


    ⟨·_e , 1 ⟩_∂Ω_e ∖Γ_D + ⟨_D ·_e, 1 ⟩_∂Ω_e∩Γ_D = 0 .


Consider the Voigt counterpart ^T _e = 0 of the aforementioned constraint (cf. Equation (<ref>)).

The resulting compatibility condition reads as


    ⟨^T _e  , 1 ⟩_∂Ω_e ∖Γ_D + ⟨^T _e _D , 1 ⟩_∂Ω_e∩Γ_D = 0   for e=1,… ,


and it is utilized to close the global problem.





 §.§ Weak form of the local and global problems




Consider the following discrete functional spaces according to the notation introduced in <cit.>:


	
    (Ω)     := { v ∈(Ω) : v |_Ω_e∈(Ω_e)  ∀Ω_e   ,   e=1,… ,} , 
    (S)     := {∈(S) : |_Γ_i∈(Γ_i)	 ∀Γ_i⊂ S⊆Γ∪∂Ω},



where 𝒫^k(Ω_e) and 𝒫^k(Γ_i) are the spaces of polynomial functions of complete degree at most k in Ω_e and on Γ_i, respectively. 

The discrete weak formulation of the local problems in Equation (<ref>) is as follows: for e=1,… ,, given _D on Γ_D and ^h on Γ∪Γ_N, find (^h_e ,^h_e,p^h_e) ∈ [(Ω_e)]^× [(Ω_e)]^×(Ω_e) such that


	
    - (,^h_e)_Ω_e + (^T    , ^h_e)_Ω_e = 
       ⟨_e^T  , _D⟩_∂Ω_e∩Γ_D + ⟨_e^T  , ^h ⟩_∂Ω_e∖Γ_D , 
    	
       
    
    
    	
    		-(, ^h_e)_Ω_e  - (^T    , p^h_e)_Ω_e
        + ⟨ ,_e^T ( ^h_e +   p^h_e )⟩_∂Ω_e =  (,s)_Ω_e , 
    	
       
    	(^T   q, ^h_e)_Ω_e  = ⟨ q , ^T _e _D ⟩_∂Ω_e∩Γ_D + ⟨ q , ^T _e ^h ⟩_∂Ω_e∖Γ_D , 
       1/|∂Ω_e|⟨ p^h_e, 1 ⟩_∂Ω_e = ρ^h_e ,



for all ( ,,q) ∈ [(Ω_e)]^× [(Ω_e)]^×(Ω_e).

The trace of the numerical normal flux in Equation (<ref>) is defined as follows


    _e^T (^h_e +   p^h_e ) := 
    
    _e^T ( ^h_e +   p^h_e ) +  (^h_e - _D)    on ∂Ω_e∩Γ_D,
    
    _e^T ( ^h_e +   p^h_e ) +  (^h_e - ^h)    elsewhere,


where the stabilization parameter  plays a crucial role in the stability, accuracy and convergence properties of the resulting HDG method <cit.>.

By plugging Equation (<ref>) into Equation (<ref>) and integrating by parts, the symmetric form of the discrete weak local problem is obtained: for e=1,… ,, given _D on Γ_D and ^h on Γ∪Γ_N, find (^h_e ,^h_e,p^h_e) ∈ [(Ω_e)]^× [(Ω_e)]^×(Ω_e) that satisfy


	
    - (,^h_e)_Ω_e + (^T    , ^h_e)_Ω_e = 
       ⟨_e^T  , _D⟩_∂Ω_e∩Γ_D + ⟨_e^T  , ^h ⟩_∂Ω_e∖Γ_D , 
    	
       
    
    
    
    	
    	(, ^T ^h_e)_Ω_e +    ⟨, ^h_e ⟩_∂Ω_e + (, ^T   p^h_e)_Ω_e = 
        (,s)_Ω_e + ⟨ , _D⟩_∂Ω_e∩Γ_D + ⟨ , ^h ⟩_∂Ω_e∖Γ_D,
    	
       
    	(^T   q, ^h_e)_Ω_e  = ⟨ q , ^T _e _D ⟩_∂Ω_e∩Γ_D + ⟨ q , ^T _e ^h ⟩_∂Ω_e∖Γ_D , 
       1/|∂Ω_e|⟨ p^h_e, 1 ⟩_∂Ω_e = ρ^h_e ,



for all ( ,,q) ∈ [(Ω_e)]^× [(Ω_e)]^×(Ω_e).


From a practical point of view, the constraint on the mean value of the pressure on the boundary of the element introduced in Equation (<ref>) is handled by means of a Lagrange multiplier.

Thus, the matrix associated with the resulting local problem has a saddle point structure <cit.>.


For the global problem, the discrete weak formulation equivalent to (<ref>) is: find ^h∈[(Γ∪Γ_N)]^ and ρ^h ∈^ such that

 

    ∑_e=1^{⟨, _e^T ^h_e ⟩_∂Ω_e∖Γ_D    + ⟨, ^T _e p^h_e ⟩_∂Ω_e∖Γ_D 
    + ⟨, ^h_e ⟩_∂Ω_e∖Γ_D. 
       . 
    - ⟨, ^h ⟩_∂Ω_e∖Γ_D}
    = -∑_e=1^⟨, ⟩_∂Ω_e∩Γ_N,
    
       
    ⟨^T _e  , 1 ⟩_∂Ω_e ∖Γ_D = - ⟨^T _e _D , 1 ⟩_∂Ω_e∩Γ_D = 0   for e=1,…,  ,



for all ∈[(Γ∪Γ_N)]^.




































































 §.§ Local post-process of the velocity field




As usual in HDG, an element-by-element post-process procedure is considered to construct an improved approximation of the velocity field.

Modifying the Brezzi-Douglas-Marini (BDM) projection operator <cit.>, in <cit.>, a technique to retrieve an -conforming and exactly divergence-free velocity field was discussed.

In this section, a simpler approach inspired by the work of Stenberg <cit.> and exploited in <cit.> is considered.

The requirement of -conformity is relaxed and the resulting local post-process problem exploits the optimal convergence rate of order k+1 of the mixed variable to construct a velocity field ^⋆ superconverging with order k+2.

Nevertheless, it is known <cit.> that using the Cauchy formulation of the Stokes equation, a loss of superconvergence is experienced by low-order approximations.

The Voigt notation introduced in Section <ref> allows to remedy this issue and to circumvent the complex mathematical framework of M-decomposition discussed in <cit.> to devise superconvergent HDG approximations with strongly and weakly symmetric stress tensors.

Following <cit.>, the space _⋆(Ω) of the polynomials of complete degree at most k+1 on each element Ω_e


    _⋆(Ω) := { v ∈(Ω) : v |_Ω_e∈𝒫^k+1(Ω_e)  ∀Ω_e   ,   e=1,… ,}


is introduced.

For each element Ω_e,  e=1,…,, the post-processed velocity ^⋆ is the solution of the problem


    {^T _e^⋆   = - ^T ^h_e             in Ω_e,
    _e^T _e^⋆   = - _e^T ^h_e              on ∂Ω_e,
    .


in the space [ _⋆(Ω) ]^.

The element-by-element problem in Equation (<ref>) is obtained by the definition of the mixed variable in Equation (<ref>) and exploits the naturally equilibrated fluxes as condition on the boundary of the element.

The solution of Equation (<ref>) is determined up to rigid motions, namely  translations and  rotations, being = 2 and = 1 in 2D and =  = 3 in 3D.

According to <cit.>, a set of + constraints is introduced to retrieve the uniqueness of the solution.

On the one hand, the indeterminacy due to the  rigid translational modes is resolved introducing the following constraint on the mean value of the velocity:


    ( ^⋆_e , 1 )_Ω_e = ( ^h_e , 1)_Ω_e .
 

On the other hand, the  rigid rotational modes are taken care of by means of a condition on the  of the velocity, namely


    ( ×^⋆_e , 1 )_Ω_e  =  ⟨^h · , 1 ⟩_∂Ω_e ,


where the right-hand side of Equation (<ref>) follows from the application of Stokes' theorem, being ^h_e = ^h on ∂Ω_e and  the tangential direction to the boundary ∂Ω_e.

It is worth noting that other conditions may be considered to resolve the indeterminacy of the problem in Equation (<ref>).

Nevertheless, in order for the post-processed velocity to be superconvergent, the quantities appearing in these constraints have to converge with order m ≥ k+2.

If this is not the case, despite the resulting system admits a unique solution, the superconvergence property is lost.

For the strategy discussed in the present work, extensive numerical experiments have shown that the right-hand sides of both (<ref>) and (<ref>) converge with order m > k+2.

A rigorous proof of this result is currently under investigation.


Recall that the  of the velocity represents the vorticity of the fluid.

Within this context, the left hand side of Equation (<ref>) may be physically interpreted as the mean value of the vorticity inside the element Ω_e.

Similarly, the right-hand side represents the circulation of the flow around the boundary ∂Ω_e.


Eventually, by exploiting the Voigt notation, Equation (<ref>) is equivalent to


    ( ^⋆_e , 1 )_Ω_e = ⟨[ ^h ]^T  , 1 ⟩_∂Ω_e .






§ NUMERICAL STUDIES




In this section, several examples with known analytical solution are considered, in two and three dimensions, to verify the optimal convergence and superconvergence properties of the error of the primal, mixed and post-processed variables, measured in the (Ω) norm and for different element types.

First, a numerical study of the influence of the stabilization parameter  on the accuracy of the proposed HDG method is performed.





 §.§ Influence of the stabilization parameter




As previously stated and extensively studied in a series of publications by Cockburn and co-workers (cf. e.g. <cit.>), the HDG stabilization parameter has an important effect on the convergence properties of the method.

For the sake of simplicity, a stabilization tensor of the form = τ, equal on all the faces of the internal skeleton Γ∪Γ_N is considered.

In what follows, a numerical study of the role of the scalar parameter τ is presented.





  §.§.§ Two dimensional example




The first example considers the well-known problem of the Wang flow in the domain Ω = [0,1]^2. 

The source term s is selected so that the analytical velocity field has the following expression


    () =
       
      2 a x_2 - b λcos(λ x_1) exp{-λ x_2}
     
      b λsin(λ x_1) exp{-λ x_2} ,


whereas the pressure is uniformly zero in the domain.

The values a = b = λ = 1 are set for the constants and the kinematic viscosity ν is taken equal to 1.

Neumann boundary conditions, corresponding to the analytical normal flux, are imposed on Γ_N = {(x_1,x_2) ∈Ω  |   x_2=0} and the analytical velocity field is enforced on Γ_D = ∂Ω∖Γ_N via Dirichlet boundary conditions.

Uniform meshes of quadrilateral and triangular elements are considered. The second level of refinement of the meshes is shown in Figure <ref>.



It is worth noting that the triangular mesh #1 has considerably more degrees of freedom than the triangular mesh #2 for a similar characteristic size. 

The components of the velocity field computed on the fourth level of refinement of the triangular mesh #2 and using a quadratic degree of approximation are depicted in Figure <ref>.



Figure <ref> shows the evolution of the error of the primal, mixed and post-processed variables, , p,  and ^⋆, in the (Ω) norm as a function of the stabilization parameter τ.

The numerical study is performed on the fourth level of mesh refinement, using polynomial approximations of complete degree 1 and 2 and values of τ spanning from 0.1 to 10,000.



It is straightforward to observe that for all the meshes under analysis, there exists a value of τ minimizing the (Ω) norm of the error of the velocity.

Nevertheless, to guarantee the accuracy of the approximation, the (Ω) norm of the error should be accounted for and consequently both  and  are considered in the choice of the optimal value of τ.

Within this context and in order for the post-processed velocity field ^⋆ to provide a gain in accuracy with respect to , the value τ = 4 is chosen for quadrilateral meshes and triangular meshes of the first type.

For triangular meshes of the second type, the minimum of the error in the primal variable is achieved for values of τ substantially larger than 10.

Despite the approximation of the mixed variable deteriorates when the stabilization parameter increases, this effect is limited for values of τ < 50.

The value of τ=40 is thus considered as it provides a good compromise for the quality of the approximation of the primal, mixed and post-processed variables.


Consider the family of meshes in Figure <ref>.

The triangular mesh #1 features one node located in the barycenter of each underlying quadrilateral. The resulting mesh provides significantly more information than the triangular mesh #2 of the corresponding refinement level.

Thus, owing to the aforementioned extra node and to the tensorial nature of the basis functions defined on the quadrilateral meshes, the behavior of the triangular meshes #1 is expected to be more similar to the quadrilateral ones than to the triangular meshes #2, as observed in the previous numerical simulations in Figure <ref>.






  §.§.§ Three dimensional example




The second example, inspired by <cit.>, is an analytical solution of the problem in Equation (<ref>) set in the domain Ω = [0,1]^3. 

The source term is selected so that the analytical velocity is 


    () =
       
      b exp{a(x_1-x_3) + b(x_2-x_3)} - a exp{ a(x_3-x_2) + b(x_1-x_2)}
     
      b exp{a(x_2-x_1) + b(x_3-x_1)} - a exp{a(x_1-x_3) + b(x_2-x_3)}
    
      b exp{a(x_3-x_2) + b(x_1-x_2)} - a exp{a(x_2-x_1) + b(x_3-x_1)}


and the corresponding pressure field is


    p() = x_1 (1-x_1) .


The values a = 1 and b = 0.5 are considered and the kinematic viscosity ν is taken equal to 1.

Neumann boundary conditions, corresponding to the analytical flux, are imposed on Γ_N = {(x_1,x_2,x_3) ∈Ω  |   x_3=0} and the analytical velocity field is enforced on Γ_D = ∂Ω∖Γ_N via Dirichlet boundary conditions.

Figure <ref> shows a cut through the third level of refinement of the uniform meshes of hexahedral, tetrahedral, prismatic and pyramidal elements considered in this study.



The velocity and pressure fields computed on the third level of refinement of the hexahedral mesh and using a cubic degree of approximation are depicted in Figure <ref>.



The evolution of the error of the primal, mixed and post-processed variables, , p,  and ^⋆, in the (Ω) norm as a function of the stabilization parameter τ is presented in Figure <ref>.



As highlighted by the theory <cit.> and confirmed by the analysis of the two dimensional case, a value of the stabilization parameter of order one (i.e. τ∈ [1,10)) guarantees stability and convergence of the HDG method.

More precisely, a value near τ=10 provides the minimum error for the primal variable but limited or no extra gain in accuracy is obtained through the post-process of the velocity field.

Thus, a value of τ=4 is selected for the following simulations.

The discussed numerical results show that the HDG discretization is robust to the choice of the stabilization parameter.

Moreover, the optimal value of τ is not dependent upon the degree of approximation or the dimensionality of the problem.

Considering the different types of elements under analysis, the triangular meshes #2 require a slightly larger value of the stabilization parameter to enter the asymptotic regime and show the optimal convergence and superconvergence properties expected from the theory.





 §.§ Optimal convergence and superconvergence of the primal, mixed and post-processed variables




Consider the optimal values of τ identified in the previous section.

The optimal convergence properties of the velocity , the pressure p and the mixed variable  representing the strain rate tensor, are tested for different element types using the (Ω) norm.

Moreover, the superconvergence of the post-processed velocity field ^⋆ is also analyzed.





  §.§.§ Two dimensional example




In Figure <ref>, the first column presents the convergence of the error of the primal and mixed variables p and , measured in the (Ω) norm, as a function of the characteristic element size h for both quadrilateral and triangular elements and for a degree of approximation ranging from k=1 up to k=3. 

In a similar fashion, the second column provides the corresponding convergence history for the primal and the post-processed velocities  and ^⋆.



It can be observed that almost the optimal or the optimal rate of convergence h^k+1 is obtained for , p and , for all the element types and degrees of approximation considered. 

As previously mentioned, the triangular mesh #1 has considerably more degrees of freedom than the other meshes: in particular, for the same characteristic element size, the triangular meshes #1 have approximately 2.5 times (respectively, 5 times) more internal faces than the triangular mesh #2 (respectively, the quadrilateral mesh). 

Thus, despite the results in Figure <ref> indicate that the triangular mesh #1 provides more accuracy than the other meshes, a comparison in terms of the global number of degrees of freedom confirms that similar results are obtained using meshes of different element types.

Concerning the post-processed variable, the rate of convergence h^k+2 is achieved and the superconvergence property is verified.

This confirms that the average of the hybrid variable  on the boundary leads to a superconvergent approximation, as observed in <cit.> for the linear elastic problem.

Beside the improved convergence rate, the discussed post-process procedure is responsible for a gain in accuracy of ^⋆ with respect to the original approximation  of the velocity field.

Hence, the information encapsulated in the primal and post-processed variables may be exploited to construct an error indicator and devise an automatic degree adaptivity strategy as discussed in <cit.>.





  §.§.§ Three dimensional example




Similarly to the previous example, the convergence of the error of p and  (Fig. <ref>) and  and ^⋆ (Fig. <ref>), measured in the (Ω) norm, as a function of the characteristic element size h is presented for hexahedral, tetrahedral, prismatic and pyramidal elements and for a degree of approximation ranging from k=1 up to k=3. 





As for the two dimensional case, almost the optimal or the optimal rate of convergence h^k+1 is obtained for , p and  in 3D, for all the element types and degrees of approximation considered (cf. Fig. <ref>-<ref>). 

In Figure <ref>, the post-processed variable is shown to superconverge with a rate of convergence h^k+2.

Beside the improved convergence rate, the discussed post-process procedure is responsible of a gain in accuracy of ^⋆ with respect to the original approximation  of the velocity field.

The presented numerical experiments in two and three dimensions confirm that exploiting Voigt notation the HDG approximation of the Stokes equation achieves optimal convergence rate h^k+1 for both the primal variables  and p and the mixed one . 

In particular, contrary to what observed in <cit.>, the convergence of the mixed variable does not deteriorate when considering the Cauchy formulation of the Stokes flow.

As discussed in <cit.> for the linear elastic problem, the post-process technique exploiting the  of  allows to construct an approximation of the primal vector field superconverging with order k+2.

Moreover, the post-process strategy provides an extra gain in accuracy with respect to the original approximation of the velocity field.

As highlighted in Figure <ref>, a solution that is almost one order of magnitude more precise than the HDG solution is obtained, even for linear approximations.





 §.§ Numerical evaluation of quantities of interest: drag force on a sphere




The last example considers the classical test case of the viscous flow around a sphere. 

The objective of this test is to show the capability of the described HDG method to provide an approximation of the pressure and the viscous forces sufficiently accurate to evaluate a quantity of interest with the precision required by industrial standards.

Consider the domain Ω = ( [-7,15] × [-5,5] × [-5,5] )∖ℬ_1,0, ℬ_1,0 being a ball of unit radius centered at the origin.



To reduce the computational effort, the symmetry of Ω is exploited and solely one fourth of the domain is taken into account to perform the numerical experiments.

Different tetrahedral meshes of the domain are considered, ranging from 3,107 to 204,099 elements.

High-order computations employ isoparametric curved meshes. The extension to high-order is performed using the solid mechanics analogy described in <cit.>.

Figure <ref> (a)-(b) shows the magnitude of the velocity with streamlines of the flow and the pressure field computed on the third level of refinement of the mesh, featuring 43,682 tetrahedrons, and using a quadratic degree of approximation.

The results in Figure <ref> (c) show the convergence of the drag force as the number of degrees of freedom is increased, i.e. for different levels of mesh refinement and for a degree of approximation ranging from k=1 up to k=3. 

The numerically computed drag is compared with the analytical value from the literature <cit.>.

In Table <ref>, a quantitative analysis of the relative error in the computation of the drag force is reported for all the mesh refinements and degrees of approximation considered. 

Using linear elements, almost 4 millions degrees of freedom are required by the method to compute the drag coefficient with a relative error of 2%.

The same level of accuracy is achieved by quadratic and cubic elements using the coarsest mesh under analysis and less than 200,000 degrees of freedom.

More precisely, moving to high-order approximations, errors lower than 0.5% are obtained using few hundreds thousands degrees of freedom.

The observed additional accuracy results from the concurrent use of high-order polynomial functions for the discretization of the unknown variables and high-order approximations of the geometry via meshes featuring curved elements.

Thus, the superiority of high-order methods with respect to low-order ones discussed in the literature (cf. e.g. <cit.>) is confirmed.







§ CONCLUSION




This paper describes a hybridizable discontinuous Galerkin method using Voigt notation, first introduced in <cit.>, for the Cauchy formulation of the Stokes equation. 

Owing to Voigt notation, the symmetry of the stress tensor is strongly enforced by storing in a vector format only half of the off-diagonal terms.

Moreover, physically meaningful tractions may be naturally imposed on the Neumann boundary.

Contrary to the existing superconvergent HDG formulations involving the symmetric part of the gradient, the proposed method does not enrich the discrete spaces of approximation and it reduces the number of degrees of freedom of the mixed variable.

Hence, the resulting local problems are smaller and computationally more efficient.

The optimal convergence order k+1 is achieved for all the unknowns, as proved for the classical HDG equal-order approximation of the velocity-pressure formulation and for the more involved discretization of the Cauchy formulation based on the M-decomposition.

The novelty and main advantage of the present approach relies on being able to exploit the same degree of approximation for both primal and mixed variables, in presence of the symmetric part of the gradient.

In addition, a velocity field superconverging with order k+2 is obtained via a local post-process procedure, exploiting the optimal convergence of the mixed and hybrid variables.



Numerical studies show the optimal convergence and superconvergence properties of the method in 2D and 3D using meshes of different element types and the robustness of the approach with respect to the choice of the HDG stabilization parameter.

Eventually, the drag force on a sphere is evaluated using different degrees of approximations to show the capability of the method to compute industrially relevant quantities of interest with an acceptable precision.





§ ACKNOWLEDGEMENTS



This work was partially supported by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska–Curie grant agreement No. 675919 and the Spanish Ministry of Economy and Competitiveness (Grant number: DPI2017-85139-C2-2-R).

The support of the Generalitat de Catalunya (Grant number: 2017SGR1278) is also gratefully acknowledged.

Finally, Alexandros Karkoulias was supported by the European Education, Audiovisual and Culture Executive Agency (EACEA) under the Erasmus Mundus Joint Doctorate Simulation in Engineering and Entrepreneurship Development (SEED), FPA 2013-0043.




plain



