Advances in Mathematical Physics

Volume 2017, Article ID 5210708, 10 pages

https://doi.org/10.1155/2017/5210708

## Numerical Simulations of the Square Lid Driven Cavity Flow of Bingham Fluids Using Nonconforming Finite Elements Coupled with a Direct Solver

^{1}Department of Mathematics, Air University, Islamabad, Pakistan^{2}Department of Mathematics, Riphah International University, Islamabad, Pakistan

Correspondence should be addressed to N. Kousar; kp.ude.ua.liam@aleeban

Received 20 December 2016; Accepted 28 February 2017; Published 20 March 2017

Academic Editor: Xin Yu

Copyright © 2017 R. Mahmood et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

In this paper, numerical simulations are performed in a single and double lid driven square cavity to study the flow of a Bingham viscoplastic fluid. The governing equations are discretized with the help of finite element method in space and the nonconforming Stokes element is utilized which gives 2nd-order accuracy for velocity and 1st-order accuracy for pressure. The discretized systems of nonlinear equations are treated by using the Newton method and the inner linear subproablems are solved by the direct solver UMFPACK. A qualitative comparison is done with the results reported in the literature. In addition to these comparisons, some new reference data for the kinetic energy is generated. All these implementations are done in the open source software package FEATFLOW which is a general purpose finite element based solver package for solving partial differential equations.

#### 1. Introduction

The study of non-Newtonian fluids gained much importance in view of its extensive industrial and technological applications. Examples of such fluids include slurries, shampoo, paints, clay coating and suspensions, grease, cosmetic products, and body fluids. Among the non-Newtonian fluids, the viscoplastic fluids are those that exhibit a yield stress and thus combine the behavior of solids and liquids in different flow regimes. The idea of yield stress was first given by Shwedov [1]. Afterwards, Bingham [2] presented the flow shear diagram showing a linear relationship between stress and strain to explain the nature of plastic. A detailed review of viscoplastic behavior of materials is carried out by Barnes [3]. The behavior of such materials is like a solid (elastic or inelastic) below the certain value of shear stress and a liquid otherwise. This critical value of stress is termed as yield stress. Based on this fact the flow field of such materials is divided into unyielded (solid) and yielded (fluid) regions.

The simple version of the constitutive equation describing the viscoplasticity is that proposed by Bingham [2]where is the yield stress, is the plastic viscosity, is stress tensor, and is the rate of strain tensor given by is the velocity vector, and superscript denotes the transpose of the velocity-gradient tensor . The symbols and denote the magnitudes of the stress and rate-of strain tensors, respectively:Despite its simplicity, this model contains all the ingredients of viscoplastic materials, namely, a yield stress and a nonlinear variation of the effective viscosity. It is also observed while dealing with Bingham Model; three different zones in the flow domain can be identified(i)Shear zone .(ii)Plug zone and .(iii)Dead zone and .

All these regions are analyzed with great detail in the flow problems considered in this paper.

To date, a large number of research papers are devoted to simulate quantitatively and qualitatively these fluids; see, for example, [4–7]. Due to their yield stress property, the numerical solution of Bingham flows is difficult because the interface between the yielded and unyielded regions is a priori not known. To circumvent this difficulty, there are two approaches. One approach includes methods which approximate equation by a regularized constitutive equation, treating the whole material as fluid of variable viscosity and it is applicable throughout the domain. For such methods a very high value is assigned to the viscosity locally in unyielded regions as discussed by Bercovier and Engelman [7] and by Papanastasiou [8]. The Papanastasiou regularization introduces an exponential term to replace the discontinuous constitutive equation(1) by a single equation as applicable throughout the material. Here is a stress growth parameter, which should be “sufficiently” large so that the ideal Bingham behavior is approximated with satisfactory accuracy. In view of (4), the viscosity is given by that can be used over the entire flow domain.

The other approach includes methods which start by deriving variational inequalities and minimizing the rate of strain or maximizing the stress [9, 10]. A good review giving comparisons of numerical methods based on variational inequality approach is given in [11].

The rest of the paper is organized as follows. In Section 2, the mathematical formulations of the governing equations are presented. In Section 3, the numerical approach is presented. The numerical results for single and double lid driven cavity flow are presented by means of velocity, viscosity, and stream function plots in Section 4. Finally, Section 5 contains our concluding remarks.

#### 2. Mathematical Formulation

The general form governing the incompressible behavior of these fluids can be written aswhere is the velocity vector, is the pressure scaled by density , is the body forces, and is the stress tensor which is responsible for categorizing different fluids. Assuming the flow as steady-state, two-dimensional isothermal, and incompressible and using nondimensional variables , , and and some reference velocity and reference length as and , respectively, the dimensionless form of the governing equations in the absence of body forces becomesin which the nondimensionalization procedure introduces the following important dimensionless numbers: the Reynolds number and the Bingham number and is the dimensionless stress growth parameter, given as

The dimensionless viscosity is now given as

The higher the value of* M *is, the better (8) approximates the actual Bingham constitutive equation, in the yielded regions of the flow field (), and the higher the apparent viscosity is in the unyielded regions, making them behave approximately as solid bodies. For practical reasons though,* M* must not be so high as to cause convergence problems to the numerical methods used to solve the above equations [4, 5].

#### 3. Numerical Approach

Our numerical approach is based on the Galerkin weighted residual finite element method for the solution of the governing equations for the velocity and pressure. The equations are discretized with the help of finite element method using the nonconforming LBB stable Stokes element of 2nd-order accuracy for velocity and 1st-order accuracy for pressure. This quadrilateral Stokes element is based on “rotated” bilinear shape functions having four local degrees of freedom (DOFs) for velocity component and one DOF for a piecewise constant pressure approximation (see [12, 13] for details). The chosen nonconforming element requires additional stabilization for handling the deformation tensor formulation due to missing Korn’s inequality [13]. To this end we employ the standard edge oriented stabilization [14, 15] in our simulations. Newton’s method is applied to solve discrete nonlinear algebraic systems and the direct solver UMFPACK [16] is used as an inner linear solver. For grid refinement, we generate a sequence of grids by uniform refinement from a coarsest mesh. Starting from the mesh level , we generate grid of mesh level by dividing each quadrilateral cell of grid level into four new quadrilaterals connecting the mid points of opposite edges. Figure 1 shows the grid on level , respectively, and the mesh size on grid level is . The libraries of the open source software package FEATFLOW [17] are used in the simulations.