An accurate finite element method for the numerical solution of isothermal and incompressible flow of viscous fluid, part I: computational analysis

Despite the numerical challenges, finite element method is used for the simulation of viscous fluid flow. A consensus on the cause of numerical problems has been reached; however, general algorithms---allowing a robust and accurate simulation for any process---are still missing. Either a very high computational cost is necessary or some limiting procedure is used by adding artificial dissipation to the system such that the accuracy is in danger. As a means of achieving a reliable and accurate method independent on the underlying process, we attack the problem by starting a discussion about the admissible governing equations leading to the necessary weak form without use of any stabilization or splitting schemes. We implement and solve the weak form by using open-source packages and present the accuracy by comparing the simulation to the closed-form solutions in the case of a linear viscous fluid. The computational method is accurate and robust for simulating isothermal and incompressible flows.


Introduction
Isothermal flow of viscous fluid is governed by balance equations of mass and linear momentum: in Cartesian coordinates, where ρ denotes the mass density, v i the velocity, σ ij the non-convective flux term (Cauchy's stress), g i the specific supply (gravitational forces); here and henceforth we apply Einstein's summation convention to repeated indices. A linear relation for stress furnishes the governing equations with accurate detail for materials response in the case of Newtonian fluids such as water. This linear relation is the well-known Navier-Stokes equation: with the material constants λ, µ; and a new parameter p called hydrostatic pressure. Consider a control volume initially filled with homogeneous water, i.e., initial mass density is constant in space. Under conventional pressure differences, we assume that the flow be incompressible, i.e., the mass density remains constant in time.
Then, from the mass balance, we obtain which is used for computing the pressure p. For an incompressible flow, the mechanical pressure − 1 3 σ ii becomes identical to the hydrostatic pressure p such that we handle p as the pressure generated in a pump. Velocity and pressure fields have to satisfy Eqs. (1) 2 and (3).
In analytical mechanics, the aforementioned equations hold true locally (in every infinitesimal point in space). For a computation we discretize the space, for example with finite sized elements in the finite element method (FEM). Within each element, the analytical functions for the unknowns v i and p are represented by form (shape) functions with a local support, i.e., by means of a discrete element. The shape functions are not smooth, they belong to C n with a finite n. In other words, the unknowns are finitely differentiable and depending on the governing equations and constitutive relations-from the mathematical analysis in [6]-we know that the correct choice of the form functions for velocity and pressure is of paramount importance for a robust computation. This so-called Ladyzhenskaya-Babuska-Brezzi (inf-sup compatibility) condition (LBB condition) tells us how to adjust the shape functions of velocity and pressure.
Within a finite element, the governing equations are written as holding globally (over the domain of the element) instead of locally meaning in every location in the element. There exists a general assumption that we can use the same local governing equations holding globally in finite elements; however, this strategy leads to several numerical problems and to various proposals in [7,13,14,25], for a review of such suggestions see [19]. These so-called stabilization methods introduce a parameter depending on element-length. This parameter induces an artificial viscosity in the direction of the velocity change for suppressing and eliminating the numerical oscillations. Otherwise, a computation is supposed to be not possible since the numerical oscillations cause to divergence of any numerical solution schema. There are several successful implementations as in [16,22,15,9,12,18]. From a practical point of view, such a stabilization method is useful; however, it is challenging to obtain the parameters necessary for the stabilization since they are not measurable quantities. In other words, we need to fine-tune the numerical parameter for every different application and we cannot assure that the computed solution is approximating the reality accurately, see references in [8] about accuracy problems of different techniques. By using small elements-smaller than the characteristic length-scale, for example Kolmogorov scale-we can overcome these numerical problems. This direct numerical simulation (DNS) is often not feasible without access to super-computers. Since small elements allow us to use the local governing equations, actually the stabilization methods become necessary as a consequence of the assumption that the local governing equations be holding globally. Therefore, we actually need to find out the suitable equations holding globally. Based on this idea, a gedankenexperiment in [2] lead to the discrete-mechanical method including terms mimicking the stabilization terms. The main difference between stabilization methods and discrete-mechanical method is that the discrete-mechanical method generates a weak form for the computation without any process depending parameters.
In this work, we discuss a special yet general case, namely an isothermal and incompressible flow. For understanding the numerical problems, often, the pressure related numerical problems and velocity related numerical problems are studied separately. We use one balance equation for each of them and expect that the coupling due to the constitutive equation for stress is capable of relating each other with sufficient accuracy. However, we know from the LBB condition that their relation is more complicated. Hence, we search for the governing equations supplying necessary information about the pressure and velocity fields globally instead of locally. Numerical problems are surpassed by incorporating the balance of angular momentum delivering the necessary smoothness for the pressure. Conventionally, the balance of angular momentum is neglected since it is fulfilled locally. We discuss this issue and explain in detail how to implement this equation and generate the weak form. Moreover, we use open-source packages and solve some academic examples in order to present the accuracy and robustness of the method. All codes are made public on the web site in [3] to be used under the GNU Public license as in [10] for promoting an efficient scientific exchange as well as further studies.

Variational formulation
Consider the following general balance equation: where the rate of the variable ψ is balanced with the supply term z acting volumetrically and with the flux terms-convective f and non-convective φ-applying on the surface ∂Ω of a domain (control volume) Ω. By using Table 1 we can obtain the balance equations of mass, linear momentum, and angular momentum. The domain may have its own velocity x • i = w i independent on the velocity of the fluid particle, v i . For a discussion of the balance equations in a control volume with the domain velocity, we refer to [20,21]. This domain velocity can be chosen arbitrarily without affecting the underlying physics. For the sake of clarity, consider an Table 1: Volume densities, their supply and flux terms in the balance equations.
where we record the motion of fluid by means of a high-speed camera. We may move the camera arbitrarily, the motion of the fluid fails to alter. Normally we fix the camera to the laboratory frame such that the domain has no velocity, w i = 0. We will use this assumption of a fixed domain having two consequences. First, the time rate reads Second, the rate of the infinitesimal volume element vanishes The domain velocity becomes important in computer methods for fluid structure interaction.
We start with the balance of mass for a fixed domain and apply the Gauss-Ostrogradskiy theorem We assume that initially the fluid rests, v i (x, t = 0) = 0, and it is a homogeneous material, ρ(x, t = 0) = const. Moreover, we assume that the flow is incompressible, i.e., the mass density remains constant in time, ∂ρ/∂t = 0. Then the mass density is constant in space and time throughout the simulation leading to We multiply the latter with an arbitrary test function with the same rank of the integrand, i.e., a scalar. According to the Galerkin approach, we will choose the test functions from the same space as the unknowns v i and p. Therefore, it is natural to use δp and δv i as possible test functions. We emphasize that the choice of the scalar test function is critical and we suggest to use which is also the common procedure in the literature. Actually by multiplying with the test function, we started to utilize the discrete representation of the continuous velocity field. We omit a clear distinction in the notation since we never use continuous and discrete functions together in one equation.
For the balance of linear momentum, the procedure is the same, we acquire for a fixed domain and incompressible flow The integrand is a vector, hence, we use the test function δv i as being the case in the literature, We refrain ourselves from inserting the mass balance, ∂v i /∂x i = 0, into the latter formulation. Since this condition is tested by δp, we cannot expect that it is fulfilled for the velocity distribution such that we choose to enforce it by leaving the formulation as it is.
The same formalism is applied for the balance of angular momentum-for an incompressible flow ρ = const| x,t within a fixed domain where ∂x i /∂x j = δ ij and the spin density s i , its flux term (couple stress) m ij , and its supply i vanish for a non-polar medium like water (furnishing a symmetric stress). Then the angular momentum is identical to the moment of (linear) momentum, hence, in analytical mechanics, we can neglect it. However, we observe this term as an important aid for resolving numerical challenges such that we want to enforce it by multiplying the balance of angular momentum by a vector test function. If we multiply the balance of angular momentum for non-polar medium by δv l , then the outcome is identical to the variational form in Eq. (11). However, if we multiply by ∂δp/∂x l we obtain a restriction about the gradient of pressure, which is indeed the missing condition for the necessary numerical smoothness for pressure. We may see this condition related to the LBB condition; however, we enforce it by using an additional integral form instead of changing the order of shape functions. The variational formulation for the balance of angular momentum reads which is one of the key contribution of this work providing stability and robustness.
We solve the transient integral forms in discrete time slices with a time step ∆t, this discretization in time can be established by using Euler backwards method where v 0 i indicates the value from the last time step. This implicit method is stable (for real valued problems) and easy to implement.

Generating the weak form
We have to bring integral forms in Eqs. (9), (11), (14) into the same unit such that they can be summed up. We divide Eq. (9) by the mass density ρ and bring it to the unit of power. Equation (11) is already in the unit of power. We divide Eq. (14) by the mass density and multiply it by the time step ∆t and bring it also to the unit of power.
We rewrite the constitutive equation, σ ij = −pδ ij + τ ij , by combining all terms with d ij into τ ij . The symmetric part of the velocity gradient, d ij , thus, the term τ ij , include already velocity's first derivative in space. Therefore, we need to have a representation of the velocity function, which is at least C 1 continuous. In the integral forms we observe another space derivative to τ ij leading to class C 2 for velocity approximation. This condition can be weakened by integrating by parts. We emphasize that the integration by parts is applied only on the terms already including a differentiation in order to switch one differentiation to the test function.
This strategy is not conventional. Often, integration by parts is used to all flux terms. After bringing to the same unit and integrating by parts where necessary, we obtain by using the usual comma notation for space derivative () ,i = ∂()/∂x i . The weak form for one finite element reads zero by inserting the correct pressure and velocity distribution. A control volume is decomposed into several elements such that the assembly over all elements, Form = e F Ω , is the weak form for the whole control volume. The sum over elements generate on the element boundaries the following term with jump brackets (·) indicating the difference between the values of the quantity computed in adjacent elements. We use continuous form elements for pressure and velocity such that n i p = 0, moreover, we enforce n i σ ij = 0 relying on the balance of linear momentum on singular surfaces. Therefore, the integral term along the element boundaries vanish within the control volume. On the boundaries of the control volume either velocity or pressure is given. For the parts, where velocity is given, we choose δv i = 0 such that the boundary term vanishes. For the parts, where the water enters the domain with a prescribed pressure p against the plane outward, t i = −pn i , we obtain n i τ ij = n i (σ ij + pδ ij ) = t j + n j p = 0. Therefore, all boundary terms vanish in the weak form for an incompressible flow.

Algorithm and computation
where [H n ] 4 is a 4-dimensional Hilbert space of class C n as defined in [11] with additional differentiability properties such that it is called a Sobolev space. The test functions, δP = {δp, δv 1 , δv 2 , δv 3 }, stem from the same spaceV = δP ∈ [H n (Ω)] 4 : δP | ∂Ω = given , which is the Galerkin approach. We use n = 1 for all simulations, in other words, we use the same linear form functions for pressure and for velocity. The spurious oscillations do not occur owing to the additional governing equation restricting the pressure gradient as well as the careful choice of terms to integrate by parts.
The weak form is nonlinear and coupled such that we need to linearize and solve it monolithically. For the linearization we follow the ideas in [17, Part I, Sect. 2.2.3] and perform an abstract linearization using Newton's method at the partial differential level. The functional Form= F (P , δP ) is an integral of the function depending on the primitive variables P and their variations (test functions) δP . We know the correct values of P at t 0 . The weak form is initially zero-we obtained it by subtracting left-hand sides from right-hand sides in the balance equations. We search for the next time step, t = ∆t + t 0 , by using the known P . This algorithm holds for every time steps, since we compute subsequently in time. We may describe the algorithm: given: P (t) for x , find: P (t + ∆t) at x , satisfying: F (P (t + ∆t), δP ) = 0 .
Now, by rewriting the unknowns P (t + ∆t) in terms of the known values P (t + ∆t) = P (t) + ∆P (t) , (22) we redefine the objective to searching for ∆P (t) instead of P (t + ∆t). If ∆t is chosen sufficiently small, then the solution is near to the known solution such that ∆P (t) is small. This condition leads to a Taylor expansion around the known values, P (t), up to the (polynomial) order one where we omit the time argument for the sake of clarity in notation. The expansion is linear in ∆P , hence we need to construct a linear in ∆P differentiation operator, ∇ P , which is established by the so-called Gateaux derivative: where is an arbitrary parameter. Since we first differentiate in and then set the parameter equal to zero, only terms of order one in ∆P remain in the solution. By introducing the so-called Jacobian: we rewrite the algorithm, given: P for x , find: ∆P at x , The last line is a linear function in ∆P such that we can solve the equation by obtaining ∆P and update the solution in an iterative manner, where ":=" is an assign operator in computational algebra. Here is the ultimate algorithm: while |∆P | > TOL.
solve ∆P , where F (P , δP ) + J (P , δP ) · ∆P = 0 P := P + ∆P (28) The term J · ∆P is computed automatically by means of symbolic differentiation implemented under the name SyFi within the FEniCS project, see [4], [5]. This automatic linearization procedure allows us to use any nonlinear constitutive equation in the code. Herein we use a linear constitutive equation in order to achieve a comparison with closed-form solutions. The geometry is constructed in Salome 1 by using NetGen algorithms 2 for the triangulation. Then the mesh is transformed as explained in [1, Appendix A.3] and implemented in a Python code using packages developed by the FEniCS project, 3 which is wrapped in C++ and solved in a Linux machine running Ubuntu. 4

Comparative analysis
The suggested weak form is compared with semi-analytical closed-form solutions for simple geometries, we call them analytical solutions. They are well-known and can be found in different textbooks, for example see [23]. This analysis is of importance to show the accuracy of the suggested weak form. A similar precision can be obtained by using DNS with an extreme computational cost. Herein we use a modest workstation equipped with Intel Xeon processor (8 cores) running at 3.4 GHz. Then we present a computational challenging problem in order to test the robustness of the proposed strategy. For all problems we use a direct solver (mumps) solving in parallel (up to 8 cores).

Steady state Hagen-Poiseuille flow
Consider a laminar flow in an infinite pipe as a result of the given pressure difference. This configuration is called the Hagen-Poiseuille flow and it has a steady-state solution obtained in cylindrical coordinates, r, θ, z, under the assumption that the flow is only along the pipe. The pipe is oriented along z, which is set as the axis of the pipe. No-slip condition, v i = 0, is applied on the outer walls, r = a. The steady-state solution reads for an incompressible flow of a linear viscous fluid like water of viscosity µ. We use this solution for an analysis of the convergence and accuracy. A finite sized pipe of length is constructed and on the inlet and outlet the pressure is given as Dirichlet boundary conditions, p(z = 0) = p in and p(z = ) = p out . The flow is driven by the pressure difference such that the gravity is neglected. Moreover, we are interested in the steady-state where the inertial term vanishes. Furthermore, in order to mimic the infinite pipe, we set the radial and circumferential velocities zero on the inlet and outlet. FEM computation is realized in Cartesian coordinates, so we basically allow v z on inflow plane x = 0 and outflow plane x = by setting v x = v y = 0 via Dirichlet conditions.
The pressure distribution is expected to be linear along z, hence, for the analytic solution dp/ dz = (p out − p in )/ . For a better comparison we use the diameter D as the characteristic length and half of the maximum velocity v M z = v z (r = 0), as the characteristic velocity for calculating the Reynolds number: For a small pipe of an inch long and a quarter inch wide, we construct the mesh by using a global element length h. We use SI units such that = 25.4 mm and D = 6.35 mm as well water as the fluid with ρ = 998.2 · 10 −6 g/mm 3 , µ = 1001.6 · 10 −6 Pa s , λ = 0.6 Pa s .
Especially the choice of λ is of importance since this parameter is not measured directly. In the constitutive equation in Eq. (2), the term λd kk vanishes for the correct velocity solution. Therefore, for the numerical sense, λ has to be great enough that d kk is enforced to vanish. Hence, we choose λ multiple times greater than µ, in reality (for compressible flows) they are independent parameters. Interestingly, we have observed that choosing λ greater than suggested value slows down the convergence. In other words, by doubling λ the same accuracy of the solution is obtained with more degrees of freedom (DOF). Its value does not change the convergence behavior, as long as it is great enough. In order to determine the value of λ, we simply decreased until the maximum Re was achieved. Less than the used value leads to numerical problems in the Newton-Raphson iterations.
By using a standard convergence analysis, we compare three different meshes starting with a global edge length of tetrahedrons, h = 0.6 mm, then reducing it by half. Every time a new mesh is generated such that the number of nodes fail to increase exactly by 2 3 times in 3D. We use an unstructured mesh and the mesh quality is nearly identical because of using the same algorithm on the relatively simple geometry. The expected monotonic convergence has been attained and can be seen in Fig. 1 for 2 different Reynold numbers. Two important facts need to be underlined. First, the parabolic distribution of the velocity along the diameter is achieved even with a coarse mesh. Second, the relative error is 1.1% at r = 0 for a low Reynolds number and 4.1% at r = 0 for a high Reynolds number in the case of a laminar flow achieved with 25 min of computational time using a direct solver.

Starting Hagen-Poiseuille flow
In order to test the accuracy in the transient simulation, we use the same configuration and solve it transiently in time. Since we have obtained the expected convergence in space discretization for the steady-state solution, we expect to have a monotonic convergence in time discretization, too. Therefore, we solve the same example with different time steps from t = 0 to t = 10 s with the initially applied pressure difference. For this case there is a closed-form solution under the same assumptions as before, with τ = tµ/(ρa 2 ) and Bessel functions (of the first kind) J 0 and J 1 with Λ n being the roots of J 0 . We compute the (semi-)analytical solution by using SciPy packages for Bessel functions as well as its roots and choose N = 50. By choosing the best mesh obtained from the convergence analysis in the steady-state case, namely the mesh with h = 0.15 mm, we compute the solution transiently in time by using different time steps. We present in Fig. 2 the maximum value v z (r = 0, t) over time for different time steps showing the expected convergence with decreasing time steps as well as the distribution at different time instants for the smallest time step. In addition to the convergence in time, the relative error remains unchanged over time. We conclude that  the suggested formalism is capable of simulating a simple, laminar, three-dimensional flow of a linear viscous fluid accurately under a monotonic loading. The velocity and pressure distributions show no artifacts or mesh dependency as presented in Fig. 3. We emphasize that the pressure difference is applied instantaneously, which is numerically challenging. For transient loading scenarios, there are various assumptions used for obtaining a closed-form solution such that we omit to examine further. Based on the presented examples, we conclude that the approach delivers an accurate solution.

Flow around a cylinder
As suggested in [24] we implement a three-dimensional pulsating flow over a cylinder with a square cross section (3D-3Q). This benchmark problem has no analytical solution, but it is seen as a numerically highly challenging problem used for testing the robustness of the developed solution strategy. We emphasize that the proposed approach makes a solution possible even with a coarse mesh, thus, we choose a relatively coarse mesh with less than 70 000 degrees of freedom-the mesh is shown in Fig. 4. The geometry is a box of 2500 mm x-length and H × H cross section on yz-plane. A through hole of 100 mm diameter is placed asymmetrically, its center is at y = 200 mm and the y-length is H = 410 mm. This configuration leads to an asymmetric flow around the cylinder. We use again water as the fluid with the aforementioned material parameters. The inflow on x = 0 is given as follows with the mean velocity of v m = 2250 mm/s. This parabolic function reads v m at the middle of the inflow plane y = z = H/2 at t = 4 s. On the outflow plane x = 2500 mm we set the pressure as the reference pressure. On the walls and around the cylinder, no-slip condition, v i = 0, is implemented as Dirichlet boundaries. A variable time stepping and relaxation parameter scheme is implemented. Every time step lasts between 5-10 s depending on the relaxation parameter leading to 2-4 Newton-Raphson iterations. Since we use a coarse mesh, we fail to capture all details of the flow; however, we visualize in Fig. 5 the velocity field at t = 4 s, when the maximum velocity on the inflow plane has reached 2250 mm/s leading to a Reynolds number of Re = 100 with the diameter of the cylinder used as the characteristic length, and also at t = 8 s, when the inflow diminishes by having decreased to zero leading to eddies.

Conclusion
A new method is proposed to compute fluid dynamics of isothermal and incompressible flows by means of FEM. We have carefully investigated the convergence and the accuracy of the proposed approach by using closed-form solutions. Convergence in space and time is difficult to achieve in FEM for flow problems. The accuracy is only possible by using extremely high degrees of freedom. With the proposed method, we have tackled both of this issues by attaining a monotonic convergence in space and time as well as accurate solutions with a coarse mesh. Therefore, at least for laminar flows, we expect the method to deliver accurate results. Moreover, the suggested implementation by using FEniCS has provided a robust code.