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

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]

[mathematical expression not reproducible] (1)

where [[tau].sub.y] is the yield stress, [mu] is the plastic viscosity, [tau] is stress tensor, and [??] is the rate of strain tensor given by

[??] [equivalent to] [nabla]u + [([nabla]u).sup.T]. (2)

u is the velocity vector, and superscript T denotes the transpose of the velocity-gradient tensor [nabla]u. The symbols [tau] and [??] denote the magnitudes of the stress and rate-of strain tensors, respectively:

[mathematical expression not reproducible] (3)

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 [parallel]D[parallel] [not equal to] 0.

(ii) Plug zone [parallel]D[parallel] = 0 and u = constant.

(iii) Dead zone [parallel]D[parallel] = 0 and m = 0.

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(l) by a single equation as

[mathematical expression not reproducible] (4)

applicable throughout the material. Here m 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

[eta] = [[tau].sub.y]/[??]{1 - exp(-m[??])} + [mu] (5)

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

[mathematical expression not reproducible] (6)

where u is the velocity vector, p is the pressure scaled by density [rho], f is the body forces, and [tau] 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.sup.*], [p.sup.*], and [[tau].sup.*] and some reference velocity and reference length as [U.sub.ref] and [L.sub.ref], respectively, the dimensionless form of the governing equations in the absence of body forces becomes

[nabla] x [u.sup.*] = 0 Re [u.sup.*] x [nabla][u.sup.*] = - [nabla][p.sup.*] + [nabla] x [[tau].sup.*] (7)

in which

[mathematical expression not reproducible] (8)

the nondimensionalization procedure introduces the following important dimensionless numbers: the Reynolds number

Re [equivalent to] [rho][U.sub.ref] [L.sub.ref]/[mu] (9)

and the Bingham number Bn

Bn [equivalent to] [[tau].sub.y][L.sub.ref]/[mu][U.sub.ref] (10)

and M is the dimensionless stress growth parameter, given as

M [equivalent to] m[U.sub.ref]/[L.sub.ref] (11)

The dimensionless viscosity is now given as

[mathematical expression not reproducible] (12)

