Assessment of Nitsche’s Method for Dirichlet Boundary Conditions Treatment

One of the big advantages of the standard finite element method is its efficiency in treating complicated geometries and imposing the associated boundary conditions. However in some cases, such as handling the Dirichlet-type boundary conditions, the stability and the accuracy of FEM are seriously compromised. In this work, Nitsche’s method is introduced, as an efficient way of expressing the Dirichlet boundary conditions in the weak formulation. It is shown that Nitsche’s method preserves the rate of convergence and gives more accuracy than the classical approach. The method is implemented for the simplest case of Poisson equation, for Stokes flow and Navier-Stokes equations, with slip and noslip boundary conditions, in the case of viscous Newtonian incompressible flows. Error norms are calculated on different meshes in terms of size, topology and adaptivity, for a fair assessment of the proposed Nitsche’s method.


Introduction
The Nitsche's method was introduced by [1] to capture the essential boundary conditions in FEM and express it in particular weak form.Nitsche's method is a particular case of Lagrange multiplier method [2], [3] and [4].We know that imposing boundary conditions is a key issue in study solutions of problems not only in fluid dynamics but in all scientific problems.
The original idea with this method was handling the Dirichlet boundary conditions weakly without using Lagrange Multipliers.We know that finite element can treat very complicated geometries with different mesh shapes; however, it has some restrictions in few cases.Like in fluid /structure interaction FSI where we have the slip/no slip boundary conditions for viscous incompressible flow.It was shown that the method will converge optimally if the finite element space satisfies the "sup-inf" condition [5], [6].In this regard, Nitsche's method has resurgence in recent years.The idea behind the Nitsche's method is to simply replace the Lagrange Multipliers arising in a dual formulation, through their physical representation.Nitsche [1] also added an additional penalty like term to restore the coercivity of bilinear form.The Nitsche's method is employed in major domains of science: In electro-hydrodynamics potential of flows [7], in plasticity [8], solid deformations [9] and other domains where boundary conditions are crucial.

Poisson equation
To introduce the Nitsche's method, let us consider the Poisson equation with Dirichlet boundary condition Where Ω is a bounded domain in two dimensional space, f is a given function and  0 is the value of u on the boundary ∂Ω.
The weak formulation of this problem with strong imposed Dirichlet boundary conditions is to find  ℎ ∈  ℎ such that: Where ( , ) Ω denote the  2 (Ω) scalar product The Nitsche's method for Poisson equation is formulated as follow: Find  ℎ ∈  ℎ ⊂  1 ( Ω) such that The classical (symmetric) Nitsche's formulation of the problem is: The method is consistent and has an optimal error estimate in the norm [10].
More details about the convergence of the method we refer to [10].Asymmetric Nitsche's formulation With  ∂Ω > 0 and constant, it's called stabilization parameter.In our case we are going to take  ∂Ω = 2,  ⃗ is the outward unit normal vector to the boundary ∂Ω.
For the numerical example let consider Ω is a bounded domain in two dimensional computational space, Ω ⊂ ℝ 2 , with boundary Ω .The normal vector to Ω is denoted by  ⃗ .In this example we will consider a unite square domain [0,1] Χ [0,1], with  ×  cells.We are going to take the number of partitions N variable in each test.And compare each numerical result obtained from the Nitsche's method with the exact analytical solution of the problem and calculate the  2 error norms in each case.The classical FEM weak formulation is coupled with the Nitsche's method.
The exact solution of the Poisson problem in 2D Cartesian space (x,y) is taken as follow Let's define the errors and norms, which ¸are equivalent in the symmetrical case and the unsymmetrical case.We plot in fig 1(a), the variation of the log (e) and log (h) with increasing the number of cells N. We can see that both are decreasing with increasing the number of cells.
Fig. 1 (b) represents the grid convergence, which estimates the order of the method numerically (the order is the slop of the curve), as we can see that the added Nitsche's terms don't degrade numerically the order of the method.In Fig. 1 (c), we represent the variation of log(e) with log(h) for different values of the stabilization parameter  ∂Ω .As we can see the curves have the same shape and more or less the same values.

