Finite Element Analysis on Bingham–Papanastasiou Viscoplastic Flow in a Channel with Circular/Square Obstacles: A Comparative Benchmarking

Abstract

A CFD (computational fluid dynamics) analysis was carried out for the Bingham viscoplastic fluid flow simulations around cylinders of circular and square shapes. The governing equations in space were discretized with the finite element approach via a weak formulation and utilizing Ladyzhenskaya–Babuška–Brezzi-stable pair $Q_2/P_1{}^{\mathrm{disc}}$ for approximation of the velocity and pressure profiles. The discrete non-linear system was linearized through Newton’s method, and a direct linear solver was iterated as an inner core solver. The study predicts the functional dependence and impact of Bingham number, $Bn$, on the drag coefficient and lift coefficient. The effect of the shape of an obstacle is also provided by providing comparative data for the hydrodynamic forces with the published results.

1. Introduction

It can be seen from the research over the last few decades that the emphasis in continuum and rheological mechanics has been completed with single-phase materials that do not require certain threshold stress for their actual deformation, such as polymer solutions and polymer melts. However, some fluid dynamic researchers, during their experiments, found substances that are recognized as viscoplastic materials, which require a certain threshold stress level for their deformation. Some of the fluids that exhibit viscoplastic features are slurries, pastes, chocolate, mayonnaise, margarine, suspension, etc. These materials are encountered in various industries, and food processes and have now become a topic of intensified interest for researchers due to their superb features and practical utilization in many scientific problems.

Shwedov [1] was the first to communicate the idea about the yield stress and its engineering reality. He also evaluated the flow behavior, quality, and composition of viscoplastic materials at different magnitudes of yield stress. Later, Bingham conducted vast experimentation in this direction and provided details about fluidity and plasticity [2]. After these initial studies, a great effort was made by researchers to investigate the viscoplastic materials in non-trivial flows, theoretically and experimentally, and several fluid models for the better physical interpretation of viscoplastic materials have been proposed. These models include the Bingham model [2], Herschel–Bulkley model [3], and Casson model [4], but the strongest one is the Bingham fluid model. Bingham did extraordinary work by proposing several paints that exhibit viscoplastic features and by elucidating the fluidity and plasticity of such materials. After that, seminal research work was conducted by Bird et al. [5,6] who provided a list of several other materials exhibiting these properties. Along those lines, several attempts [7,8,9,10] to modify the modeling of Bingham plastic fluids and the construction of numerical techniques for the solutions were made. Bercovier and Engelmann [7] identified the discontinuity in the Bingham model and gave a solution for eradication by linearizing the viscosity of the fluid. Papanastasiou [8] performed phenomenal work by adding an exponential term to the expression of yield stress. This addition to the model also helps in the description of yielded and non-yielded regions. A detailed review on the viscoplastic behavior of fluids is provided by Barnes [11], where he asserted that the real viscoplastic fluids behave more like the Papanastasiou regularized viscoplastic fluid rather than the ideal Bingham fluid. Some more studies for such yield-stress fluids are given in [12,13,14,15,16,17].

The flow around obstacles is the classical problem in fluid mechanics and is the topic of interest from the computational, experimental, and analytical point of view. In recent years, due to advancements in technology, the issues regarding the computational and time cost involved during the simulations around obstacles have been resolved. Researchers are now effectively working on the understanding of and physical interpretations of fluid flow around obstacles. In this regard, Schaefer et al. [18] performed exclusive work by investigating the flow features of Newtonian fluid around obstacles. Some results related to Newtonian fluid across the cylinder along with various imperative physical situations are available in [19,20,21,22]. The literature contains little work related to the flow of non-Newtonian fluids around obstacles. The theoretical work concerning this subject was studied by Adachi and Yoshioka [23]. Tokpavi et al. [24,25] investigated Bingham fluids around a circular cylinder by incorporating inertial effects. They experimentally analyzed the flow features of Bingham fluid and found excellent agreement between previously conducted theoretical findings with their experimental observations. Nirmalkar et al. [26] explicated forced thermal convection past the square cylinder by capitalizing the Bingham fluid. Moreover, they came up with the idea of ‘effective Reynolds number’, $R{e}^{*}=\frac{Re}{Bn+1}$, and ‘effective Bingham number’, $B{n}^{*}=\frac{Bn}{Bn+1}$, since the viscosity varies throughout the flow and effective viscosity might be more representative of the viscous stress. They provided correlations for the pressure and total drag as functions of the modified Reynolds number, $R{e}^{*}$. Both Mossaz et al. [27] and Syrakos et al. [28] computed the magnitude of drag and lift coefficients for Re = 45 and with variations in the Bingham number. Further, Syrakos et al. [29] tested the concept of effective Reynolds number for Bingham fluid, in other words, $R{e}^{*}=\frac{Re}{Bn+1}$ as discussed in Nirmalkar et al. [26], and showed that $R{e}^{*}$ is a better indicator of inertia in viscoplastic flows than the usual Reynolds number.

