#### Abstract

The numerical modeling and simulation for the stationary Bingham fluid flow around two confined circular cylinders with various gap ratios are studied. The singularity in the model’s apparent viscosity is dealt by Papanastasiou’s regularization. The model equations are discretized by adopting the methodology based on finite element method (FEM) by choosing a mixed higher order LBB-stable finite element pair. The direct solver PARADISO has been utilized to solve the linearized system of equations. Hydrodynamic forces represented by drag and lift coefficients are computed, and a correlation coefficient is calculated for the gap ratios and for several values of the Bingham number . Line graphs for horizontal and vertical velocities are drawn. Moreover, velocity and pressure profiles are plotted for pertinent values of the parameters. Plug and shear zones are revealed via velocity snapshots in the domain. Pressure is nonlinear in the vicinity of the obstacles and becomes linear downstream in the cylinders as expected in channel flows.

#### 1. Introduction

The theoretical hydrodynamics prospered in nineteenth century but was a bit distant from hydraulic science. The introduction of boundary layer concept in the 20^{th} century reduced this gap to zero. The fluid mechanics suffering due to the complexity of fluid models in analytical solutions have crossed this barrier due to numerical procedures and invention of computer and technological advancement in the recent past. An important question is what should be the next step? The intellect is the one who looks beyond the senses. Non-Newtonian fluids, high Reynolds numbers, and turbulent flows [1–3] captured the attention of mathematicians and engineers in the past few decades. The mysteries of turbulence are still being revealed in order to understand the flow phenomenon. The notable contributions of Euler, Bernoulli, and d’ Alembert have been explained in detail by Darrigol and Frisch [3].

Over the past two decades, the flows of both Newtonian and non-Newtonian fluids have been investigated thoroughly [4–10]. For various values of Reynolds numbers, a variety of flow patterns can be observed even for Newtonian fluids [4, 5]. The impact of fluid interaction with bluff bodies is very complex. So far, less information is available in literature about such complexities; a few studies can be found in [4–7].

The hydrodynamic forces such as drag and lift are the key factors to be investigated in studying solid-liquid interaction. While investigating the flow characteristics when more than one obstacle is there, the arrangement and placement of obstacles are of particular interest, as it has meaningful impact on hydrodynamic forces. Two cylinders with different arrangements is the simplest of such cases [8]. The shear and dead zones are dependent not only on the placement of obstacles but on their sizes as well. Such zones are important to be investigated as they have practical significance in many fields such as placement of means for bridges, fluid flow in channels with submerged obstacles, and of course in food industry while dealing with mixing of different products. Tandem and side-by-side are the two basic configurations of obstacle arrangements, of which tandem arrangement captured more attention of researchers. The drag and lift have different impacts on up- and downstream obstacles in tandem arrangements. The gap spacing of obstacles is also of keen interest to control the impact of fluid forces especially drag and lift.

While studying the rheological aspects of nonlinear viscosity fluids, there are two types of fluid models to deal with: shear rate-dependent models [11–17] and yield stress models [18, 19]. The materials like toothpaste and tomato ketchup are yield stress oriented. The rheology of such is investigated by many researchers [18–24]. There are some recent discussion and hypothesis about the existence of yield stress. One statement is that the fluid behaves like solid below a certain threshold value of yield stress (of course fluid flows beyond that threshold value). The second statement summarizes that, all solids and fluids can flow over a long period of time [25, 26]. The first statement looks more practical as described in [27]. Beris et al. [28] were the first to investigate yield stress of the Bingham fluids using FEM. Adachi and Yoshioka [29] also established some results of viscous fluids past a cylinder. Mitsoulis [30] also used FEM to investigate the hydrodynamic forces of the Bingham fluid flowing around the cylinder. Some studies are available in literature which confined to Newtonian fluid past a couple of cylinders; one of these is proposed by Vakil and Green [31] and produced some results for Reynolds number between 1 and 20. Jossic and Magnin produced some results [32] of two cylinders by fixing the range of the Bingham number up to 40. Recently, Koblitz et al. [33] investigated the hydrodynamic forces for very high values of Bn. Calder and Yezzi [34] considered a rate of convergence and application for obstacle problems. In some studies mentioned therein [35–38], authors have numerically investigated the flows of incompressible non-Newtonian fluids in various computational domains. Ameur and Menni [36] performed simulations at high Reynolds numbers giving rise to a turbulent flow in a backward facing step configuration, and a comparison with the experimental studies is performed. Houssem and Mohamed [39–41] analyzed a mixed convective thermal flow over a single cylinder and two cylinders. They have measured the effects of wall proximity of obstacles and gap ratios for varying Reynolds and Richardson number. In addition, they have also considered the effects of the Reynolds and Richardson numbers on fluid forces and the Nusselt number.

