Advances in Mathematical Physics
Volume 2017 (2017), Article ID 5210708, 10 pages
https://doi.org/10.1155/2017/5210708

Research Article

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

R. Mahmood,¹ N. Kousar,¹ M. Yaqub,² and K. Jabeen²

¹Department of Mathematics, Air University, Islamabad, Pakistan
²Department of Mathematics, Riphah International University, Islamabad, Pakistan

Correspondence should be addressed to N. Kousar; nabeela@mail.au.edu.pk

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 as

where 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 becomes

in 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.

Figure 1: Sequence of grids on space mesh level: 1, 2, and 3 (from left to right).

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.

Table 1: Mesh statistics at different refinement levels.

Refinement level
Number of elements
Degrees of freedom

3
16
96

4
64
352

5
256
1344

6
1024
5248

7
4096
20736

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.

Figure 2: Stream function contours for different values of Bingham number .

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.

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

(a)

(b)

Figure 4: Vertical cut-lines at for velocity and velocity for different values of Bingham number.

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.

Figure 5: Viscosity contours for different values of Bingham numbers .

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.

Table 2: Kinetic energy at different refinement levels for different Bingham numbers.

Level

2

3

4

5

6

7

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.

Figure 6: Stream function contours for double lid driven cavity at different values of Bingham number .

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

(a)

(b)

Figure 8: Vertical cut-lines at for velocity and velocity for different values of Bingham number.

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.

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

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.

Table 3: Kinetic energy for double lid driven cavity at different Bingham numbers Bn.

Level

2

3

4

5

6

7

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

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

Level 2
Level 5
Level 7
Level 2
Level 5
Level 7
Level 2
Level 5
Level 7

0.00
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000
0.000000

0.05
−0.013251
−0.0010289
−0.0010116
−0.013028
−0.0005843
−0.0005535
−0.012885
−0.0003243
−0.0002947

0.10
−0.026501
−0.0020266
−0.0019455
−0.026055
−0.0011590
−0.0010759
−0.025771
−0.0006451
−0.0005783

0.15
−0.039752
−0.0033673
−0.0032162
−0.039083
−0.0020363
−0.0018627
−0.038656
−0.0012314
−0.0010568

0.20
−0.053003
−0.0059797
−0.0057659
−0.052111
−0.0044433
−0.0041157
−0.051542
−0.0035617
−0.0033013

0.25
−0.066254
−0.0113470
−0.0122350
−0.065138
−0.0105230
−0.0112940
−0.064427
−0.0103240
−0.0109420

0.30
−0.079504
−0.0237110
−0.0240660
−0.078166
−0.0234810
−0.0238490
−0.077312
−0.0234520
−0.0238150

0.35
−0.092755
−0.0389840
−0.0395680
−0.091194
−0.0391150
−0.0397450
−0.090198
−0.0392330
−0.0398780

0.40
−0.106010
−0.0557910
−0.0567500
−0.104220
−0.0561110
−0.0571140
−0.103080
−0.0563100
−0.0573210

0.45
−0.119260
−0.0728560
−0.0739620
−0.117250
−0.0732490
−0.0743940
−0.115970
−0.0734660
−0.0746230

0.50
−0.132510
−0.0889990
−0.0899600
−0.130280
−0.0894080
−0.0904000
−0.128850
−0.0896280
−0.0906270

0.55
−0.019256
−0.1023800
−0.1037600
−0.017249
−0.1027800
−0.1041800
−0.015968
−0.102980
−0.1044000

0.60
0.093994
−0.1128200
−0.1146200
0.095779
−0.113200
−0.1150200
0.096917
−0.113400
−0.1152300

0.65
0.207240
−0.1197600
−0.1218100
0.208810
−0.120130
−0.122200
0.209800
−0.120320
−0.1224000

0.70
0.320500
−0.1225400
−0.1244700
0.321830
−0.122910
−0.124860
0.322690
−0.123100
−0.1250500

0.75
0.433750
−0.1203500
−0.1212400
0.434860
−0.120740
−0.121650
0.435570
−0.120950
−0.1218600

0.80
0.547000
−0.1078300
−0.1082700
0.547890
−0.108310
−0.108750
0.548460
−0.108550
−0.1089800

0.85
0.660250
−0.0405360
−0.0578980
0.660920
−0.041105
−0.058487
0.661340
−0.041380
−0.0587660

