A benchmark problem to evaluate implementational issues for three-dimensional flows of incompressible fluids subject to slip boundary conditions

We consider ﬂows of an incompressible Navier–Stokes ﬂuid in a tubular domain with Navier’s slip boundary condition imposed on the impermeable wall. We focus on several implementational issues associated with this type of boundary conditions within the framework of the standard Taylor-Hood mixed ﬁnite element method and present the computational results for ﬂows in a tubular domain of ﬁnite length with one inlet and one outlet. In particular, we present the details regarding variants of the Nitsche method concerning the incorporation of the impermeability condition on the wall. We also ﬁnd that the manner in which the normal to the boundary is numerically implemented inﬂuences the nature of the computational results. As a benchmark, we set up steady ﬂows in a tube of ﬁnite length and compare the computational results with the analytical solutions. Finally, we identify various quantities of interest, such as the dissipation, wall shear stress, vorticity, pressure drop, and provide their precise mathematical deﬁnitions. We document how well these quantities are computationally approximated in the case of the benchmark. Although the geometry of the benchmark is simple, the correct computational results require careful selection of numerical methods and surprisingly non-trivial computational resources. Our goal is to test, using the setting with a known analytical solution, a robust computational tool that would be suitable for computations on real complex geometries that have relevance to problems in engineering and medicine. The model parameters in our computations are chosen based on ﬂows in large arteries.


