Numerical Investigation of Two-Dimensional Oldroyd-B Fluid Flows Over A Straight Rectangular Domain

In this paper, non-Newtonian viscoelastic Oldroyd-B fluid flows in two-dimensional rectangular domain is numerically investigated, where the flow between two rigid walls is driven by a pressure difference along -direction (horizontal). The numerical results of the nonlinear system of partial differential equations are obtained by decoupling the system into Navier-Stokes system and tensorial transport equation. Computational Fluid Dynamics (CFD) simulations are done by using the finite element method. The numerical simulations are presented in terms of the contours of velocity, pressure and extra stress tensor. The Hood-Taylor finite element method is used for the approximation of the velocity and the pressure while the discontinuous Galerkin method is used to approximate the stress tensor. All the meshes and simulations are carried out by the general finite element solver FreeFem++, which has been found as a potential tool to provide a reasonably good numerical simulations of complicated flow behavior.


Introduction
The aim of this work is to analyze and simulate the incompressible non-Newtonian Oldroyd-B fluid flows [1][2][3][4][5][6], in a straight rectangular domain in case of two dimension. The constitutive equations for Oldroyd-B fluid is highly nonlinear system of partial differential equations (PDE), and therefore it is decoupled into two auxiliary problems, namely, Navier-Stokes equations and the tensorial transport equation. CFD describes the Oldroyd-B fluid flow in terms of these two auxiliary problems. In the present paper, we consider the two-dimensional flow of an Oldroyd-B fluid over a straight rectangular domain. The auxiliary problems are studied separately in our previous works [7][8][9].
The Oldroyd-B model is a composed problem with three unknowns: the velocity components , the pressure , and the viscoelastic extra stress tensor . We use iterative scheme to solve this system. If is fixed, the model defines a Navier-Stokes systems in the variables and , which is analyzed using the Hood-Taylor finite element method for the approximation of the velocity and the pressure field , [8,9]. On the other hand, if (and ) is fixed, then the model is transformed into a transport equation in the variable . The approximation of is done by using discontinuous Galerkin finite element method [10][11][12], as discussed [7].
We briefly discussed the problem formulation, and mathematical and numerical analysis of Oldroyd-B model in the context of Hood-Taylor [13][14][15] and discontinuous Galerkin finite element method [10][11][12]. The variational formulation is discussed in terms of the auxiliary problems. To employ numerical simulations in Oldroyd-B model, Navier-Stokes problems is discretized using the Hood-Taylor finite elements and the discontinuous ℙ 1 elements (ℙ 1 ) are used to discretize the transport problem. The finite element solver FreeFem++ [16] is used to obtain the approximation of the composed problem. The numerical results of the coupled Oldroyd-B problem are obtained both computationally and graphically by the numerical simulations of these two auxiliary problems implemented on the general finite element solver FreeFem++ [16].

Nomenclature
Before discussing the mathematical analysis of Oldroyd-B fluid flows with boundary conditions and variational formulation, we introduced some notations of different function spaces in the following table, details of which can be found [17,18]

The Constitutive Equations for Oldroyd-B Model
Let is a bounded, open and connected Lipschitz domain of ℝ , = 2,3. Oldroyd-B model consists of the system of non-linear partial differential equations formed by the law of conservation of mass, the momentum equations and the form of transport equations, which can be written in the domain as [19]   Assuming the Oldroyd-B constitutive equations (1) can be written as The above set of equations describes the behavior of an incompressible viscoelastic fluid of Oldroyd-B type, in the open subset , where the fluid is homogeneous.
We observed that the conservation of momentum leads the symmetry properties of the tensor , i.e., = . The problem (2) is a mixed problem. The first two equations form a parabolic system for ( , ) which is in the form of Navier-stokes equation. The last equation has a hyperbolic characteristic which is in the form of Transport equation.

Non-dimensional governing equations
To obtain a system of dimensionless variables, we discuss some scaling properties of three equations of the system (1) to introduce Reynolds number and Weissenberg number that respectively measures the effect of viscosity and elasticity on the flow.
Let be the characteristic length, represents a characteristic velocity of the flow and = + be the viscosity coefficient. We transform the system (1) into dimensionless form by changing variables and by introducing the following dimensionless quantities: where the symbol is attached to dimensional parameters.
The dimensionless form of the system (1) can be written as where = 1 − 2 1 = + , is the retardation parameter (0 < < 1), and 2 ≥ 0 is the retardation time of fluid. Reynolds and Weissenberg numbers are the dimensionless numbers. Small values of means that the fluid is little elastic and small values of means that the fluid is very viscous. In case of stationary motion, the problem can be written as follows : Find the quantities , and defined in such that The above equation is composed of a Navier-Stokes like system for ( , ) and a transport equation for extra stress tensor .