The Stokes equation
Let Ω be a bounded domain in ℝ 2 with a smooth boundary Ω.The stationary Stokes equation for a viscous incompressible flow, in dimensionless variables, can be written as: With  is the velocity,  is the dynamic viscosity,  the pressure and  the external forces acting on the fluid.The 2D exact solution of the Stokes problem can be taken We are going to solve the Stokes problem, associated to Dirichlet Boundary conditions, numerically using the finite element method FEM, coupled with the Nitsche's method for expressing the Dirichlet boundary conditions in Nitsche's form.After that we compare the numerical results with the exact solution of the problem and calculate the error between them.
Let's take first the weak formulation of the problem, with a viscosity equal to the unity.Dirichlet conditions are now enforced weakly using Nitsche's method by penalty-like terms and additional terms that maintain the skew-symmetry of the Stokes operator.
We use the Taylor Hood [11]setting and choose biquadratic finite element functions for the discrete velocity space and bilinear finite element functions for the discrete pressure space .Adding the skew symmetrisation, the penalty term  1 ℎ   < ,  > ∂Ω +  2 ℎ  < ., . > ∂Ω and the inflow stabilization [11].With  1 and  2 are the stabilization parameters.
For the simulation a rectangle domain with quadratic mesh is taken with length of 20 and width of 4 , [0,20] × [0,4], the number of cells  ×  is going to be different in each case.Reynolds number will be constant during the simulation and will be equal to  = 300.
Symmetric formulation for the volumetric part For this part we add terms forcing Dirichlet boundary conditions and then obtaining the Nitsche's symmetric formulation of our problem.176-5 With ( ⃗ , ) trial functions and ( , ) test functions of our problem, with  is the constant of parameterization and it's taken equal to  = 10.ℎ  denote the local mesh size and u  the exact solution of the Stokes equation.
Skew-Symmetric formulation is: We couple the classical formulation and the Nitsche's formulation of the Stokes problem in the same running program.We will use the well-known stabilizing technique [12] to handle the divergence constraint.The slip condition  ⃗ . ⃗ = 0 is taken on the wall boundary of ∂Ω.This seems to be more practical than imposing it in strong form as that would involve both local coordinate transformations and the approximation of the normal at corners.
As we can see fig.
(2) represents the profile of velocity over vertical section in the middle of the domain and correspond to the shape of the given exact solution   .We can see also that the free slip condition over the wall boundaries is respected and the value of the tangential velocity is not equal to zero.
In fig (3) we can see the increasing value of velocity when the flow is involving in the horizontal direction.And the graphic (4) represents the increasing value of the velocity in the central line of the domain.
Fig. 5 represents the graphic of pressure in a vertical section in the middle of the domain, the curve is in saddle shape.We summarize in Table 2 the variation of the error and the rate related to the value of N and the local mesh size h.We obtain the following graphics Figs. 6 (a) which represent the variation of the logarithmic value of the error and the logarithmic value of the mesh size in function of number of cells N.Then, we can see that the value of log(e) and log(h) are decreasing with increasing the number of cells N, which is first a confirmation of the fact that the Nitsche's method is a mesh dependent method.Secondly, the error between the Nitsche's method and the exact solution of the problem is really small and starts to be negligible if we increase the number of cells sufficiently.
In Fig. 6 (b) we give the representation of the variation of log(e) with log(h) for different values of the stabilization parameter .
The internal stress tensor is defined as: It is composed from the pressure () and a product of the dynamic viscosity denoted as  and the strain rate tensor () given by As Boundary conditions of the problem we take the inlet pressure constant and equal to   = 1000  .The inlet velocity is  ⃗  = (sin() , 0) .
And  is the external volumetric forces acting on the fluid and its taken equal to 1, ‖ ‖ = 1 .The density is taken  = 1000/ 3 and the dynamic viscosity  = 1// Our domain Ω is a rectangle with quadratic mesh and it is taken with a length of 10 and a width of 2, Ω = [0,10] × [0,2] and a number of cells equal to 2500.
In this section we are going to study the free slip boundary condition and the no slip boundary condition near to the wall boundary of the stationary Navier Stokes equation for an incompressible flow, using the Nitsche's formulation for the Dirichlet boundary condition in each case.

2-3-A-Free slip case
Free slip case correspond to:  ⃗⃗⃗ . ⃗ = 0 and  ⃗⃗⃗ . ≠ 0 at the wall boundary.It seems the normal component of the velocity is equal to zero and the tangential component of the velocity is not equal to zero.We have kind of a drift flow at the wall boundary.
Nitsche's formulation of the Dirichlet boundary conditions for the Navier Stokes equation in the free slip case is: Volumetric part To the volumetric terms we need to add Nitsche's terms for all the Dirichlet boundaries of the domain with imposing the condition of free slip.Table 3 gives us different terms added to the bilinear and linear form for each boundary nature (Inlet, outlet and wall) of the domain to obtain the final Nitsche's formulation of the stationary Navier stokes equation.Terms added to a(u,v) Terms added to L(v)