The higher the value of M is, the better (8) approximates the actual Bingham constitutive equation, [tau] = [Bn/[??] + 1][??] in the yielded regions of the flow field ([tau] > 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].

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 [Q.sub.1]/[Q.sub.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 [Q.sub.1] 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 l = 1, we generate grid of mesh level l +1 by dividing each quadrilateral cell of grid level I into four new quadrilaterals connecting the mid points of opposite edges. Figure 1 shows the grid on level l = 1,2,3, respectively, and the mesh size on grid level l is h = [2.sup.-l+1]. 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 [OMEGA] = [0,1] x [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 U. 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 [U.sub.ref] and [L.sub.ref], 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.

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 u = 0 from y = 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

E = [1/2] [[integral].sub.[OMEGA]] [[parallel]u[parallel].sup.2] dx. (13)

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 u = U = 1, v = 0 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 cutlines along the center line at x = 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 x-component of the velocity (u) is tabulated at the chosen points in Tables 4 and 5 for selected values of Bingham numbers to show the influence of the parameter M and refinement level on the solution. The present choices of M 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 M. 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.

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

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

References

[1] F. N. Shwedov, "La Rigidite de liquides," Rapport Congr. Intern. Phys. Paris, vol. 1, pp. 478-486, 1900.

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

[3] H. A. Barnes, "The yield stress--a review or '[phrase omitted]' --everything flows?" Journal of Non-Newtonian Fluid Mechanics, vol. 81, no. 1-2, pp. 133-178, 1999.

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

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

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

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

[8] T. C. Papanastasiou, "Flows of materials with yield," Journal of Rheology, vol. 31, no. 5, pp. 385-404, 1987.

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

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

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

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

[13] P. Knobloch, "On Korn's inequality for non-conforming finite elements," Technisshe Mechanik, vol. 20, pp. 205-214, 2000.

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

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

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

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

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

[19] C.-H. Bruneau and M. Saad, "The 2D lid-driven cavity problem revisited," Computers and Fluids, vol. 35, no. 3, pp. 326-348, 2006.

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

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

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

R. Mahmood, (1) N. Kousar, (1) M. Yaqub, (2) and K. Jabeen (2)

(1) Department of Mathematics, Air University, Islamabad, Pakistan

(2) 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

[Please Note: Some non-Latin characters were omitted in this article.]

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

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

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

Caption: Figure 4: Vertical cut-lines at x = 0.5 for u velocity and v velocity for different values of Bingham number.

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

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

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

Caption: Figure 8: Vertical cut-lines at x = 0.5 for u velocity and v velocity for different values of Bingham number.

Caption: Figure 9: Viscosity contours for double lid driven cavity at different values of Bingham number Bn.
```Table 1: Mesh statistics at different refinement levels.

Refinement     Number of      Degrees
level          elements       of freedom

3                   16             96
4                   64            352
5                  256            1344
6                  1024           5248
7                  4096          20736

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

Level   Bn = 0            Bn = 10           Bn = 20

2       5.285538E -02     4.285017-E -02    4.277462E -02
3       3.841554E -02     2.417936E -02     2.257387E -02
4       3.480833E -02     1.750518E -02     1.536717-E -02
5       3.365933E -02     1.579714E -02     1.218378E -02
6       3.359881E -02     1.540207-E -02    1.152689E -02
7       3.358246E -02     1.529102E -02     1.137576E -02

Level   Bn = 50           Bn = 100          Bn = 500

2       4.275709E -02     4.275291E -02     4.274996E -02
3       2.206525E -02     2.202012E -02     2.198380E -02
4       1.401530E -02     1.342773E -02     1.284529E -02
5       8.984129E -03     8.415416E -03     8.015731E -03
6       7.721842E -03     5.705813E -03     4.536025E -03
7       7.450573E -03     5.383305E -03     2.441167-E -03

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

Level   Bn = 0            Bn = 10           Bn = 20

2       9.422139E -02     8.561376E -02     8.560186E -02
3       7.213333E -02     5.170203E -02     4.909028E -02
4       6.531628E -02     3.784316E -02     3.265451E -02
5       6.357059E -02     3.438146E -02     2.649637E -02
6       6.314135E -02     3.359117-E -02    2.531207-E -02
7       6.302763E -02     3.337230E -02     2.502079E -02

Level   Bn = 50           Bn = 100          Bn = 500

2       8.559770E -02     8.559568E -02     8.559568E -02
3       4.675238E -02     4.580014E -02     4.505332E -02
4       2.988453E -02     2.879644E -02     2.773203E -02
5       1.926139E -02     1.748698E -02     1.667230E -02
6       1.686371-E -02    1.266649E -02     9.243272E -03
7       1.632972E -02     1.171277-E -02    5.139763E -03

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

y                       M = 100
Level 2        Level 5        Level 7

0.00      0.000000       0.000000       0.000000
0.05     -0.013251      -0.0010289     -0.0010116
0.10     -0.026501      -0.0020266     -0.0019455
0.15     -0.039752      -0.0033673     -0.0032162
0.20     -0.053003      -0.0059797     -0.0057659
0.25     -0.066254      -0.0113470     -0.0122350
0.30     -0.079504      -0.0237110     -0.0240660
0.35     -0.092755      -0.0389840     -0.0395680
0.40     -0.106010      -0.0557910     -0.0567500
0.45     -0.119260      -0.0728560     -0.0739620
0.50     -0.132510      -0.0889990     -0.0899600
0.55     -0.019256      -0.1023800     -0.1037600
0.60      0.093994      -0.1128200     -0.1146200
0.65      0.207240      -0.1197600     -0.1218100
0.70      0.320500      -0.1225400     -0.1244700
0.75      0.433750      -0.1203500     -0.1212400
0.80      0.547000      -0.1078300     -0.1082700
0.85      0.660250      -0.0405360     -0.0578980
0.90      0.773500       0.1480600      0.1287800
0.95      0.886750       0.4934200      0.4848400
0.97      0.932050       0.6960500      0.6724000
0.98      0.954700       0.7973700      0.7766300
0.99      0.977350       0.8986800      0.8862300
1.00      1.000000       1.000000       1.000000

y                       M = 200
Level 2        Level 5        Level 7

0.00      0.000000       0.000000       0.000000
0.05     -0.013028      -0.0005843     -0.0005535
0.10     -0.026055      -0.0011590     -0.0010759
0.15     -0.039083      -0.0020363     -0.0018627
0.20     -0.052111      -0.0044433     -0.0041157
0.25     -0.065138      -0.0105230     -0.0112940
0.30     -0.078166      -0.0234810     -0.0238490
0.35     -0.091194      -0.0391150     -0.0397450
0.40     -0.104220      -0.0561110     -0.0571140
0.45     -0.117250      -0.0732490     -0.0743940
0.50     -0.130280      -0.0894080     -0.0904000
0.55     -0.017249      -0.1027800     -0.1041800
0.60      0.095779      -0.113200      -0.1150200
0.65      0.208810      -0.120130      -0.122200
0.70      0.321830      -0.122910      -0.124860
0.75      0.434860      -0.120740      -0.121650
0.80      0.547890      -0.108310      -0.108750
0.85      0.660920      -0.041105      -0.058487
0.90      0.773940       0.147500       0.128200
0.95      0.886970       0.493040       0.484440
0.97      0.932180       0.695820       0.672130
0.98      0.954790       0.797220       0.776440
0.99      0.977390       0.898610       0.886130
1.00      1.000000       1.000000       1.000000

y                       M = 400
Level 2        Level 5        Level 7

0.00      0.000000       0.000000       0.000000
0.05     -0.012885      -0.0003243     -0.0002947
0.10     -0.025771      -0.0006451     -0.0005783
0.15     -0.038656      -0.0012314     -0.0010568
0.20     -0.051542      -0.0035617     -0.0033013
0.25     -0.064427      -0.0103240     -0.0109420
0.30     -0.077312      -0.0234520     -0.0238150
0.35     -0.090198      -0.0392330     -0.0398780
0.40     -0.103080      -0.0563100     -0.0573210
0.45     -0.115970      -0.0734660     -0.0746230
0.50     -0.128850      -0.0896280     -0.0906270
0.55     -0.015968      -0.102980      -0.1044000
0.60      0.096917      -0.113400      -0.1152300
0.65      0.209800      -0.120320      -0.1224000
0.70      0.322690      -0.123100      -0.1250500
0.75      0.435570      -0.120950      -0.1218600
0.80      0.548460      -0.108550      -0.1089800
0.85      0.661340      -0.041380      -0.0587660
0.90      0.774230       0.147230       0.1279300
0.95      0.887110       0.492860       0.4842500
0.97      0.932270       0.695720       0.6720100
0.98      0.954850       0.797140       0.7763500
0.99      0.977420       0.898570       0.8860800
1.00      1.000000       1.000000       1.000000

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

y                       M = 100
Level 2        Level 5        Level 7

0.00      0.00000        0.00000        0.00000
0.05     -0.01276       -0.0004015     -0.00040
0.10     -0.02552       -0.0007936     -0.00078
0.15     -0.03828       -0.0012214     -0.00120
0.20     -0.05104       -0.0017347     -0.00170
0.25     -0.06379       -0.0023872     -0.00238
0.30     -0.07655       -0.0034638     -0.00341
0.35     -0.08931       -0.0054841     -0.00524
0.40     -0.10207       -0.0100510     -0.00902
0.45     -0.11483       -0.0176790     -0.01558
0.50     -0.12759       -0.0277930     -0.02323
0.55     -0.01483       -0.0369840     -0.03037
0.60      0.09792       -0.0448550     -0.03652
0.65      0.21069       -0.0512220     -0.04149
0.70      0.32345       -0.0559550     -0.04517
0.75      0.43621       -0.0589840     -0.04743
0.80      0.54896       -0.0595660     -0.04808
0.85      0.66172       -0.0579270     -0.04688
0.90      0.77448       -0.0594230     -0.04293
0.95      0.88724        0.1492400      0.02744
0.97      0.93234        0.4895400      0.26696
0.98      0.95490        0.6597000      0.47095
0.99      0.97745        0.8298500      0.71853
1.00      1.000000       1.000000       1.000000

y                       M = 200
Level 2        Level 5        Level 7

0.00      0.000000       0.000000       0.0000000
0.05     -0.012632      -0.0002256     -0.0002200
0.10     -0.025265      -0.0004461     -0.0004284
0.15     -0.037897      -0.0006895     -0.0006578
0.20     -0.050530      -0.0009877     -0.0009424
0.25     -0.063162      -0.0013773     -0.0013403
0.30     -0.075795      -0.0020713     -0.0020025
0.35     -0.088427      -0.0037121     -0.0034697
0.40     -0.101060      -0.0088888     -0.0077991
0.45     -0.113690      -0.0178200     -0.0156190
0.50     -0.126320      -0.028356      -0.0237720
0.55     -0.013692      -0.037644      -0.0309950
0.60      0.098941      -0.045507      -0.0371250
0.65      0.211570      -0.051834      -0.0420640
0.70      0.324210      -0.056536      -0.0457100
0.75      0.436840      -0.059545      -0.0479510
0.80      0.549470      -0.060137      -0.0486150
0.85      0.662100      -0.058539      -0.0474550
0.90      0.774740      -0.060227      -0.0436020
0.95      0.887370       0.148400       0.0257970
0.97      0.932420       0.489040       0.2649700
0.98      0.954950       0.659360       0.4693400
0.99      0.977470       0.829680       0.7176000
1.00      1.000000       1.000000       1.000000

y                       M = 400
Level 2        Level 5        Level 7

0.00      0.0000000      0.0000000      0.0000000
0.05     -0.012567      -0.0001239     -0.0001164
0.10     -0.025134      -0.0002451     -0.0002267
0.15     -0.037701      -0.0003803     -0.0003491
0.20     -0.050268      -0.0005491     -0.0005036
0.25     -0.062835      -0.0007745     -0.0007269
0.30     -0.075402      -0.0012027     -0.0011303
0.35     -0.087969      -0.0026646     -0.0023786
0.40     -0.100540      -0.0083352     -0.0073288
0.45     -0.113100      -0.0181320     -0.0158290
0.50     -0.125670      -0.0287160     -0.0241450
0.55     -0.013103      -0.0380350     -0.0313710
0.60      0.099464      -0.0458770     -0.0374760
0.65      0.212030      -0.0521730     -0.0423870
0.70      0.324600      -0.0568550     -0.0460110
0.75      0.437170      -0.0598500     -0.0482400
0.80      0.549730      -0.0604460     -0.0489040
0.85      0.662300      -0.0588640     -0.0477610
0.90      0.774870      -0.0606400     -0.0439500
0.95      0.887430       0.1479700      0.0250270
0.97      0.932460       0.4887800      0.2640800
0.98      0.954970       0.6591900      0.4686200
0.99      0.977490       0.8295900      0.7171900
1.00      1.000000       1.000000       1.000000

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

Bn      # NL

10      63
20      72
50      78
100     83
500     88
```
COPYRIGHT 2017 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.