We organize the rest of the manuscript as follows. “Mathematical Modeling” is concerned with the mathematical modeling of flow equations for a regularized Bingham fluid. The nondimensional reformulation of model equations is also presented in the same section. Problem setup and numerical approach including element type and solvers are described in “Problem Setup.” “Results and Discussions” is devoted to display the results of numerical simulations for the involved parameters and the analysis of gap ratio on hydrodynamic forces. The correlation coefficient between the forces versus gap between the cylinders is computed. The conclusion of the present study is revealed in “Conclusions.”

#### 2. Mathematical Modeling

The set of coupled partial differential equations for stationary incompressible flows in dimensional form are defined generally as below where the all terms have their specific meanings. A simplified rheological relationship is proposed by Bingham [4] fluid for viscoplastic materials as where , , , and , represent the yield stress, stress tensor, plastic viscosity, and shear rate, respectively. The strain tensor is defined as

Here, velocity vector is denoted as . The shear rate and stress magnitude are defined as

Unyielded regions produced by (3) do not create any problem as long as , but apparent viscosity gives rise to a singularity in the limit . Computational schemes such as finite element method implemented in our study do not allow such singularities. To circumvent this issue, some regularization schemes are present in the literature. The purpose of regularization schemes is to replace discontinuous apparent viscosity for with such an expression that remains bounded for arbitrarily small and approximating the rheological behavior at the same time. In this direction, we removed the singularity with the help of the regularization introduced by Papanastasiou [20] as

Here, the parameter is showing the stress growth. Owing to Equation (4), the viscosity can be written as which is applicable in the whole flow domain. In order to introduce the nondimensional numbers for the simulation purposes, we choose as reference length and velocity, respectively, and nondimensional variables and obtain the following formulation: in which where and are the Reynolds and Bingham number. The stress growth parameter is now given by . The viscosity in dimensionless form is

is the nondimensional analogue of .

#### 3. Problem Setup

##### 3.1. Flow Configuration

Two equally sized cylinders with diameter are placed in a rectangular channel at several relative positions. The dimension of the computational domain is . The circular obstacle is fixed at , and is placed at different locations to alter spacing between the cylinders . The up and down walls of channel are set with no slip conditions. An inflow parabolic profile with is exposed to the inlet of the channel, and a do-nothing boundary condition at the outlet is chosen (see Figure 1).

##### 3.2. Numerical Approach and Solvers

A wide range of flow problems can be described with the incompressible Navier-Stokes equations provided in (2). These equations describe real-life processes and help to understand nature. However, there are many difficulties associated with the numerical treatment of the system given in (1) and (2). Firstly, the nonlinearity of the convective term gives rise to numerical instabilities at larger values of Reynolds numbers. The second problem is caused by the incompressibility constraint that gives rise to checkerboard pressure modes resulting in pressure oscillations. The third problem is associated with the nonlinear algebraic system that results after the discretization of these equations. The size of the global matrix is very large and in general has bad condition number, so to tackle these issues, efficient iterative solvers with preconditioning techniques [42] are required.

Keeping in view all these difficulties, the methodology was adopted in the finite element method. The LBB-stable FEM pair is chosen to approximate the discrete velocity and pressure. A hybrid coarse mesh is first created and then is successively refined to meet the grid independent results. For the solution of nonlinear algebraic systems, Newton’s method is selected to linearize the system, and a direct linear solver PARADISO with special reordering of unknowns is applied for the linearized system. To stop the nonlinear iterative process, we adopt the following criteria:

Here, represents the iteration number, and denotes a component of solution. Figure 2 shows the computational grids at refinement level 1.

**(a)**

**(b)**

**(c)**

Following this strategy of refinement, Table 1 shows the number of elements and degrees of freedom at different refinement levels.

Domain discretization of channel with pair of circular obstacles in tandem arrangement (several gap spacing) at different refinement levels is provided in Table 1. From the evaluated data about the degrees of freedom at , it is concluded that for the high refinement levels, the degree of freedom is 190551 at 190941 at , and 218349 at . Table 1 also provides the number of domain elements as well as boundary elements for all the generated meshes.

##### 3.3. Quantities of Interest

The important quantities are values of drag and lift which represent the forces in the horizontal and vertical directions, respectively, defined by the surface integration over the cylinders:

These drag and lift forces combine pressure and viscous forces as evident from the above surface integrals. Once the hydrodynamic forces are available, then their nondimensional analogue namely the drag coefficient and lift coefficient are readily available as postprocessing by

Here, is the average velocity, and *D* is the diameter of the obstacle.

#### 4. Results and Discussions