This study was aimed at investigating the features of Bingham–Papanastasiou fluids past circular and square cylinders in a channel. The fluid-governing rheological expressions were discretized by the finite element method. Code validation was also present to check the reliability of the computed results. To the author’s knowledge, the inertial flow of Bingham fluids has not been widely reported, and this work will serve as a reference study for upcoming research in this field. The remainder of this paper is structured as follows: The flow-governing equation and the expression of the features of the Bingham plastic model are presented in Section 2. The procedure of the implemented computational scheme is outlined and debated in Section 3. The grid independence and the code validation are reported in Section 4, as are numerical findings by way of sketches of the various involved quantities. Section 4 also contains plots for lift and drag coefficients. Conclusions are presented in Section 5.

2. Mathematical Formulation

The stationary incompressible flows are governed by the set of coupled partial differential equations given by:

$$\nabla .u=0,$$

(1)

$$\rho (u\cdot \nabla u)=-\nabla p+\nabla .\mathsf{\tau },$$

(2)

where the symbols have their usual meanings. Bingham [4] proposed a simplified rheological relationship for viscoplastic materials:

$$\{\begin{array}{ll}\dot{\mathsf{\gamma }}=0,& \tau \le {\tau }_y,\\ \mathsf{\tau }=\left(\frac{{\tau }_y}{\dot{\gamma }}+{\mu }_p\right)\dot{\mathsf{\gamma }},& \tau >{\tau }_y, \end{array}$$

(3)

where the ${\tau }_y$, ${\mu }_p$, $\mathsf{\tau }$, and $\dot{\mathsf{\gamma }}$, represents the yield stress, plastic viscosity, stress tensor, and the rate of strain tensor, respectively. The strain tensor is defined as

$$\dot{\mathsf{\gamma }}\equiv \nabla u+{\left(\nabla u\right)}^T,$$

(4)

here, $u$ denotes the velocity vector. The stress magnitude and strain rate are defined as

$$\tau \equiv {\left[\frac{1}{2}\mathsf{\tau }:\mathsf{\tau }\right]}^{\frac{1}{2}},\dot{\gamma }\equiv {\left[\frac{1}{2}\dot{\mathsf{\gamma }}:\dot{\mathsf{\gamma }}\right]}^{\frac{1}{2}}$$

(5)

A key finding for viscoplastic fluid flows is that the computational domain can be partitioned into three zones, the one where $\dot{\gamma }\ne 0$ represents the shear zone, while the plug zone remarks $\dot{\gamma }\equiv 0 \mathrm{and} u\ne 0$ and for the dead zone, we have $\dot{\gamma }\equiv 0 \mathrm{and} u=0$. Equation (3) inherits a discontinuity which was dealt with by Papanastasiou [8] by introducing an exponential function as

$$\mathsf{\tau }=\left[\frac{{\tau }_y}{\dot{\gamma }}\left\{1-\exp (-m\dot{\gamma })\right\}+{\mu }_p\right]\dot{\mathsf{\gamma }},$$

(6)

here, m is the parameter showing the stress growth. Owing the Equation (4), the viscosity can be written as

$$\eta =\frac{{\tau }_y}{\dot{\gamma }}\left\{1-\exp (-m\dot{\gamma })\right\}+{\mu }_p,$$

(7)

which is applicable in the whole flow domain.

Using non-dimensional variables $u^{*}, p^{*}, {\mathsf{\tau }}^{*}$ and by choosing $L_{ref}\hspace{0.17em}\hspace{0.17em}\mathrm{and} \hspace{0.17em}\hspace{0.17em}{U}_{ref}$ as reference length and velocity, respectively, so we obtained:

$$\nabla .u^{*}=0,$$

(8)

$$Re{u}^{*}.\nabla u^{*}=-\nabla p^{*}+\nabla .{\mathsf{\tau }}^{*},$$

(9)

in which

$${\mathsf{\tau }}^{*}=\left[\frac{Bn}{{\dot{\gamma }}^{*}}\left\{1-\exp (-M{\dot{\gamma }}^{*})\right\}+1\right]{\dot{\mathsf{\gamma }}}^{*}.$$

(10)

where $Re=\frac{\rho U_{ref}{L}_{ref}}{{\mu }_p}$ is the Reynolds number and $Bn=\frac{{\tau }_y{L}_{ref}}{{\mu }_p{U}_{ref}}$ is Bingham number. The stress s growth parameter m is now given by $M=\frac{m{U}_{ref}}{L_{ref}}$. The viscosity in non-dimensional form is

$${\eta }^{*}=\frac{Bn}{{\dot{\gamma }}^{*}}\left\{1-\exp (-M{\dot{\gamma }}^{*})\right\}+1$$

(11)

M is a non-dimensional analog of m. The value of M is taken as 500 in our simulations.

3. Numerical Procedure

Due to the lack of an analytical solution for the present configuration, we opted for a numerical procedure based on finite element discretization. The quadrilateral Stokes element $Q_2/P_1{}^{\mathrm{disc}}$ is applied from the finite element library. This element possesses shape functions that have nine degrees of freedom (DOFs) for velocity component locally and three local DOF for a piecewise linear discontinuous pressure approximation (see [30,31] for details). After linearizing the discrete non-linear algebraic system, the direct solver [32] is iterated in the inner loop. Grid refinement is achieved by forming four sub-elements from one parent element by connecting midpoints of opposite edges. In order to make a good comparison for the coefficients $C_D$ and $C_L$, the coarse grids are constructed in such a way that at a certain refinement level the DOFs are compatible (see

Table 1. Meshing scheme at different refinement levels.

Lev	Circular Cylinder		Square Cylinder
	#EL	#DOFs	#EL	#DOFs
1	130	1534	112	1344
2	520	5928	448	5152
3	2080	23,296	1792	20,160
4	8320	92,352	7168	79,744
5	33,280	367,744	28,672	317,184
6	133,120	1,467,648	114,688	1,265,152

).

4. Flow Configuration and Numerical Results

The simulation domains are schematically shown in Figure 1a,b. To be more specific, Figure 1a provides the domain as a rectangular channel $\Omega =[0,2.2]\times [0,0.41]$ with a circular cylinder and Figure 1b is specified with the square cylinder as an obstacle. The center of both the cylinders is set at C (0.2, 0.2). The structure of both meshes at the coarse level is shown in Figure 2. The velocity of the fluid at the upper and lower walls is taken as zero. A parabolic profile is exposed at inlet while the Neumann condition is carried at an outlet.

The reference velocity is taken $U_{\infty }=\frac{2}{3}{u}_{\max }$, where $u_{\max }$ is the maximum value of the parabolic input. The length of the diameter $D=0.1$ (for circular cylinder case) and side of the square of the same length as of D (for square cylinder case) is taken as reference length. The involved Bingham number is defined as $Bn=\frac{{\tau }_y.D}{{\mu }_p{U}_{\infty }}$ and $Re=\frac{\rho U_{\infty }D}{{\mu }_p}$ is Reynolds number. The dimensionless drag coefficient $C_D$ and lift coefficient $C_L$ are given as

$$C_D=\frac{2F_d}{\rho U_{{}_{\infty }}^{{}^2}D},$$

(12)

$$C_L=\frac{2F_l}{\rho U_{{}_{\infty }}^{{}^2}D}$$

(13)

4.1. Code Validation and Grid Independency

We performed some test simulations in order to validate the solver by setting $Bn=0$ and $Re=20$ and compared the important quantities with the results reported in [18]. The numerical values of the drag coefficient $(C_D)$ and lift coefficient $(C_L)$ are listed in

Table 2. Lift and drag coefficients for Newtonian case (Bn = 0).

Lev	${\mathit{C}}_{\mathit{L}}$	${\mathit{C}}_{\mathit{D}}$
1	4.811786 × 10−3	5.330670
2	9.006015 × 10−3	5.500019
3	1.042621 × 10−2	5.557633
4	1.042621 × 10−2	5.557633
5	1.060405 × 10−2	5.578171
6	1.061492 × 10−2	5.579198

and the good agreement is found with existing literature. The velocity distribution at $Re=20$ for $Bn=0$ is shown in Figure 3. The reference data is

$$\begin{array}{l}{C}_D=5.579535\\ \\ C_L=0.010618\end{array}$$

A comparison of values at level 5 and level 6 in Table 2 also confirms the grid independency of results.

4.2. Further Parametric Study: Graphical and Tabular Results

After having done the code validation, further parametric study for Bingham number is performed for the case of square cylindrical obstacle and quantities of interest are computed and a quantitative comparison is tabulated against circular obstacle case given in [33].

From

Table 3. Lift and drag coefficients towards Bn.

$\mathit{B}\mathit{n}$	Circular Obstacle		Square Obstacle
	${\mathit{C}}_{\mathit{D}}$	${\mathit{C}}_{\mathit{L}}$	${\mathit{C}}_{\mathit{D}}$	${\mathit{C}}_{\mathit{L}}$
1	6.237896	3.471506 × 10−2	7.802527	1.066414 × 10−1
5	8.956400	9.303962 × 10−2	11.23487	2.200490 × 10−1
10	12.76176	1.392485 × 10−1	16.00524	2.989600 × 10−1
15	16.76099	1.304087 × 10−1	20.95338	3.052176 × 10−1
20	20.81272	1.029614 × 10−1	25.93892	2.867756 × 10−1
25	24.86881	7.090231 × 10−2	30.91964	2.579676 × 10−1
30	28.91413	3.833286 × 10−2	35.88314	2.221996 × 10−1
35	32.94459	6.344714 × 10−3	40.82713	1.821432 × 10−1
40	36.95979	−2.481294 × 10−2	45.75258	1.430765 × 10−1
45	40.96059	−5.530762 × 10−2	50.66154	1.064771 × 10−1
50	44.94819	−8.510983 × 10−2	55.55585	7.492779 × 10−2

we can conclude that for various values of yield stress parameter represented by Bingham number, $C_D$ increases with Bn and $C_L$ show a decreasing trend. The is an increasing trend of drag as a function of Bn for both types of cylinders, however, the drag values are higher for the case of the square obstacle as the boundary of the square is not smooth as compared with the circular one.

Figure 4, Figure 5, Figure 6 and Figure 7 provide the variation in pressure, velocity, viscosity, and streamlines, respectively when Bn = 5, 10, 20 and 50. In detail, Figure 4 pressure profiles in the channel are presented, where it is observed that the pressure shows non-linear variations in the vicinity of the obstacle and it becomes linear downstream the obstacle. The velocity magnitude is displayed in Figure 5 for the selected values of Bn. The parabolic profile given at the inlet gets a speedy bifurcation at the cylinder and then it declines in the center of the channel due to plasticity effects increased by augmentation in the Bingham numbers.

From the viscosity plots given in Figure 6, the growth in the unyielded region is manifested. As we increase the yield stress parameter Bn, the unyielded region grows in the center of the channel downstream the channel and where the fluid moves like a rigid body showing the existence of plug zone. The flow pattern represented by streamlines is demonstrated in Figure 7 only for the selected values of Bn = 10 and Bn = 20, respectively, as the flow pattern remains almost unchanged.

For both the square and circular cylinders, the variation in u-velocity at x = 1 of channel is inspected towards different Bn and the observations are provided in Figure 8 and Figure 9. In detail, Figure 8 depicts the u-velocity when the square cylinder is placed as an obstacle towards ongoing Bingham fluid. It is noticed for the iterations $0\le Bn\le 50$ the u-velocity varies significantly. Figure 9 depicts the variation in u-velocity when the circular cylinder is placed as an obstacle towards Bingham fluid. From both figures, one can conclude that for the higher values of Bn, the peak of u-velocity in the center line declines significantly reflecting a plug zone in the flow and the boundary layer region becomes thinner at higher values of Bn.

Figure 10 and Figure 11 show the drag and lift coefficients at Re = 20 and for $0\le Bn\le 50$ both for circular and square obstacles. The trend of drag as a function of Bn is the same for both types of cylinders, however, the drag values are higher for the case of the square obstacle as the boundary is non-smooth compared with the circular one. The lift values show an increasing trend for $0\le Bn\le 10$ and then show a decreasing trend throughout for further augmentation in Bingham number for both types of obstacles.

5. Conclusions and Outlook

In this article, the interaction of Bingham viscoplastic fluid with both the square and circular cylinders as an obstacle is examined at a fixed Re = 20 and for $0\le Bn\le 50$. The finite element method is utilized to report flow field properties. It is observed that both the drag and lift coefficients subject to each obstacle are a function of Bingham number. The values of the drag coefficient increase significantly by increasing yield stress while the opposite trend is noticed for the case of lift coefficient confirming the resistive quality of yield stress. The drag approaches constant values at fixed Bingham number towards improved meshing which justifies the grid independency, hence, confirms the obtained results. Further, both the drag and lift coefficient are found higher for the square obstacle as compared to the circular one. In future, the idea of effective Reynolds number floated by Nirmalkar et al. [26] will be implemented and the study will be enhanced for a wide range of Bn and Re numbers.