Mathematical Modelling and Applications
Volume 1, Issue 2, December 2016, Pages: 26-35

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

Email address:

(A. Elhanafy)
(A. Guaily)
(A. Elsaid)

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


1. Introduction

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 [1]. 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 [2].

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 [1] 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) [6].

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 [7]. The first type, the penalty formulation, is an application of the penalty function concept of Courant [8] to the finite element analysis of incompressible flows. It is proposed by Temam [9] and Zienkiewicz [10]. Hughes [11] 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 [12] 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 [7]. 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 [13]. There are many schemes for simulating viscoelastic flows such as EEME [14], EVSS [15] and DEVSS scheme [16]. 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 [17]. 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 [1], [18] and [19].

2. Governing Equations

The conservation of mass and linear momentum for unsteady, incompressible, isothermal, flow leads to

(1)

(2)

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 [20]. 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

(3)

(4)

(5)

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

(6)

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

(7)

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

(8)

(9)

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 [21]. 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 [22] 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 [22], the incompressibility constraint is replaced by

(10)

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 [23]. The linear momentum balance equation (1) could be rewritten as follows

(11)

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 [24] and the extra variable  can be calculated from the weighted residual formulation of the following non-dimensional equation

(12)

We can rewrite Eq. (11) in the non-dimensional form as follows

(13)

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 [17].

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

(14)

and the weak form of the constitutive equation is

(15)

and the weak form of the modified continuity equation is

(16)

For Eq. 12, the weak form is

(17)

where  is the standard Galerkin test function and  is the modified test function for convective terms and that is given by

(18)

(19)

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 [25] has proposed the following relation for unsteady flow

(20)

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

(21)

(22)

(23)

(24)

(25)

(26)

with the following matrices

(27)

(28)

(29)

(30)

(31)

(32)

(33)

(34)

(35)

(36)

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 [26]. 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 [20]. The channel has length of 4H where H is the channel width. The curved boundary of the lower wall is given by

(37)

where ,  and  is the height-to-length ratio.

Fig. 1. Geometry and grid for the channel flow problem.

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

(38)

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 [27]

At the solid walls, the stress components are obtained by solving the constitutive equations [27].

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.

Fig. 2. Axial velocity contours for We =0.3 and Re = 1.0.

Fig. 3. Normal velocity contour for We = 0.3 and Re = 1.0.

Fig. 4. Pressure contour for We = 0.3 and Re = 1.0.

Fig. 5. Axial stress contours for We= 0.3 and Re = 1.0.

Fig. 6. Shear stress contours for We = 0.3 and Re = 1.0.

Fig. 7. Normal stress contours for We = 0.3 and Re = 1.0.

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.

Fig. 8. Axial stress contours for We = 0.7 and Re = 1.0.

Fig. 9. Shear stress contours for We = 0.7 and Re = 1.0.

Fig. 10. Normal stress contours for We = 0.7 and Re = 1.0.

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 [19]. 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 [1] which is given by the following relation

(39)

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 [28]. 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 [1], [18] and [19]. 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. 11. Lower wall shear stress (LWSS) for different values of We and Re =1.0.

Fig. 12. The grid for the cavity problem.

Fig. 13. The steady streamlines for  and Re =0.5.

Fig. 14. Axial stress contours for We= 0.5 and Re =0.5.

Fig. 15. Shear stress contours for We= 0.5 and Re =0.5.

Fig. 16. Normal stress contours for We= 0.5 and Re =0.5.

Fig. 17. Pressure contours for We= 0.5 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 [1]. 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 [27] and experimental studies by Pakdel et al. [28] 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 [27].

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.