Introduction
Most fluids that can be described by the Navier-Stokes model are assumed to meet the no-slip boundary condition, that is adherence of the fluid to the boundary, in flows that take place in pipes, channels, and other simple domains when the flow is reasonably slow. The application of this particular boundary condition is attributed to Stokes who is supposed to have advocated its use. It is not really clear that Stokes believed in its general validity. He even had doubts concerning its accuracy in pipes and channels, as he remarked "Du Buat found by experiment that when the mean velocity of water flowing through a pipe is less than about one inch in a second, the water near the inner surface of the pipe is at rest. If these experiments may be trusted, the conditions to be sat-these materials. This is definitely a formidable task. For materials such as blood (which is a complex mixture of plasma 1 , red and white blood cells, platelets, lipo-proteins, and a variety of complex ions) and blood vessel (which has a very complex layered structure) it is far from clear that the no-slip boundary condition ought to hold, see ( Rajagopal and Rajagopal, 2020 ) for more details.
Experiments concerning the flow of blood in capillaries by Bennett (1967) and Bugliarello and Hayden (1963) report the slipping of red blood cells that come into contact with the boundary. In Nubar (1967) , Nubar (1971) , Misra and Shit (2007) , the possible connection of slippage on the wall and shear rate dependent blood viscosity measured in the rheometer is suggested. The above flows are reasonably slow flows and even in such flows certain components of blood seem to be slipping at the boundary. Casting further doubt concerning the appropriateness of the no-slip boundary condition is the fact that the flow of blood can be turbulent at some locations in the cardiovascular system and far from the slowness presumed by Stokes when the no-slip boundary condition is experimentally observed. Thus, it would be worthwhile to study the flow of a material like blood by assuming the slip boundary condition.
Numerous types of boundary conditions that go beyond the no-slip boundary condition were proposed by the early pioneers in the field of fluid dynamics ( de Coulomb, 1800;Du Buat, 1786;Girard, 1813;Navier, 1823;de Prony, 1804 ), and Poisson (1831) . See ( Neto et al., 2005 ) for a more recent overview of experimental studies on boundary slip measurements.
In this study, we focus on Navier's slip boundary condition (proposed by Navier, 1823 in agreement with molecular arguments) characterized by the parameter and involving, as the limiting cases, no-slip boundary condition at one end ( = 1 ) and (perfect) slip boundary condition on the other end ( = 0 ). This extended class of boundary conditions, taking into account the slipping mechanisms of various degree, represents the most commonly used boundary conditions imposed on the impermeable parts of the boundary of the flow domain. Interestingly, it was shown (see ( Hron et al., 2008 )) that for both plane and threedimensional Poiseuille flow varying the slip parameter can influence the flow more markedly than the change of the constitutive equation for the fluid.
We consider flows of an incompressible Navier-Stokes fluid in a special three-dimensional domain, namely a tubular domain of finite length with one inlet and one outlet together with the above described class of slip boundary conditions imposed on the impermeable wall. The reason for choosing this special geometrical setup is twofold. First, knowing the analytical solutions for steady flows subject to different slippage at the wall of cylinder, we are able to set a benchmark that can serve as a basic test for various numerical methods and their implementations. Second, this geometrical setting can be viewed as a very simplified, yet reasonable approximation of flows in a large vessels, applicable for instance when investigating/assessing flow through aortic valve (native or its replacement), or a pathologically narrowed segment of large artery (e.g. coarctation of aorta or pulmonary artery stenosis).
Next, we present energetic considerations that might provide support for the choice of boundary conditions that we enforce at a rigid boundary and at the outlet and provide the function spaces in which the weak formulations of the problem is formulated. These formulations differ depending on how the incompressibility constraint and the impermeability condition are incorporated. In particular, we present the details regarding variants of the Nitsche method ( Nitsche, 1971 ) concerning (in this study) the incorporation of the impermeability condition on the wall. We provide a literature overview concerning the Nitsche method at the end of Section 3 .
Then, we focus on several issues associated with the implementation of this type of boundary conditions within the framework of the stan-1 Plasma is oftentimes described by a Navier-Stokes fluid while flowing in large blood vessels. dard Taylor-Hood finite elements. We discuss how the manner in which the normal to the boundary is numerically implemented influences the nature of the computational results. We also identify various quantities of interest, such as the dissipation, pressure drop, vorticity, wall shear stress, and provide their precise mathematical definitions. The objective is to document how well these quantities are computationally approximated in the case of the proposed benchmark.
We would like to emphasize that although the geometry of the benchmark is simple, the correct computational results require a careful selection of numerical methods and non-trivial computational resources. Our goal is to test, using the setting with a known analytical solution, a robust computational tool that would be suitable for computations on real complex geometries that have relevance to problems in engineering and medicine. The model parameters in our computations are chosen based on flows in large arteries; this corresponds to flows at Reynolds number of approximately 1050.
The structure of the paper is the following. In Section 2.1 , we formulate the problem describing unsteady flows of incompressible Navier-Stokes fluid in three-dimensional tube with Navier's slip boundary condition on the impermeable wall, a given velocity at the inlet and special boundary conditions at the outlet that takes into account the given normal stress (i.e. the pressure) as well as the control (boundedness) of the energy of the whole system. The energy considerations also lead to the choice of the function spaces in which we state the weak formulation of the problem. In Section 3 , we take a slightly different viewpoint with regard to the weak formulation of the problem and we discuss two variants of the Nitsche method that are used to numerically incorporate, in a weak sense, the impermeability condition. Section 4 provides a brief description of the numerical discretization of weak formulations connected with the Nitsche method. The emphasis is devoted to the way in which the normal vector is implemented. In Section 5 , the benchmark with known analytical solution under the simplifying assumption of the simple shear flow is formulated. We also define two forms of mechanical dissipation of the system and present its connection to the kinetic energy and the pressure drop along the length of the tube. Pressure drop is often used in lieu of the dissipation rate in cases where the real dissipation is unknown. This is clarified in Section 5.2 . Section 6 brings together all the results connected with our numerical approach and computational tools and compares them with the analytical (benchmark) solution. We also compare the efficiency of the variants of the Nitsche method. The concluding Section 7 summarizes the main results that have been established.

