Numerical simulation of Oldroyd-B fluid with application to hemodynamics

Oldroyd-B viscoelastic fluid is numerically simulated using the stabilized Galerkin least squares finite element method. The instabilities due to the connective nature of the Oldroyd-B model is treated using the discrete viscous elastic split stress method. The model is used to study the behavior of the flow of blood through an abdominal aortic segment. The results show that the viscoelastic nature of the blood should be considered specially at low shear rates.


Introduction
One of the most common problems which occur in the abdominal aortic artery is the abdominal aortic aneurysms (AAAs) where the artery has a balloon-shaped expansion, and hence, the increase in the lumen diameter reach to 50% of its normal diameter. 1 One of the causes of the AAAs is the loss of arterial wall integrity. 2 In men, it is found that AAA happens commonly in those above the fifth decade. 3 In this work, the geometry is considered as a twodimensional (2D) planar channel with non-deformable walls. Although this assumption is less realistic than other assumptions such as axisymmetric or threedimensional, this assumption is acceptable for geometries with spatial symmetry. There are many works for blood flow simulations that consider the assumption of planar channel for arteries. [4][5][6][7] Many studies 8-14 related to AAAs considered the blood as a purely viscous fluid with constant viscosity, that is, Newtonian fluid. In this work, the blood is modeled as a viscoelastic fluid where the Oldroyd-B constitutive model is used to represent the viscoelastic property of the blood.
The stabilized Galerkin least squares (GLS) method is used 15 since high Reynolds numbers are expected in blood flows. Due to lack of enough diffusion in the Oldroyd model, the discrete viscous elastic split stress (DEVSS) method 16 is used to suppress the associated instabilities. An in-house built FORTRAN code is built to solve the governing equations, and the results are visualized using the Techplot software. This article is organized as follows: the problem's geometry is described in ''Problem geometry'' section while the governing equations and stabilization techniques are described in ''Mathematical modeling,''''The incompressiblity constraint,'' and ''The DEVSS method'' sections. Numerical parameters of the problem are stated in ''Numerical simulation'' section. Finally, results are presented and discussed in ''Results and discussion'' section. Figure 1 describes the upper half of geometry used in the current work in which a 2D channel with single aneurysm is considered. Equation (1) is the mathematical description of the upper wall geometry 13

Problem geometry
in which D is the artery diameter at the inlet section and L T is the artery length with L 1 = 2:5D, L 2 = 5D, L T = 7:5D, and D 1 = 2D. In this study, the diameter of the artery is taken to be 8 mm. 13 Mathematical modeling

Balance laws
Considering non-polar incompressible, viscous fluids, the balance of linear momentum is and the continuity equation reads where r is the fluid density, V is the velocity vector, p is the pressure, m s is the solvent zero-shearrate viscosity, and t p is the polymeric contribution to the extra stress tensor. For Newtonian fluids, t p vanishes.

Constitutive laws
For viscoelastic fluids, the extra stress tensor is related to the velocity gradient via a differential equation. In the Oldroyd-B model, polymeric contribution is where m p is the polymer zero-shear-rate viscosity and l is the relaxation time. 17 Nondimensionalization scheme The following nondimensionalization scheme is used in which V m is the inlet average velocity and p exit is the exit pressure. The velocity vector is V = (u, v) and the polymeric stress tensor is in which S, Q, and T are the axial, shear, and normal stresses, respectively. The non-dimensional form of equations (2) and (4) are (the asterisk is dropped for clarity) where Re = (rV m D)=m is the Reynolds number, m = m s + m p is the total viscosity, We = l V m =D ð Þis the Weissenberg number, and b = m s =m is the viscosity ratio.

The incompressiblity constraint
To recover the link between the continuity and momentum equations, a modified form of the continuity equation is used using the pressure stabilized technique 18 in which, e is a controlling parameter used to stabilize the pressure. In this study, the value of stabilization parameter is chosen to be 0.0001, where the time step values range from 0.0001 to 0.00005 according to the value of Reynolds number. The main advantage with this technique is that equal order basis functions could be used for all primitive variables. 19 Figure 1. Geometry of the artery.

The DEVSS method
Since the additional polymeric stress induces numerical instabilities, 20 the DEVSS method is used in which the discrete form of the momentum equation is modified by adding extra diffusion. Introduce the Newtonian-like stress tensor D defined as D = m p (rV + (rV) T ) to the momentum equation which leads to a symmetric influence matrix. The momentum equation, equation (2), is reformulated as follows The extra variable D can be calculated from the following non-dimensional equation We can rewrite equation (9) in the non-dimensional form as follows This is the stabilized form of the balance of momentum equation.

Numerical simulation
The set of equations (7), (8), (10), and (11) are solved numerically using the finite elements method. The GLS technique is used in which the basis function is modified which leads to stable solutions. 10 The unsteady terms are discretized using the first-order explicit Euler scheme. The time step is chosen small enough to obtain accurate and stable solutions. Zero values for the flow variables are used as an initial guess.

Numerical parameters
The physiological condition 21-23 described in Table 1 is used in the current work.

Boundary conditions
The no-slip boundary condition is imposed over the walls. A fully developed flow is assumed at the artery inlet. At the artery exit, the normal component of the Cauchy stress tensor is set to zero, namely, the pressure, p, is set to zero as well as the normal stress T.
According to the theory of characteristics, two components of the stress have to be specified at the inlet. 17 Using the fully developed flow assumption at the inlet, the axial and shear stresses are 24,25

Results and discussion
Results are obtained for both cases, namely, the Newtonian case (t p = 0) and the viscoelastic case. In both cases, different values for the blood volumetric flow rate are considered as shown in Table 2. The corresponding Reynolds and Weissenberg numbers are described in Table 2 as well.

Grid-independent solution
Different mesh sizes are used to reach the meshindependent solution. The computational domain shown in Figure 2 is discretized using 4000 bilinear quadrilateral iso-parametric elements. Table 3 includes the number of elements for each grid. The wall shear stress (WSS) is calculated using all meshes for the volumetric flow rate Q = 3 cm 3 /s. Figure 3 shows that the last two solutions are very close to each other assuring that the finer the mesh, the better the solution.

Viscoelastic flow patterns
Parametric study is performed to investigate the effect of the blood volumetric flow rate and, consequently, the Re and We on the velocity contours and the formulation of recirculation zones. In this study, the steady state solutions are considered. It is found that the time required for obtaining the steady solution is dependent on the Reynolds number. For example, the number of time steps for steady state for Q = 0.1 cm 3 /s is 10,000 and the time step in this case is 0.0001 while the number of time steps is 50,000 for Q = 1 cm 3 /s at the same time step. Figure 4 shows low Re flows in which there are no vortices as indicated in cases (a) and (b), and the main stream fills the aneurysm. The vortices begin to appear for higher Re, and it is noted that the vortex size increases with the increase of flow rate till it fills the entire aneurysm. The formulation of vortices proves to have a great effect on the integrity of the arterial wall. 26

WSS distribution
The WSS, defined as the sum of the Newtonian contribution and the polymeric contribution, is one of the most important factors that affects the integrity of the arterial wall is the WSS. Boyd et al. 26 found that the minimum values of the WSS have more substantial effect than large values. Figure 5 shows the spatial variation of the WSS for different flow rates. The minimum values for WSS occur in the region of the   aneurysm due to the convective deceleration of the flow which results in small velocity gradients at the wall. 13 The WSS changes its sign from negative values within the aneurysm to positive sign at the distal end of the aneurysm at which it reaches its maximum values.

Wall pressure distribution
Another important factor in hemodynamics is the blood pressure distribution at the arterial wall shown in Figure 6 for different flow conditions. The pressure is decreasing downstream till it reaches the aneurysm region in which it starts to build up until it reaches it maximum, which is proportional to the flow rate; at the distal end, it starts to decrease again toward the artery exit.

Axial velocity
To have more insights into the blood flow through the abdominal aneurysms, it is convenient to calculate the axial velocity on a vertical line at the mid-section of the aneurysm for different flow rate conditions as shown in Figure 7. It is noted that the axial velocity increases with the increase of the Reynolds number due to the increase of the flow rate. The axial velocity takes negative values near the arterial wall, and these values appear clearly for high flow rates. This is due to the formation of the vortices which assert again that the aneurysm zone is a region which is worth to be studied.

Viscoelastic versus Newtonian flows
To show the importance of the inclusion of the viscoelasticity effects, the WSS is compared for both cases for different flow conditions. Figure 8 shows a difference of about 13% for the low flow rate case, while this difference reaches 40% in the distal end region of the aneurysm for higher flow rates as indicated in Figure 9.
This asserts that the effect of red blood cell in the blood flow simulation is important and should not be ignored, although there are many studies that treated the blood as a Newtonian fluid. This conclusion is consistent with Berger and Jou. 27

Summary and conclusion
In this contribution, mathematical modeling and finite element analysis are presented for viscoelastic Oldroyd-B fluid. The model is used to simulate blood flow in a rigid artery. In order to simulate low shear rate conditions as well as high shear rates, different flow conditions are considered ranging from volumetric flowrate of 0.1-3 cm 3 /s. Results show that there are no vortices observed at high shear rates (low volumetric flow rates), and hence, the difference between the Newtonian and viscoelastic solution is about 13% and the Newtonian assumption could be safely used keeping in mind the associated difference. Vortices are observed at higher volumetric flow rates which correspond to low shear rate situations. The relative error in this case is of order 40% and cannot be neglected, and hence, the effect of red blood cells sounds at the macroscale. However, the Oldroyd-B model is limited for low Reynolds number flow and still needs some modifications to account for shear dependent viscosity or shear-thinning property.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.