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

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 Q1/Q0 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.


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] γ = 0,  ≤    = (   γ + ) γ ,  >   , where   is the yield stress,  is the plastic viscosity,  is stress tensor, and γ is the rate of strain tensor given by γ ≡ ∇u + (∇u)  .
(2) u is the velocity vector, and superscript  denotes the transpose of the velocity-gradient tensor ∇u.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 (ii) Plug zone ‖‖ = 0 and  = constant.
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][5][6][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.

Mathematical Formulation
The general form governing the incompressible behavior of these fluids can be written as where u is the velocity vector,  is the pressure scaled by density , f 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 u * ,  * , and  * and some reference velocity and reference length as  ref and  ref , respectively, the dimensionless form of the governing equations in the absence of body forces becomes in which  (10) 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,  = [Bn/ γ + 1] γ in the yielded regions of the flow field ( > Bn), 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].

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 Q1 / 0 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 Q1 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  = 1, we generate grid of mesh level  + 1 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  = 1,2,3, respectively, and the mesh size on grid level  is ℎ = 2 −+1 .The libraries of the open source software package FEATFLOW [17] are used in the simulations.9)- (11).The mesh statistic for this benchmark problem is provided in Table 1. Figure 2 depicts the stream function contour snapshots which indicate the effect of yield stress on the primary vortex.It can be visualized that the main vortex shrinks and approaching to the upper lid with an increase in the value  of Bingham number Bn indicating that as the yield limit increases, the dead region increases and it occupies more of the cavity and the shear region is moved to be close to the moving lid.

Results and Discussions
Effect of Bingham number Bn on velocity profile is shown in Figure 3 which also confirms that at higher values of Bingham number Bn the velocity is nonzero only in the region close to the upper moving lid.The numeric data for these snapshots of velocity along a vertical centerline is presented in Figure 4.
In Figure 4(a), we plot horizontal velocity along a vertical center line.The curve Bn = 0 corresponds to Newtonian case.For all other cases, the yielded and unyielded zones can be identified by decomposing each curve into two segments.The horizontal segment with  = 0 from  = 0 up to a certain height corresponds to unyielded zone due to increase in Bingham number.The velocity is extremely low throughout the lower portion of the cavity for higher Bn numbers.The other segment corresponds to the upper yielded zone.Such behavior is also observed in [4,5] and is in good agreement concerning qualitative analysis.The vertical velocity components for all cases are plotted in Figure 4(b).
The dimensionless viscosity as a function of Bingham number Bn is displayed in Figure 5.The progressive growth of unyielded regions can be seen in the bottom of the cavity due to an increase in the plasticity effects produced at higher Bingham numbers Bn.
In addition to the local quantities like velocity, pressure, and viscosity we have also computed the kinetic energy in the cavity which is one of the global quantities of interest.Kinetic energy is defined by To see the Newtonian results of this benchmark quantity we refer to the results of Bruneau and Saad [19].Table 2 shows kinetic energy for the single lid driven cavity.The results are grid independent after level 7.It is also noted that an increase in Bingham number results in the decrease of kinetic energy due to the enhanced plasticity effect producing unyielded regions.geometry of the problem is again a unit square but now the upper and lower walls of the cavity are moving with the constant speed  =  = 1, V = 0 from left to right.Further studies on this problem can be found in [20][21][22].

Double Lid Driven Cavity
In case of double lid driven cavity the two primary vertices are generated as shown in the stream function plots in Figure 6.These vertices move along to the upper and lower walls with increasing the Bingham numbers.The velocity profile snapshots are presented in Figure 7 and vertical cutlines along the center line at  = 0.5 in Figure 8.
Figure 8(a) reflects the evolution of horizontal velocity along the vertical centerline and in Figure 8(b) the plots for vertical velocity along the same centerline are given.The profiles for horizontal velocity become more and more flat in region near to the center showing that the motion of lids cannot penetrate towards the center of cavity due to the enhanced plasticity effects at higher Bn numbers.
In Figure 9, the viscosity contours are presented to show the effect of Bingham numbers on the viscosity.In contrast to the single lid driven cavity, we now see the creation of unyielded regions in the center of the cavity due to distance from the source of motion.These snapshots also reveal the presence of side eddies in the center attached to the stationary walls.An increase in the Bn number limits the yielded region closer to the moving lids.
The kinetic energy is tabulated in Table 3.The convergence of the results can be seen by moving up to down in each column.After level 7 the results are the same which shows grid independency.It is also noted that an increase in Bingham number results in the decrease of kinetic energy due to the expansion in the dead zones.The -component of the velocity () is tabulated at the chosen points in Tables 4 and 5 for selected values of Bingham numbers to show the influence of the parameter  and refinement level on the solution.The present choices of  do not have a significant effect on the solution.It is also worth mentioning that the coarser grid level is a larger source of error than the smallness of .This observation is also noted in [4].
We close this section by showing the dependence of the maximum nonlinear iterations on Bingham number in Table 6.An increase in the number of nonlinear iterations (# NL) is observed with an increase in the non-Newtonian dimensionless parameter Bn.

Conclusions
This work presents an insight into the behavior of numerical simulations for Bingham viscoplastic fluid flows in benchmark configurations.The implementations are done via finite element methods in the framework of a monolithic approach.The results obtained for the Bingham Models are able to describe the viscosity function accurately for the configurations of single and double lid driven cavity flows.Results have been obtained for Bingham numbers in the range of 0-500 and are presented by means of velocity, viscosity, and stream function plots.Beside these local quantities, the tabular data for the kinetic energy in the cavity is also generated.It is also noted that number of iterations for

2
Advances in Mathematical Physics dealing with Bingham Model; three different zones in the flow domain can be identified (i) Shear zone ‖‖ ̸ = 0.

4. 1 .
Lid Driven Cavity Flow.Using the methodology described in the previous section, the lid driven cavity problem is simulated.This important benchmark problem is considered by many researchers[5,6,18].Consider a square cavity domain Ω = [0, 1]×[0, 1], filled with a Bingham fluid which is set to motion by the upper lid of the cavity which moves with a uniform horizontal velocity .The zero Dirichlet boundary conditions are given on all other walls.The velocity of the upper lid and length of the cavity are taken as  ref and  ref , respectively, appearing in (

Figure 3 :
Figure 3: Velocity profiles for different values of Bingham numbers Bn.

5 −
Flow.This section presents numerical results for double lid driven cavity flow.TheAdvances in Mathematical Physics

Figure 4 :Figure 5 :Table 2 :Figure 6 :
Figure 4: Vertical cut-lines at  = 0.5 for  velocity and V velocity for different values of Bingham number.

Figure 7 :Figure 8 :
Figure 7: Velocity profiles for double lid driven cavity at different values of Bingham number Bn.

Figure 9 :
Figure 9: Viscosity contours for double lid driven cavity at different values of Bingham number Bn.

Table 1 :
Mesh statistics at different refinement levels.

Table 4 :
Values of the -component of the velocity along the vertical center line for Bn = 10.

Table 5 :
Values of the -component of the velocity along the vertical center line for Bn = 100.