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.
4. Results and Discussions
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 , 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 and , respectively, appearing in (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.
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.
(a)
(b)
In Figure 4(a), we plot horizontal velocity along a vertical center line. The curve 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 from 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 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 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.
4.2. Double Lid Driven Cavity Flow
This section presents numerical results for double lid driven cavity flow. The 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 , from left to right. Further studies on this problem can be found in [20–22].
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 cut-lines along the center line at in Figure 8.
(a)
(b)
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.
5. 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 Newton’s solver increased at larger values of Bingham number due to enhanced nonlinearity.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.