Formulation of the initial-boundary-value problem
We consider the flow of an incompressible Navier-Stokes fluid over time interval , in a three-dimensional tube of the finite length representing a bounded, open and connected set Ω ⊂ 3 , see Fig. 1 . The unsteady flows of such fluids described in terms of the velocity = ( 1 , 2 , 3 ) ∶ (0 , ) × Ω → 3 and the mean normal stress − ∶ (0 , ) × Ω → ( being usually called the pressure) are then governed by the following set of equations: where is the Cauchy stress and * > 0 , * > 0 are given constants: the density and the dynamic viscosity. We also use the notation where the incompressibility constraint (1) has been incorporated. Furthermore, for any matrix , the symbol denotes its transpose. The boundary Ω of the tubular domain Ω consists of three nonoverlapping parts Γ in , Γ out and Γ wall and we impose the following initial and boundary conditions involving three given functions: Here ∈ [ 0 , 1 ] , * ∈ (0 , ∞) and is a unit normal vector at the boundary. Furthermore, for ∈ we set + ∶= max {0 , } and − ∶= min {0 , } , i.e. = + + − , and for an arbitrary vector , ∶= ( ⋅ ) is the normal component of and ∶= − is the projection of to the tangent plane.
Note that instead of (6) we could impose the condition By taking the scalar product of (6) with , we observe that (6) and (8) give the identical description in the normal direction. On the other hand, by projecting (6) into the tangent plane, we obtain 1 2 * ( ⋅ ) − = ( ) , which is different than the first condition in (8) . The condition (8) is more appropriate to the situations where the velocity field at the outflow is unidirectional (which is, in fact, the case of the analytical solution used in Section 5). Although the boundary conditions (8) are different than those stated in (6) , in the few computational tests made so far we have not observed any significant difference. This is why the conditions (8) are not discussed anymore below and we consider merely (6) at Γ out in what follows.
The boundary conditions (6),(7) and the function spaces for the weak formulation are in keeping with basic energy estimates that can be established for the problem (1) -(7) , see below. We wish to remark that the presence of the term 1 2 * ( ⋅ ) − in the outflow boundary condition (6) seems to be essential not only for obtaining the energy estimates (see below) but also for the robustness and stability of the numerical methods.