The incompressible fluid is considered in the present study to explore some new dynamics of hydrodynamic forces. For this purpose, a pair of obstacles of circular shapes is placed in a channel to observe the impact of the Bingham fluid flow. The upstream cylinder is fixed at the position (0.2, 0.2), and the position of downstream cylinder is changed at three different locations to establish the gap spacing . For the fixed value of Reynolds number , the flow of the Bingham fluid is observed. Starting from (the Newtonian case) up to , the drag and lifts are calculated at both the cylinders (see Tables 2 and 3). Fluid velocity, pressure, and viscosity are observed for gap spacing () and for different values of . For all values, it is observed that velocity decreases with the increase in (see Figures 3, 4, and 5(a)–5(c)) and plug zone extends from the center of channel towards the solid walls, and shear zone is restricted to the neighborhood of obstacles.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

Figures 6, 7, and 8(a)–8(c) show the pressure profiles representing the stagnation zones. Stagnation point appears as the fluid hits upstream cylinder. There is lesser impact of pressure on downstream cylinder for (Figure 6(a)). However, it increases with the increase in (Figure 6(c)), and it is also observed that the pressure increases on downstream cylinder with the increase of gap spacing (Figure 8).

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

The impact of pertinent parameters on viscosity is plotted in Figures 9, 10, and 11. It is observed that viscosity increases with the increase in . Higher values of viscosity observed even between the two cylinders as increases. Moreover, small islands of high viscosity appear around the cylinders for all gap spacing and their size grows with an increase in .

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

The steady flow in two-dimensional space is considered for analysis. Cut lines are marked exactly at the middle of two cylinders () to observe the velocity impact in the center of the two cylinders, and and components are considered for minute analysis. Plots show the dependence of stream lines on , , and. For , the visual impact of velocity is observed for the component in Figure 12(a). For (non-Newtonian case), the fluidity is much higher than the rest of the . The gradual increase in nonlinearity due to increase in can also be observed. The increase in gap ratio also has an impact on flow pattern in the middle of cylinders that can be seen in Figures 12(b) and 12(c).

**(a)**

**(b)**

**(c)**

The impact of spacing on vertical velocity is negligible; however, for Newtonian case , it gets some peak as described in Figure 13.

**(a)**

**(b)**

**(c)**

For a better analysis, the data presented in Tables 2 and 3 have been plotted in Figure 14 to see the trend of gap ratio versus . A linear increasing profile for drag on both left and right cylinders is observed. Furthermore, the drag coefficient is higher for left cylinder as compared with the right one since the fluid forces are dominant on the upstream cylinder. The lift coefficient varies in a nonlinear fashion for all gap ratios; however, for left cylinder, it increases with an augmentation in whereas for right cylinder, the lift coefficient drops for after a certain threshold of .

**(a)**

**(b)**

**(c)**

**(d)**

A correlation analysis has also been done to check the strength of relationship between drag and lift coefficients at various gap ratios and for all the values of ranging from 0 to 50. The software package SPSS 23.0 is used for this purpose.

In Table 4, the correlation between drag and lift is evaluated for for left and right cylinders, respectively, with different spacing. For left cylinder, the drag and lift are positively correlated, and the relationship is very strong. It gets even stronger as the gap increases. For right cylinder, the drag and lift are inversely related for all gap spacing. For downstream cylinder at , the drag and lifts are almost independent. But as the gap increases, the correlation becomes strong.

#### 5. Conclusions

We have simulated the flow of the Bingham fluid around obstacles using the Papanastasiou regularization and the finite element method. The effects of the Bingham number and gap spacing on the drag and lift coefficient of the cylinders have been investigated in detail. It has been observed that the separation plays a key role in determining the drag and lift coefficients, and it influences the shape of unyielded zones near cylinders. An increment in the Bingham number results in the enhancement of plug zone towards the walls of the channel downstream the obstacles. Increasing the separation between the cylinders gives rise to more stagnation pressure on the downstream cylinder. The drag and lift coefficients are positively correlated for the upstream cylinder for all gap spacing while the correlation becomes negative for the downstream cylinder.

#### Nomenclature

: | Yield stress |

: | Plastic viscosity |

: | Rate of strain tensor |

: | Shear rate |

: | Stress tensor |

: | Gap ratio |

: | Dimensional velocity vector |

: | Dimensionless velocity vector |

: | Reference velocity |

: | Dimensional pressure |

: | Dimensionless pressure |

: | Dimensional stress growth parameter |

:_{M} | Dimensionless stress growth parameter |

: | Reynolds number |

: | Bingham number |

#EL: | Number of elements |

DOF: | Number of degrees of freedom |

: | Drag coefficient |

: | Lift coefficient |

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors are very thankful for the anonymous referees whose suggestions and comments helped to improve the quality of the manuscript.