References

  1. S.zou, X.Xu, J.Chen, X.Guo and Q.Wang, Benchmark numerical simulations of viscoelastic fluid flows with an efficient integrated lattice Boltzmann and finite volume scheme, Advances in mechanical engineering, Hindawi Publishing Corporation 2, 2014.
  2. K.Kwak, C.Kiris, J. Chang, Computational challenges of viscous incompressible flows, Journal of computers and fluids.34, 283-299, 2005.
  3. G.X.Xu, E.Li, V.Tan and G.R.Liu, Simulation of steady and unsteady incompressible flow using gradient smoothing method(GSM), Journal of computers and structures 90-91, 131-144, 2012.
  4. D.Choi, C.Merkle, The application of preconditioning in viscous flows, Journal of ComputPhys 105, 203–226, 1993.
  5. E.Turkel, Preconditioned methods for solving the incompressible and low speed compressible equations, Journal of ComputPhys.72, 277–375, 1987.
  6. S.V.Patankar, Numerical heat transfer and fluid flow, New York: Hemisphere Publishing, 1980.
  7. H.P.Langtangen, K.Mardal, R.Winther, Numerical methods for incompressible viscous flow, Journal of advances in water resources25, 1125-1146, 2002.
  8. R. Courant, Calculus of Variations and Supplementary Notes and Exercise. New York University, New York, 1956.
  9. R. Temam, Navier-Stokes Equations, North-Holland, Amsterdam, 1984.
  10. R. S. Marshall, J. C. Heinrich, and O. C. Zienkiewicz, Natural convection in a square enclosure by a finite element penalty function method using primitive fluid variables, Journal of Numerical Heat Transfer.1, 315–330, 1987.
  11. T. J. R. Hughes, W. K. Liu, and A. Brooks, Finite element analysis of incompressible viscous flows by the penalty function formulation, Journal of Computational Physics.30, 1–60, 1979.
  12. A. J. Chorin. A numerical method for solving incompressible viscous flow problems, Journal of Computational Physics, 212–26, 1967.
  13. X.Li, X.Han, X. wang, Numerical modeling of viscoelastic flows using equal low-order finite elements, J.Comput. Methods Appl. Mech. Engrg.199, 570-581,2010.
  14. S.R.Burdette, P.J.Coates, R.C.Armstrong, R.A.Brown, Calculations of viscoelastic flow through an axisymmetric corrugated tube using the explicitly elliptic momentum equation formulation, J.Non-Newtonian fluid Mech.33, 1-23, 1989.
  15. D.Rajagopalan, R.C.Armstrong, R.A.Brown, Finite element methods for calculation of steady viscoelastic flow using constitutive equations with a Newtonian viscosity, J.Non-Newtonian fluid Mech.36,159-192, 1990.
  16. R.Guenette, M.Fortin, A new mixed finite element method for computing viscoelastic flow, J.Non-Newtonian fluid Mech.60, 27-52,1995.
  17. T. J. R. Hughes, L.P.Franca and G.M.Hulbert, A new finite element formulation for computational fluid dynamics: VIII. The Galerkin/Least-squares method for advective-diffusive equations, Journal of computer methods in applied mechanics and engineering.73, 173-189, 1989.
  18. J.Hao, T.Pan, Simulation for high Weissenberg number viscoelastic flow by a finite element method, Applied mathematics letters 20,988-993, 2007.
  19. J.Su, J.Ouyang, X.Wang, B.Yang, W.Zhou, Lattice Boltzmann method for the simulation of viscoelastic fluid flows over a large range of Weissenberg number, J.Non-Newtonian fluid Mech.194,42-59, 2013.
  20. A.Guily, Mathematical modeling and numerical simulation of viscoelastic liquids, PhD. Thesis, Faculty of graduate studies, Calgary University, 2010.
  21. D.Kuzmin,J. Hmlinen, Finite Element Methods for Computational Fluid Dynamics: A Practical Guide, SIAM, 267-270, 2015.
  22. F. Brezzi, J. Pitkäranta, On the stabilization of finite element approximations of the Stokes problem, Efficient Solutions of Elliptic Systems(W. Hackbusch, ed.), Notes on Numerical Fluid Mechanics10, 11-19,Vieweg, Braunschweig, 1984.
  23. M.A.Hulsen, R.Fattal, R.Kupferman, Flow of viscoelastic fluids past a cylinder at high weissenberge number: stabilized simulations using matrix logarithms, J.Non-Newtonian fluid Mech.127, 27-93, 2005.
  24. Frank.P.T.Baaijens, Mixed finite element methods for viscoelastic flow analysis: a review, J.Non-Newtonian fluid Mech.79, 361-385, 1998.
  25. T.E.Tezduyar Stabilized finite element formulations for incompressible flow computations, Journal of Advances in applied mechanics 28, 1-44, 1991.
  26. A.N.Brooks, T.J.R. Hughes, Streamline upwind/ petrov-Galerkin formulations for convective dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Journal of computer methods in applied mechanics and engineering 32, 199-259, 1982.
  27. K.Yapici, B.Karasozen, Y.Uludag, Finite volume simulation of viscoelastic laminar flow in a lid-driven cavity, J.Non-Newtonian fluid Mech.164, 51-65, 2009.
  28. P. Pakdel, S.H. Spiegelberg, G.H. McKinley, Cavity flows of elastic liquids: two-dimensional flows, Phys. Fluids 9 (11) 3123–3140, 1997.

Article Tools
  Abstract
  PDF(1724K)
Follow on us
ADDRESS
Science Publishing Group
548 FASHION AVENUE
NEW YORK, NY 10018
U.S.A.
Tel: (001)347-688-8931