Energy balance and energy estimates
For clarity, we proceed in two steps. First, we simplify the setting by assuming that in = . (9) A general inflow function will be considered later.
To obtain the energy identity, assuming (9) , we form the scalar product of (2) with , use (1) and arrive at ( * where ( ) is the symmetric part of the velocity gradient defined as . Then we integrate this equation over Ω, apply the Gauss theorem, use the impermeability condition on Γ wall and the homogeneous Dirichlet condition on Γ in and conclude that Due to the symmetry of , we have ⋅ = ⋅ , and by decomposing the vectors and into their normal component and projection into the tangent plane and using the impermeability condition on Γ wall we conclude from (11) that This structure is consistent with the boundary conditions (6) and (7) to be considered in this study. Taking into account the incompressibility condition (1), (12) can be expressed as (see Braack and Mucha, 2014 ) As the fluid is incompressible there is an indeterminate part to the stress that is spherical. Setting − ∶= 1 3 ( 11 + 22 + 33 ) , i.e. the mean normal stress, we have = − + , where ∶= − 1 3 ( 11 + 22 + 33 ) . Neglecting, for simplicity, the last (non-negative) term on the left-hand side of (13) and applying the integration w.r.t. time over (0 , ) , ∈ (0 , ] , we get The above equation states that the total kinetic energy at time and the sum of the dissipated energy due to internal friction in the fluid considered over (0 , ) × Ω and the dissipated energy due to friction caused by interaction of the fluid with the inner part of the boundary is controlled (bounded) by the initial kinetic energy. In order to be consistent with the second law of thermodynamics, the total dissipated energy has to be non-negative. This is for example automatically fulfilled by requiring that where * and ̃ * are positive constants representing the dynamic viscosity and the friction coefficient, respectively. We note that Eqs. (15) and (16) are constitutive relations and are linear. In order to incorporate a parameter that measures the various degrees of slippage and includes also no-slip ( = ) and (perfect) slip ( ( ) = ) boundary conditions we introduce the equation where 0 ≤ ≤ 1 instead of (16) . To conclude, upon inserting (15) and (17) into (14) , we obtain the following energy inequality for ∈ (0 , 1) : For = 0 or = 1 the energy inequality takes a simpler form, namely Next, we consider the complete formulation (1) - (7) , i.e. with general inflow velocity in . We assume the existence of an extension * in of in , i.e. we require that there exists a continuous function To obtain the energy equality, we modify the above procedure and form the scalar product of (2) with − * in . Doing so, we obtain Then we integrate this equation over Ω and use Gauss's theorem. Using also the properties of * in given in (20) , we obtain Upon inserting the constitutive equations (15) and (17) (noticing that the boundary integral on Γ wall vanishes for = 0 and = 1 prior to inserting) and neglecting the last (non-negative) term on the right-hand side of (22) we arrive at where the last term has been added to the left-hand side. Next, integrating (23) w.r.t. time and using the standard notation for the norms in the Lebesgue spaces ( (Ω) , ‖ ⋅ ‖ ) for 1 ≤ ≤ ∞, and using also the standard Cauchy-Schwartz, Young and Hölder inequalities, we obtain where the constant (20) . Finally, moving the second term from the right-hand side to the left and applying the Gronwall lemma we conclude that * sup

Definition of weak solution to the problem (1) -(7)
Let Motivated by the above energy estimates, by the concept of Leray-Hopf weak solution for the Navier-Stokes equations with no-slip boundary conditions, and by the results concerning the -maximal regularity for the evolutionary Stokes system, we say that a couple ( , ) is a weak solution to (1) - (7) provided that Alternatively, instead of (27) , one can consider valid for all ∈ 1 , 2 div , bc and a.a. ∈ (0 , ) .

The Nitsche method
In this section we modify the weak formulation (27) so that the incompressibility constraint (1) and the impermeability condition (7) 1 on the wall are treated in a manner suitable for numerical implementation. Thus, the function space 1 , 2 div , bc introduced in the previous section is modified in such a way that it takes into account only the inflow boundary condition, i.e. instead of using 1 , 2 div , bc the velocity is now sought in the space such that Then we proceed as follows. First, the incompressibility constraint (1) is multiplied by the test function ∈ 2 (Ω) and integrated over Ω. Then, the balance of linear momentum (2) is multiplied by the test function ∈ integrated over Ω followed by the integration by parts with respect to the term involving − div . Next, the impermeability condition written with respect to the form − = on Γ wall is multiplied by the test function ∈ 2 (Γ wall ) 3 and integrated over Γ wall . Finally, we sum the established identities and obtain where we use the notation ( , ) ∶= − + 2 ( ) and where we also used the observation that The formulation (30) requires special test functions defined on the boundary Γ wall which one needs to enlarge the problem. It is tempting to replace this new test function by a combination of and , which is possible as acts on which is already the unknown.
Inspired by the structure of the last two terms in the middle line of (30) we use the test function = ( ( , ) ) and obtain the original Nitsche method. Nitsche (1971) showed that such a formulation is not stable, but when one adds the stabilization term the weak form becomes stable for sufficiently large parameter > 0 . The local size of mesh edge is denoted by ℎ . Since the extra term − ∫ Γ wall ⋅ ( ( , ) ) d has the same sign as the term − ∫ Γ wall ( ( , ) ) ⋅ d , this method is called the symmetric Nitsche method . The weak formulation to our problem using the symmetric Nitsche method then reads: for all ( , ) ∈ × 2 (Ω) and a.a. > 0 .
The main disadvantage of adding stabilization (32) is that the parameter is problem-dependent and it has to be manually adjusted. Later, Burman (2012) showed that the original Nitsche method can be improved (in the sense of dropping ) by changing the sign in front of the boundary integral which corresponds to taking = −( ( , ) ) in (31) . This method is called the non-symmetric Nitsche method and its weak formulation reads: for all ( , ) ∈ × 2 (Ω) and a.a. > 0 .
In our computational results presented below we use both formulations (33) and (35) . We conclude this section with a brief literature overview regarding the Nitsche method. A finite element approximation of incompressible Navier-Stokes equations with slip on the boundary condition is investigated in Verfürth (1986Verfürth ( , 1991 . These papers provide the relation of the Stokes flow with slip on the boundary to the well known Babu š ka paradox for boundary supported shell plate, (see. Babu š ka, 1959;Rieder, 1974 ), where a discrete problem with piece-wise linear boundary approximation converges to a solution that is different from the continuous one. As observed in Verfürth (1986) the shell problem studied by Babu š ka is in fact equivalent to a stream function formulation of Stokes flow with slip on the boundary. To work around this problem they show that mixed formulation with the Lagrange multiplier method to impose the slip boundary condition for Navier-Stokes equations is stable and converges to the desired solution. In Dione et al. (2013) ; Urquiza et al. (2014) , a comparison of penalty method, Lagrange multiplier method and several variants of the Nitsche method to impose slip boundary condition on curved boundary in 2D is presented with observation that some variants do not exhibit convergence predicted by the theory. Juntunen and Stenberg (2009) derived stability and error estimates for the symmetric Nitsche method applied to a Poisson equation  with Dirichlet boundary data. In Layton (1999) it is proposed that weak imposition of no-slip boundary condition is superior to strongly enforced conditions and in Becker et al. (2015) the Nitsche method is used to apply general boundary conditions for Euler and Navier-Stokes equations. Finally, Mekhlouf et al. (2017) presents comparison of the symmetric and non-symmetric Nitsche methods for enforcing Dirichlet conditions in a weak sense for several simple 2D flow problem.

The finite element discretization
To solve numerically the variants (33) and (35) we discretize the systems by Taylor-Hood finite elements 2 ∕ 1 in space ( Boffi et al., 2013;Wieners, 2003 ). For the finite element discretization the FEniCS library ( Alnaes et al., 2015 ) is used. Since the given parameters in Table 1 require a sufficiently fine mesh to obtain stable solution, we add the standard SUPG stabilization ( Boffi et al., 2013 ) to obtain approximate solution on coarse meshes. For the time discretization we use simple implicit Euler method and perform the iteration process until a steady state solution is reached, see TSPSEUDO in PETSc ( Abhyankar et al., 2018;Balay et al., 2020 ).

The importance of definition of discrete boundary normal vector
It is clear that the discrete implementation of such a method can depend on the choice of boundary discretization and the definition of normal and tangent vectors to the boundary can influence the computed quantities. See Sime and Wilson (2020) , Engelman et al. (1982) , Stokes and Carey (2011) for similar discussion.
To give an interested reader a warning regarding this issue we give one striking example below. Before, we need to fix some notation. We approximate our smooth domain Ω by a piece-wise linear (planar) polyhedral domain with boundary Ω ℎ . The corresponding discrete normal vector ℎ is denoted as facet normal, see Fig. 2 , and its piece-wise constant function on the surface mesh Ω ℎ . It is the normal provided by the used FEniCS numerical library as ℎ = FacetNormal(mesh) . We can also introduce a vertex normal by taking the normal vector in vertices, i.e. ℎ . This can be constructed in several ways, see for example ( Dione et al., 2013 ). We chose to define a vertex normal as a 2 -projection of ℎ to space of continuous piece-wise linear functions on Ω ℎ , see again Fig. 2 . This is realized by computing the 2 projection as Finally, in the case of the straight pipe we can use the radial vector scaled to unit length as an analytic normal to the surface, this variant is denoted by ℎ , that means There are some additional possibilities with regard to defining the boundary normal, for example using the distance function ( ) = distance ( , Ω ℎ ) to the boundary Ω ℎ , we have ‖ ∇ ‖ = 1 . This might be a priori known or can be obtained by solving the Eikonal equation, or its regularized version. Then the normal vector could be obtained as ℎ = ∇ , however, in our simple situation, this performs almost identically to the vertex normal ℎ and this is why we do not report any results for this approach here. Now we return to the promised example. Let us consider a simple calculation of Poiseuille flow with full slip boundary condition. The analytical solution for this problem is trivial: a constant pressure and a constant velocity field = (0 , 0 , ) . Fig. 3 shows the numerical results for different choices of normal vector. The left panel shows the numerical result if we use the radial vector ℎ as the normal vector. Then we recover the continuous solution. However if we use the normal ℎ corresponding to the facets of the computational domain Ω ℎ , shown in the middle panel, the approximation of the velocity field becomes very bad ( ≈ 40% error). On the other hand taking the vertex based normal ℎ reduces the error to about 4% , which is shown in the right panel. Yet this is a very simple problem.

Description of benchmark
The problem presented in Section 2.1 will be now studied in the regime of steady flows in the cylindrical domain shown in Fig. 1 . It means that for given in ∶ Γ in → 3 and ∈ we study the following problem: to find ( , ) ∶ Ω → 3 × satisfying In the rest of this section, we proceed as follows. First, in Section 5.1 , starting from the assumption that the cylinder is infinite, the flow is steady and takes the form of a simple shear flow, we derive the analytical solution for steady flows of the Navier-Stokes fluid satisfying  (0 , 0 , ) . The 3 discrete solutions correspond to the symmetric Nitsche method with 3 different normal vectors -analytic, facet and vertex normals ( ℎ , ℎ , ℎ ). The color shows the relative error between the discrete solution and the analytical one, i.e. the boundary condition (40) on Γ wall . This solution will serve both as the inflow boundary condition in in the problem (1) -(7) but also as a benchmark used for testing the numerical methods and their implementations. Then, we will define the dissipation, vorticity and pressure drop in Section 5.2 and compute their values for the analytical solution found in Section 5.1 . The description of the benchmark and the values of parameters are provided in Section 5.3 .

Steady flow in the cylinder
For a steady flow in a cylinder with radius located along the − axis, see Fig. 4 , one can look for the velocity in the form of simple shear flow, i.e., Then, the incompressibility constraint (36) is fulfilled and the system of Eqs. (37) simplifies to Here, has the meaning of the pressure gradient and can be specified as where and are defined as the pressure average over two cross sections Γ in and Γ out , i.e., Note that from (42) it follows that = − 1 2 + and = 1 2 + for the cross-sections Γ in and Γ out located at = ∓ ∕2 . Denoting further the mean inflow velocity as it is well known that the analytical solutions for the limit cases of noslip ( = 1 ) and of full-slip ( = 0 ) take the parabolic or constant velocity profiles, i.e., ( , ) = 2 ( 2 − ( 2 + 2 )) 2 for = 1 , Therefore, for ∈ (0 , 1) , we can use the ansatz 2 with general parameters 1 , 2 that may depend on . With this ansatz, we obtain From (46) and (40) 2 we get the following restrictions: This leads to Convergence of relative errors with mesh refinement -Naviers slip = 0 . 5 , non-symmetric Nitsche, analytic normal. DOFs: number of degrees of freedom; EOC: estimated order of convergence; Ξ Ω : bulk dissipation; Ξ Γ : wall dissipation; drop : pressure difference between inlet and outlet.  Table 3 Convergence of relative errors with mesh refinement -Navier's slip = 0 . 5 , non-symmetric Nitsche, facet normal.  Table 4 Convergence of relative errors with mesh refinement -Navier's slip = 0 . 5 , non-symmetric Nitsche, vertex normal.  Convergence of relative errors with mesh refinement -Navier's slip = 0 . 5 , symmetric Nitsche, analytic normal.  Table 6 Convergence of relative errors with mesh refinement -Navier's slip = 0 . 5 , symmetric Nitsche, facet normal.  Table 7 Convergence of relative errors with mesh refinement -Navier's slip = 0 . 5 , symmetric Nitsche, vertex normal. Furthermore, and we get the following explicit formula for the dissipation and pressure drop, namely

Benchmark parameters
The geometrical and flow parameters were chosen to be as close as possible to those for blood flow in the aortic root or a segment of thoracic aorta (motivated by the assessment of clinical significance of a stenosed aortic valve). The geometry considered is that of a 12 mm radius (24 mm diameter) cylinder. This is consistent with typical aortic dimensions. The geometrical data together with the density, the dynamic viscosity and the inflow velocity were taken from Š vihlová et al. (2017) . For geometrical data see also Madukauwa-David et al. (2018) .
There are two slip parameters prescribed on the wall, * and , and they are related to each other, see (40) . We use the parameter * = * * as a characteristic length of the geometry, specified as radius of the inlet . The values of the constants and their units are shown in Table 1 (7) .

Results
Next we can look at the numerical results for the model described in Section 2.1 as compared to the analytical solutions derived in Section 5.2 . This gives us the opportunity to asses several aspects of the discretization method. First of all, the mesh convergence can be compared for the symmetric Nitsche formulation with different penalization parameters and the non-symmetric Nitsche formulation. Then, with a selected formulation we compare the quantities of interest for the full scale of boundary conditions from full slip ( = 0 . 0 ) to no slip ( = 1 . 0 ).

Symmetric vs. non-symmetric Nitsche formulation
Tables 2 -7 summarize the behavior of the discrete solution with respect to the mesh refinement. The computations were carried out on the sequence of 4 regularly refined meshes. Fig. 5 collects the convergence behavior for relative errors of velocity and pressure in 2 norms . Each graph compares the behavior of the error for the 3 variants of discrete normal vectors. It is clear that the analytic normal gives the most accurate solution, however this is available only in special cases and can not be used in a general geometry. Behavior for the two other normals now depends on the variant of the Nitsche method used. The non-symmetric Nitsche method gives almost the same results for both normals -facet based and vertex based. In the case of symmetric Nitsche method the vertex based normal gives errors comparable with the non-symmetric method but the facet based normal results in much larger errors and if the penalty parameter is not chosen correctly it seems to converge at a lower rate of convergence. Fig. 6 portrays the same comparison for the dissipation. Here the observation is similar -namely that the combination of the symmetric Nitsche method with the facet based normal leads to much larger errors and slower convergence compared to all the other combinations. We note that the selection of the value for the penalty parameter was given by the fact that for smaller values the convergence is not achieved, probably due to lesser stability, and for higher penalty parameter the resulting accuracy degrades.
The final observation is that the results for the non-symmetric variant of the Nitsche method are not sensitive to the choice of the normal vector, the only significant difference is in the wall impermeability satisfaction, see the column "flux " in Tables 4 and 3 , where the facet based normal results in a significantly better enforcement of the impermeability boundary condition on Γ wall . On the other hand, the results for the symmetric variant of the Nitsche method exhibit significant dependence on the choice of normal and choice of the penalty parameter.
To summarize the observations: 1. The non-symmetric Nitsche method (with = 0 ) is comparable to the symmetric Nitsche method with carefully selected penalty parameter. 2. The influence of normal selection on the results is much smaller for the non-symmetric Nitsche method than for the symmetric one. 3. Facet normal ℎ leads to better conservation properties (as observed in Engelman et al., 1982 ). However, the derived quantities (dissipation, pressure) are more accurate with vertex normal ℎ .

The results for hemodynamic quantities
Based on previous observations we present the results only for the non-symmetric Nitsche method since this method proved to be more robust and more precise. We focus on the hemodynamic quantities introduced in Section 5.2 .
First, we compare the values for these parameters on two different meshes for 11 slip parameters , = 0 , 0 . 1 , 0 . 2 , … , 1 . Coarser mesh denoted by c0 has 181,715 tetrahedra and mean edge length of 0.82 mm resulting in 925,365 degrees of freedom after finite element discretization. Finer mesh denoted by c1 has 245,571 tetrahedra, 1,206,092 degrees of freedom and mean edge length 0.72 mm. The relative errors between the computed and analytical solutions are summarized in Table 8 .
Graphs of the dissipations Ξ Ω , Ξ Γ , pressure drop flux and vorticity computed on finer mesh depending on the slip parameter are shown in Fig. 7 . The computed cases for 11 different are displayed as red dots. The blue line represents the analytical solution. We observe very good fit for all parameters as the relative error is less than 0.4 % , see Table 8 . For the dissipation quantities, the errors are less than 0.07 % .
Finally, the graphs of the selected hemodynamic quantities (namely dissipation, pressure difference, vorticity and wall shear stress) depending on the parameter are shown for different * in Fig. 8 . Recall that both * and are connected to each other through (7) . In the literature, they are sometimes replaced by one parameter = − * (1− ) , where = 0 for = 1 and → ∞ for = 0 . Fig. 8 demonstrates the importance of the proper choice of the * .

Conclusion
The main purpose of this study was to provide a benchmark for threedimensional flows of the Navier-Stokes fluid subject to slip on a part of the boundary that can be used for any computational code developed for solving such problems to verify its validity and correctness in a relatively simple test. The benchmark is based on analytical solutions for flows in a tube (Poiseuille flow) subject to slip boundary conditions on the wall. The benchmark is used to check not only velocity/pressure solution but also global quantities such as dissipation, pressure drop, wall shear stress and vorticity.
We have also provided a novel viewpoint on the Nitsche method used to numerically incorporate the impermeability condition into the weak formulation. We presented two variants of the method: the symmetric and non-symmetric one. It turns out that the non-symmetric version performs in a better way without the need for the stabilization term.
We also wish to point out how important the choice of the method used for approximating the normal vector is with regard to obtaining the solution. In some cases, even for a simple problem, the simplest choice of the facet normal leads to very high approximation error.
In parallel with the formulation of the benchmark, we have developed a robust solver that proved capable of solving the benchmark prob-lem on hand, and can be employed for solving more complex cardiovascular problems while using more realistic geometry and velocity data.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.