0.90
0.773500
0.1480600
0.1287800
0.773940
0.147500
0.128200
0.774230
0.147230
0.1279300

0.95
0.886750
0.4934200
0.4848400
0.886970
0.493040
0.484440
0.887110
0.492860
0.4842500

0.97
0.932050
0.6960500
0.6724000
0.932180
0.695820
0.672130
0.932270
0.695720
0.6720100

0.98
0.954700
0.7973700
0.7766300
0.954790
0.797220
0.776440
0.954850
0.797140
0.7763500

0.99
0.977350
0.8986800
0.8862300
0.977390
0.898610
0.886130
0.977420
0.898570
0.8860800

1.00
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000

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

y

Level 2
Level 5
Level 7
Level 2
Level 5
Level 7
Level 2
Level 5
Level 7

0.00
0.00000
0.00000
0.00000
0.000000
0.000000
0.0000000
0.0000000
0.0000000
0.0000000

0.05
−0.01276
−0.0004015
−0.00040
−0.012632
−0.0002256
−0.0002200
−0.012567
−0.0001239
−0.0001164

0.10
−0.02552
−0.0007936
−0.00078
−0.025265
−0.0004461
−0.0004284
−0.025134
−0.0002451
−0.0002267

0.15
−0.03828
−0.0012214
−0.00120
−0.037897
−0.0006895
−0.0006578
−0.037701
−0.0003803
−0.0003491

0.20
−0.05104
−0.0017347
−0.00170
−0.050530
−0.0009877
−0.0009424
−0.050268
−0.0005491
−0.0005036

0.25
−0.06379
−0.0023872
−0.00238
−0.063162
−0.0013773
−0.0013403
−0.062835
−0.0007745
−0.0007269

0.30
−0.07655
−0.0034638
−0.00341
−0.075795
−0.0020713
−0.0020025
−0.075402
−0.0012027
−0.0011303

0.35
−0.08931
−0.0054841
−0.00524
−0.088427
−0.0037121
−0.0034697
−0.087969
−0.0026646
−0.0023786

0.40
−0.10207
−0.0100510
−0.00902
−0.101060
−0.0088888
−0.0077991
−0.100540
−0.0083352
−0.0073288

0.45
−0.11483
−0.0176790
−0.01558
−0.113690
−0.0178200
−0.0156190
−0.113100
−0.0181320
−0.0158290

0.50
−0.12759
−0.0277930
−0.02323
−0.126320
−0.028356
−0.0237720
−0.125670
−0.0287160
−0.0241450

0.55
−0.01483
−0.0369840
−0.03037
−0.013692
−0.037644
−0.0309950
−0.013103
−0.0380350
−0.0313710

0.60
0.09792
−0.0448550
−0.03652
0.098941
−0.045507
−0.0371250
0.099464
−0.0458770
−0.0374760

0.65
0.21069
−0.0512220
−0.04149
0.211570
−0.051834
−0.0420640
0.212030
−0.0521730
−0.0423870

0.70
0.32345
−0.0559550
−0.04517
0.324210
−0.056536
−0.0457100
0.324600
−0.0568550
−0.0460110

0.75
0.43621
−0.0589840
−0.04743
0.436840
−0.059545
−0.0479510
0.437170
−0.0598500
−0.0482400

0.80
0.54896
−0.0595660
−0.04808
0.549470
−0.060137
−0.0486150
0.549730
−0.0604460
−0.0489040

0.85
0.66172
−0.0579270
−0.04688
0.662100
−0.058539
−0.0474550
0.662300
−0.0588640
−0.0477610

0.90
0.77448
−0.0594230
−0.04293
0.774740
−0.060227
−0.0436020
0.774870
−0.0606400
−0.0439500

0.95
0.88724
0.1492400
0.02744
0.887370
0.148400
0.0257970
0.887430
0.1479700
0.0250270

0.97
0.93234
0.4895400
0.26696
0.932420
0.489040
0.2649700
0.932460
0.4887800
0.2640800

0.98
0.95490
0.6597000
0.47095
0.954950
0.659360
0.4693400
0.954970
0.6591900
0.4686200

0.99
0.97745
0.8298500
0.71853
0.977470
0.829680
0.7176000
0.977490
0.8295900
0.7171900

1.00
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000
1.000000

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.

Table 6: The number of nonlinear iterations (level 7) required to reduce nonlinear defect up to 10⁻⁶.

Bn
# NL

