A Hybrid Stabilized Finite Element/Finite Difference Method for Unsteady Viscoelastic Flows
Ahmed Elhanafy1, Amr Guaily2, Ahmed Elsaid1
1Mathematics & Engineering Physics Department, Faculty of Engineering, Mansoura University, Mansoura, Egypt
2Engineering Mathematics & Physics Department, Faculty of Engineering, Cairo University, Giza, Egypt
To cite this article:
Ahmed Elhanafy, Amr Guaily, Ahmed Elsaid. A Hybrid Stabilized Finite Element/Finite Difference Method for Unsteady Viscoelastic Flows. Mathematical Modelling and Applications. Vol. 1, No. 2, 2016, pp. 26-35. doi: 10.11648/j.mma.20160102.11
Received: September 10, 2016; Accepted: October 14, 2016; Published: October 21, 2016
Abstract: The Oldroyd-B constitutive equation is used for the numerical simulation of unsteady incompressible viscoelastic flows. A novelty treatment is presented for the incompressibility constraint of the incompressible viscoelastic flow by using the modified continuity equation which allows using equal-order interpolation polynomials for all variables. The proposed technique circumvents the so-called LBB compatibility condition without pressure checkerboard and the solution instabilities with less computational costs compared with the traditional techniques. The discrete elastic-viscous stress-splitting method (DEVSS) is used to treat the instabilities resulting from the numerical simulation of viscoelastic flows. Two benchmark problems are simulated, namely, the flow through a channel with a bump and the flow inside a square cavity. Solutions are obtained for different Weissenberg number values and the results are compared with the published works.
Keywords: Unsteady Incompressible Viscoelastic Flow, Oldroyd-B Model, Pressure Stabilization Technique, The DEVSS Method, Galerkin Least Squares
Numerical simulation of incompressible viscoelastic flows is considered one of the important subjects of research for the manufacturing industry of polymer materials in the last few decades. However, there are many challenges in constructing numerical techniques for solving incompressible viscoelastic flow problems in complex geometry especially at high Weissenberg number . One of the most numerical challenges of the incompressible viscoelastic flow simulation is the incompressibility constraint since the continuity equation becomes a constraint equation for the velocity field rather than being an evolution equation for the density field. Differences in various methods of solving the incompressible flow equations originate from differences in strategies for satisfying the incompressibility constraint .
In general, these methods can be classified into two types: the first is the weakly incompressible solvers, in which the incompressibility constraint is satisfied by extending the compressible flow solvers for low Mach number where the incompressible flow is the limiting case of the compressible flow when the Mach number reaches zero [3-5]. The second is to satisfy the incompressibility constraint directly using the pressure as a mapping parameter to satisfy the continuity equation. This is called the pressure-based approach  and it can be grouped into two sub-categories. The first category depends on deriving a Poisson equation for the pressure calculation by taking the divergence of the momentum equation. A typical example of this approach is the semi-implicit method for pressure-linked equations (SIMPLE) .
The second category for the pressure-based approach is the stabilization technique which depends on modifying the continuity equation by adding an auxiliary term for the pressure. There are three types for this category, namely, the penalty formulation, the artificial compressibility formulation and the pressure stabilization formulation. The first two types were not derived as stabilization methods, but were initiated from alternative mathematical and physical concepts . The first type, the penalty formulation, is an application of the penalty function concept of Courant  to the finite element analysis of incompressible flows. It is proposed by Temam  and Zienkiewicz . Hughes  has suggested this formulation for both steady and unsteady flow cases. The second type is the artificial compressibility formulation, it is originally introduced by Chorin  where an artificial partial-time derivative is added to the continuity equation and the pressure in the additional term is re-scaled by a parameter called artificial compressibility. The third type is the pressure stabilization technique which depends on adding a Laplacian term for pressure to the continuity equation . This formulation is adopted only for Newtonian flow simulations in literature. The novelty of this study is the adopting of the pressure stabilization technique for incompressible viscoelastic flow with the Oldroyd-B model.
Another challenge for the numerical simulation of the incompressible viscoelastic flow is that the stress tensor gradient in the momentum equation cannot be expressed, anymore, in terms of second-order derivatives of the velocity field which means that the momentum equation loses the elliptic term which appears in the Newtonian case. This leads to the appearance of instabilities in the numerical solution . There are many schemes for simulating viscoelastic flows such as EEME , EVSS  and DEVSS scheme . In this study, the DEVSS is used. In this scheme, the extra stress tensor is split into the viscous and elastic parts. For stability reasons, an extra elliptic term is introduced in the discrete version of the momentum equation. The convective terms presented in the governing equations are treated using the Galerkin least squares method . Numerical simulation is performed for two benchmark problems, namely, the lid-driven cavity problem using the Oldroyd-B model then the channel flow with a bump is considered. Solutions are obtained for different values of the Weissenberg number and the results are compared with the work of ,  and .
2. Governing Equations
The conservation of mass and linear momentum for unsteady, incompressible, isothermal, flow leads to
where is the density, is the velocity vector, is the isotropic pressure and is the extra stress tensor. For viscoelastic fluids, the extra stress tensor is related to the velocity gradient via a differential equation rather than an algebraic equation . In the Oldroyd-B model, the extra stress tensor is a linear superposition of two components: a viscoelastic component and a purely Newtonian contribution as follows
where is the solvent zero-shear-rate viscosity, is the polymer zero-shear-rate viscosity, is the relaxation time, and is the rate of the deformation tensor defined
From numerical point of view, the dimensionless formulation of the governing equation is preferable, so we introduce the following non-dimensional variables
where is the characteristic length and is the free stream velocity. For 2-D planar flow, the velocity vector has the form and the polymeric stress tensor is given by
where is the axial stress, is the shear stress and is the normal stress. Using the above non-dimensionalization scheme, Eqs. (1) and (5) can be rewritten in a non-dimensional form as follows
where Re is the Reynolds number defined as , is the total viscosity given by , We is the Weissenberg number defined as and is the viscosity ratio defined by .
2.1. Incompressibility Constraint
When the Galerkin finite elements method is applied to the governing equations (1-2), for Newtonian flow for simplicity, zeros appear in the main diagonal of the resulting algebraic system and hence the number of the velocity unknowns must be greater than the pressure unknowns to obtain a unique solution . For finite difference and finite volume spatial discretization, staggered grids where pressures and velocities are calculated at different points are used to obtain stable solutions. For finite elements discretization, the approximation spaces must satisfy the so-called LBB condition  in which the order of the pressure interpolation functions must be lower than the order of the velocities interpolation functions. Although the elements that satisfy the LBB condition provide stable solutions, it has some complications in the finite elements programming. To avoid this problem and in the same time avoid the pressure spurious oscillations, equal-order interpolation functions for the pressure and the velocities are used by modifying the continuity equation by adding additional stabilization terms. In this method , the incompressibility constraint is replaced by
Where is the stabilization parameter.
2.2. The DEVSS method
In this method, the discrete form of the linear momentum balance equation is modified by introducing an auxiliary variable as a discrete approximation of the rate of deformation tensor where which leads to a symmetric matrix system . The linear momentum balance equation (1) could be rewritten as follows
Note that if the exact solution is recovered, then the elliptic operator vanishes. However, in the finite element calculation this is generally not the case  and the extra variable can be calculated from the weighted residual formulation of the following non-dimensional equation
We can rewrite Eq. (11) in the non-dimensional form as follows
The above equation is the stabilized form of the momentum conservation equation re-formulated by the DEVSS method. Eqs. (13), (10) and (9) constitute the governing equations needed for the simulation of the viscoelastic flow. The boundary conditions are discussed separately under each problem.
3. Finite Elements Formulation
In this study, the standard Galerkin formulation is used for all terms except for the convective terms in the conservation momentum and constitutive equations where the Galerkin least squares (GLS) formulation is used, in which the shape function is modified by adding a stabilized term .
3.1. Spatial Discretization
Since the modified continuity equation is of elliptic type, all variables, namely, can be interpolated by equal-order interpolation shape functions. The weak form of the governing equations over each element is achieved by considering the weighted residual approach and can be written as follows
The weak form of the momentum equation is
and the weak form of the constitutive equation is
and the weak form of the modified continuity equation is
For Eq. 12, the weak form is
where is the standard Galerkin test function and is the modified test function for convective terms and that is given by
where is a differential operator, the quantities and are the average values of the velocity components calculated at the center of the element and is the stabilization parameter. There are many suggestions for the stabilization parameter . Tezduyar  has proposed the following relation for unsteady flow
Where is the velocity over an element and is the characteristic length of the element and is the kinematic viscosity. Approximating the variables over each element with the proper function and applying integration by parts to the second order derivative terms and the first order derivative terms of the pressure, yield the following elemental equations
with the following matrices
where are the average values of polymeric stress components over the element. All integrals are evaluated numerically using 3-points Gauss-Legender quadrature method.
3.2. Time Marching and Solution Algorithm
The time derivative terms are discretized by the first-order explicit scheme yielding a system of linear algebraic equations. The time step is chosen in such a way to satisfy the CFL number to obtain a stable solution. Lumping approximation is applied to the consistent mass matrix converting it to a diagonal matrix . The solution scheme is summarized in the following steps:
1. Set zero initial values for the variables at time level n.
2. Impose the proper boundary conditions.
3. Solve Eqs. (22-26) to calculate at time level n+1.
4. Calculate the pressure at time level n+1 from Eq. (22).
5. Calculate the auxiliary variable at time level n+1 from Eq. (17).
4. Numerical Results
To validate the proposed algorithm, it is applied to the following benchmark problems. Solutions are obtained at different values for the Weissenberg number and the results are compared with published works.
4.1. Flow Through a Channel with a Bump
First we consider the unsteady incompressible viscoelastic flow through a planar channel with a circular arc bump as shown in Fig. 1. This problem is selected due to the presence of high stress gradients especially in the vicinity of the bump . The channel has length of 4H where H is the channel width. The curved boundary of the lower wall is given by
where , and is the height-to-length ratio.
4.1.1. Boundary Conditions
No-slip boundary conditions for velocity components are imposed at the solid walls. At inlet the velocity profile is assumed to be parabolic and the axial velocity component is given by
and the normal velocity component is zero. Pressure is specified to be zero at exit and Neumann boundary condition for the pressure is imposed for the continuity equation as it is of elliptic type. Because of the hyperbolic nature of the constitutive equations the stress components must be specified at the entrance. They are calculated using the assumed velocity profile at the inlet from the following relations 
At the solid walls, the stress components are obtained by solving the constitutive equations .
4.1.2. Results and Discussion
We first solve this problem for the Weissenberg number using a grid of 15x60 elements as shown in Fig. 1. In this simulation, we take , the time step , and the pressure dissipation parameter is . Axial velocity, normal velocity and pressure contours are shown in Fig. 2, Fig. 3 and Fig. 4, respectively, whereas the axial stress, shear stress and normal stress are shown in Fig. 5, Fig. 6 and Fig. 7, respectively.
More results are obtained for different values for the Weissenberg number to show its effect especially on the stress components. A finer grid (20x80) is used for Weissenberg number higher than 0.3 to capture the high stress gradients. Fig. 8, Fig. 9 and Fig. 10 show the axial stress, shear stress and normal stress for We = 0.7, respectively.
The effect of the Weissenberg number is clearly shown in Fig. 11 which shows the effect of the Weissenberg number on the lower wall shear stress especially in the vicinity of the bump. We note that the wall shear stress increases at the entrance of the bump. This increase is due to the decrease of the flow area. With the increase of the flow area, the wall shear stress increases again to reach its constant values after the bump. The increase of the Weissenberg number increases the extreme values of the wall shear stress in the vicinity of the bump.
4.2. Lid-driven Cavity Flow
In this standard problem the flow is confined in a unite square bounded by solid walls with moving the top boundary to the right. The main challenge associated to this problem is the level of stress growth near the corners between the lid and the solid walls especially for high Weissenberg numbers . The singularity of the flow field at the upper corners causes the numerical scheme to diverge, therefore we chose the regularized boundary condition for the upper wall proposed in  which is given by the following relation
where is a preset threshold to smooth the time derivative and is the distance of the current node from the upper left corner.
With this condition, the discontinuity of the velocity field at the top upper corners is removed. All boundary conditions for velocity components are of Dirichlet type. Pressure is specified to be zero at the lower left corner and its normal derivative is set to zero at the walls. Boundary conditions for stress components are obtained by solving the constitutive equations at the solid walls as discussed in . The Reynolds number is based on the side length of the cavity and the speed of the lid. In this simulation, , and the time step to compare the results ,  and . Quadrilateral elements are chosen for finite element discretization and all variables are interpolated by equal order shape functions. The mesh used for our simulation is shown in Fig. 12.
Fig. 13. The steady streamlines for and Re =0.5.
First simulation is performed at . We take the pressure dissipation parameter is . The flow reaches the steady state after
Fig. 13 shows the steady state streamlines and the center of the vortex is located at (0.465, 0.795) which is consistent with . The axial stress, shear stress, normal stress and pressure contours are shown in Fig. 14, Fig. 15, Fig. 16 and Fig. 17, respectively.
To show the effect of Weissenberg number on the solution, results are obtained at different values of Wiesenberg number from to . For Newtonian fluids, i.e. zero Weissenberg number, symmetric streamlines are obtained and the primary vortex is located at (0.495, 0.78) as shown in Fig. 18. These results are consistent with the theoretical studies in  and experimental studies by Pakdel et al.  using laser Doppler velocimetry (LDV) and digital particle image velocimetry (DPIV). However, for viscoelastic fluids, as the value of increases the primary vortex shifts in the direction of the upstream towards the left corner.
Fig. 18. The steady streamlines for at Re = 0.5.
Fig. 19. The steady streamlines for at Re =0.5.
This result is indicated in Fig. 19 where the primary vortex is located at (0.45, 0.8). Fig. 20 and Fig. 21 show the effect of on the horizontal velocity component at and vertical velocity component and , respectively. With the increase in , the maximum value of the velocity components decreases as expected .
Fig. 20. Horizontal velocity component for different values of We at .
Fig. 21. Vertical velocity component for different values of We at .
5. Summary and Conclusion
Numerical simulation of unsteady incompressible viscoelastic Oldroyd-B fluids using the finite elements method is presented. Numerical challenges of incompressible viscoelastic flow simulation, namely, the loss of the elliptic term in the momentum equation and the incompressibility constraint as well as the convective nature of the momentum equation are successfully treated. For the incompressibility constraint problem, the continuity equation is modified by adding a Laplacian term and hence it is of elliptic type. With this modification, equal-order interpolation functions are used for all variables without pressure checkerboard problem. This is not the only advantage for this formulation; it is used also for the steady and unsteady case as well with lower computational costs compared with the traditional formulations such as the artificial compressibility method. Due to the convective nature of the momentum equation, Galerkin least squares method is used by adding a stabilized term to the shape function. Also, the DEVSS method is used as a stabilized method recovers the elliptic term in the momentum equation.
For validation purposes, the proposed algorithm is applied to two benchmark problems. The first problem is the flow through a planar channel with a bump. Solutions are obtained at different values for the Weissenberg number to show its effect on the stress components with acceptable degree of success. We note that the wall shear stress extreme values increase with the increase of the Weissenberg number in the vicinity of the bump while this increase has a weak effect on the wall shear stress after the bump. This remark is useful especially for the simulation of the blood flow in stenosed arteries where the lower wall shear stress values are in the stenosis region. The second test case is the lid-driven cavity flow. Although this problem is simple in its geometry, it has some numerical challenges due to the singularities in the top corners. Regularized boundary condition for the upper wall is used to avoid these singularities. The results for this problem are also compared with experimental and theoretical studies. Solutions are obtained for different values of the Weissenberg number to show the viscoelastic effect. For low Weissenberg number, i.e., the flow likes to be Newtonian and this is clearly seen in the streamlines and the main vortex location. When the Weissenberg number increases, the stress gradients increases especially at the corners which needs finer grids to capture the viscoelasticity effect. We note also the viscoelastic effect on the streamlines and the main vortex location where it moves towards the left corner.