Boundary conditions
The system of equations (3) has to be used with some boundary conditions. For a connected flow domain ⊂ 2 , the required boundary conditions are the following: (i) Dirichlet boundary conditions for the velocity on the boundary = on Ω with compatibility condition where is the unit outward normal vector to at the boundary . For homogeneous case, = 0.
(ii) For the stress, a condition on the upstream boundary section With the homogeneous Dirichlet boundary conditions, the Oldroyd-B fluid model problem is well-posed [20].
For non-dimensional case, the problem (6) can be read as the form of (5) as follows: Find the non-dimensional quantities, still denoted by , and , defined in defined in such that subject to the boundary conditions in (i) and (ii).

Algorithm to solve the discrete Oldroyd-B problem
To solve this elliptic-hyperbolic system, we applied the decoupled technique. The extra-stress tensor is computed separately from the kinematic equations. From a fixed value for the velocity (and pressure) the extra-stress tensor is evaluated solving the tensorial transport equation with fixed-point method. Then the velocity field and pressure are updated with the current extra-stress tensor whose components are treated as known body forces, solving the resulting Navier-Stokes equation by the Newton-Raphson method. This procedure is iterated.

Numerical Results and Discussions
This section is concerned with the application of the finite element method to obtain the numerical results for non-Newtonian viscoelastic Oldroyds-B fluid flows. By the implementation of the finite element method in our own script in FreeFem++, we obtained the numerical solutions of the Oldroyd-B problem.
Here, we considered the Oldroyd-B flow between two rigid walls where the flow is driven by a pressure difference along -direction. This flow is laminar and referred as Poiseuille flows.
The velocity is uniaxial and has a parabolic profile and we supposed that = 1.
The analytic solution for the kinematic is given by [22], The no-slip conditions on the two rigid walls are given by 1 = 0, 2 = 0. We assumeed that 2 = 0 everywhere at the inlet and 1 ( , ) = (1 − ) as the exact solution. At outlet we imposed 2 = 0. The condition for stress tensor on the upstream boundary section − agrees with the exact solution (14). The problem has been solved using four grids obtained by successive refinements dividing each triangle into four new triangles starting with a coarse mesh with 344 elements (Fig. 1). We consider the problem with = 1, = 1 and = 0.1. The Table 1 characterizes the mesh through the diameter , number of elements, degree of freedoms. In each case, we evaluated the error of fluid velocity in 1 -norm and the error of the components of the tensor in 2 -norm which are respectively defined by , , = 1, 2.
The results obtained for and , , = 1, 2 over the different meshes are present in the following table (Table 2). The good convergence of results for all kinematic can be confirmed by the slope value, which gives us the rate of convergence to the exact solution with respect to the corresponding norms. We used the least squares approximation to find the slope of the log-log plot of the error of the velocity and the tensor components. The following plots (Figs. [2][3][4][5] show the error curves:   The values for the rate of convergence (the slope) which we obtained guarantees the errors for all the variables approaches to zero as h tends to zero, which we expected theoretically. We are in conditions to affirm that the numerical solution converges to the exact solution. So, the algorithm is convergent, and we observed that the behavior of exact and numerical solutions is approximately same.
The approach solution obtained with 22654 elements is illustrated graphically (with color scales) in Figs. 6-11.      So, from all the above numerical and graphical results, we observed that the approximate solution converges to the exact solution with respect to the corresponding norms.

Conclusion
The main goal of this work was the mathematical and numerical study of the nonlinear system of partial differential equations that model the motion of incompressible non-Newtonian fluids of Oldroyd-B type, in dimension 2, in case of steady flow over a rectangular domain. The numerical simulations to the Oldroyd-B problem were obtained computationally by the implementation of the finite element methods (continuous for kinematic and discontinuous for the extra stress tensor) in a script of FreeFem++. This mixed problem of elliptic-hyperbolic type was decoupled into two auxiliary problems, namely, the Navier-Stokes system and the tensorial transport problem. Based on the numerical techniques described for both the auxiliary problems, the approximation of the solution of the Oldroyd-B problem was obtained.
Here, the numerical results have been obtained by considering the benchmark problem over a rectangular domain. We have observed the solutions both graphically and computationally, and also have observed the behavior of the solutions in the rectangular domain. We have found a good agreement between the approximate and exact solutions which confirms the convergency of the algorithm suggesting the capability of providing the better numerical approximations of Oldroyd-B model using FreeFem++.