10
63

20
72

50
78

100
83

500
88

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.

References

F. N. Shwedov, “La Rigidite de liquides,” Rapport Congr. Intern. Phys. Paris, vol. 1, pp. 478–486, 1900.

E. C. Bingham, Fluidity and Plasticity, McGraw-Hill, New York, NY, USA, 1922.

H. A. Barnes, “The yield stress—a review or ‘πανταρει’—everything flows?” Journal of Non-Newtonian Fluid Mechanics, vol. 81, no. 1-2, pp. 133–178, 1999.

A. Syrakos, G. C. Georgiou, and A. N. Alexandrou, “Solution of the square lid-driven cavity flow of a Bingham plastic using the finite volume method,” Journal of Non-Newtonian Fluid Mechanics, vol. 195, pp. 19–31, 2013.

A. Syrakos, G. C. Georgiou, and A. N. Alexandrou, “Performance of the finite volume method in solving regularised Bingham flows: inertia effects in the lid-driven cavity flow,” Journal of Non-Newtonian Fluid Mechanics, vol. 208-209, pp. 88–107, 2014.

E. Mitsoulis and T. Zisis, “Flow of Bingham plastics in a lid-driven square cavity,” Journal of Non-Newtonian Fluid Mechanics, vol. 101, no. 1-3, pp. 173–180, 2001.

M. Bercovier and M. Engelman, “A finite element method for incompressible non-Newtonian flows,” Journal of Computational Physics, vol. 36, no. 3, pp. 313–326, 1980.

T. C. Papanastasiou, “Flows of materials with yield,” Journal of Rheology, vol. 31, no. 5, pp. 385–404, 1987.

M. Fortin and R. Glowinski, Augmented Lagrangian Methods: Applications to the Numerical Solution of Boundary-Value Problems, North-Holland, Amsterdam, The Netherlands, 1983.

R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer, New York, NY, USA, 1984.

E. J. Dean, R. Glowinski, and G. Guidoboni, “On the numerical simulation of Bingham visco-plastic flow: old and new results,” Journal of Non-Newtonian Fluid Mechanics, vol. 142, no. 1–3, pp. 36–62, 2007.

S. Turek and R. Rannacher, “A simple nonconforming quadriletral stokes element,” Numerical Methods for Partial Differential Equations, vol. 8, no. 2, pp. 97–111, 1992.

P. Knobloch, “On Korn's inequality for non-conforming finite elements,” Technisshe Mechanik, vol. 20, pp. 205–214, 2000.

S. Turek, A. Ouazzi, and R. Schmachtel, “Multigrid methods for stabilized nonconforming finite elements for incompressible flow involving the deformation tensor formulation,” Journal of Numerical Mathematics, vol. 10, no. 3, pp. 235–248, 2002.

S. Turek and A. Ouazzi, “Unified edge-oriented stabilization of nonconforming FEM for incompressible flow problems: numerical investigations,” Journal of Numerical Mathematics, vol. 15, no. 4, pp. 299–322, 2007.

T. A. Davis, “Algorithm 832: UMFPACK V4.3—an unsymmetric-pattern multifrontal method,” ACM Transactions on Mathematical Software, vol. 30, no. 2, pp. 196–199, 2004.

S. Turek, FEATFLOW Finite Element Software for the Incompressible Navier-Stokes Equations: User Manual, Release 1.2, University of Dortmund, 2000, http://www.featflow.de.

U. Ghia, K. N. Ghia, and C. T. Shin, “High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method,” Journal of Computational Physics, vol. 48, no. 3, pp. 387–411, 1982.

C.-H. Bruneau and M. Saad, “The 2D lid-driven cavity problem revisited,” Computers and Fluids, vol. 35, no. 3, pp. 326–348, 2006.

J. Zhang, “Bifurcation of bingham streamline topologies in rectangular double-lid-driven cavities,” Journal of Applied Mathematics and Physics, vol. 2, no. 12, pp. 1069–1072, 2014.

G. R. Kefayati, “FDLBM simulation of magnetic field effect on non-Newtonian blood flow in a cavity driven by the motion of two facing lids,” Powder Technology, vol. 253, pp. 325–337, 2014.

S. Arun and A. Satheesh, “Analysis of flow behaviour in a two sided lid driven cavity using lattice Boltzmann technique,” Alexandria Engineering Journal, vol. 54, no. 4, pp. 795–806, 2015.