2-3-B-No Slip Case
The no slip case corresponds to:  ⃗⃗⃗ . ⃗ = 0 and  ⃗⃗⃗ . = 0 at the wall boundary.We impose the zero velocity for both normal and tangential components.The Nitsche's formulation of the Dirichlet boundary conditions will be different from the previous case.
Nitsche's formulation for the no slip case: The following table (Table 4) gives us different terms added to the bilinear and linear form for each boundary nature (Inlet, Outlet and Wall) of the domain to obtain the final Nitsche's formulation of the stationary Navier-Stokes equation in the no slip case.

Conclusion
In this work, we gave a presentation of the Nitsche's method and its application in different cases starting by the simple Poisson equation then the Stokes flow and finally the Navier-Stokes equation in both cases: slip and no slip boundary conditions.But, this is by no means the only field of application for the approach.
The Dirichlet boundary conditions were enforced weakly by using the Nitsche's method.The error was calculated and we demonstrate through several numerical examples the performance and the robustness of the proposed variation form.The Nitsche's method has a very wide range of application and it's a very physical method.It employs only continuity conditions regarding the primal variable and the fluxes and simplify the variable coupling.
The only problem of the Nitsche's method is to find the weak form.The generalization of the implementation for other problems is not as straightforward as for Lagrange multiplier method.
The choice of the stabilization parameters depends not only on the partial differential equation but also on the essential boundary condition.The choice of the stabilization parameters needs to be carefully made.This condition ensures the coercivity of the bilinear form in the interpolation space.Even in the case where there are not chosen rightly, the convergence is assured with a certain degree of convergence not far from the exact solution.
Regarding the choice of parameters, Nitsche's method proved that if these parameters are taken with specific choice, then the discrete solution converges to the exact solution with optimal order.The advantage of this method is evident when we take into account the local mesh sizes.
An extension of the work will be done for the three dimensional problems, and implementation challenges are associated with the computational geometry and algorithmic aspects of implementation.

Fig. 1 :
(a) Representation of the variation of log(h) and log (e) in function of N, (b) Convergence plot for the method, (c) Variation of log(e) for different values of the stabilization parameter in function of log(h).

Fig. 2 :
Fig. 2: Velocity profile in vertical section at x=10 (middle of the domain).

Fig. 7 (
Fig. 7 (a) represents the graphic of velocity, as we can see the velocity is not equal to zero near to the wall boundaries as imposed.Fig.7(b) shows the variation of pressure in the center line along the geometry and the Figs. 8 (a) represent the distribution of velocity along the geometry, and the values are matching to the previous graphic.

Fig. 8 :
Distribution of velocity and pressure in the geometry.

Figs. 8
Figs. 8 (b) represent the value of the pressure in the complete geometry.Negative values of the pressure represent the recirculation zones were we have tourbillons.

{
= ( ∇()., ) Ω + ( ∇(), ∇()) Ω − ( p, ∇. (v)) Ω − ( q, ∇. (u)) Ω − ( ∇()., v) Ω  = (,  ) Ω (17) As we can see, Fig.9(a) represent the distribution of velocity in the domain, the value of the velocity at the wall is equal to zero, we can also see the important value of the velocity in opposite imposed flow direction near to the wall due to the recirculation of the flow in the opposite direction, we can see more details about the velocity direction in the Fig.9(b) .Fig. 10(a) represents a plot of velocity magnitude in a vertical section.Fig. 10(b) is a graphical representation of pressure values in a center line of the geometry (a) (b) Fig. 9: (a) Value of velocity in the no slip condition case, (b) Recirculation zone near to the wall boundary.: (a) Velocity distribution in vertical section at the middle of the geometry x=5,(b) Pressure in center line of the domain.

Table 1 .
In each case, we have the result for h, the error 'e', the rate and the two dimensional representation of  ℎ .After running the simulation, we obtain numerical results in Table1.Values Of H, Error, The Rate And The Solution  ℎ For Different Values Of N.

Table 2 :
Values of h, Error and the Rate for Different Values of N. We consider a two dimensional computational domain Ω ⊂ ℝ 2 with a boundary Ω .The boundary is composed only on Dirichlet boundary condition.The fluid velocity () and the pressure () of the flow problem is governed by the incompressible and isothermal Navier-Stokes equation.The Navier-Stokes equation and the mass conservation equation for an incompressible and isothermal flow are:

Table 3 :
Nitsche's Terms Added to the Linear and Bilinear Parts of Navier-Stokes Equation in Free Slip Case.

Table 4 :
Nitsche's Terms Added to the Linear and Bilinear Parts of Navier-Stokes Equation in the No Slip Case.