Flow of the Bingham-Papanastasiou Regularized Material in a Channel in the Presence of Obstacles: Correlation between Hydrodynamic Forces and Spacing of Obstacles

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 P2 − P1 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 ð0:1 ≤Gp ≤ 0:3Þ and for several values of the Bingham number ð0 ≤ Bn ≤ 50Þ. 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.


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][2][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][5][6][7][8][9][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][5][6][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 Figure 1: Schematic diagram of the problem.   1  686  1806  686  1806  694  1818  2  1230  2988  1270  3048  1260  3033  3  1970  4434  1990  4464  2024  4515  4  3982  8169  4032  8244  4034  8247  5  5958  11655  5972  11676  6018  11745  6  10802  19584  10570  19236  10800  19581  7  30544  53034  25162  44961  25004  44724  8  63143  108654  63371  108996  68975  117402  9  117698  190551  117958  190941 136230 218349    While studying the rheological aspects of nonlinear viscosity fluids, there are two types of fluid models to deal with: shear rate-dependent models [11][12][13][14][15][16][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][19][20][21][22][23][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 [   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."

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 τ y , τ, μ p , and _ γ, represent the yield stress, stress tensor, plastic viscosity, and shear rate, respectively. The strain tensor is defined as

Modelling and Simulation in Engineering
Here, velocity vector is denoted as u. The shear rate and stress magnitude are defined as Unyielded regions produced by (3) do not create any problem as long as _ γ > 0, but apparent viscosity gives rise to a singularity in the limit _ γ ⟶ 0. 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 _ γ ⟶ 0 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 m 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 L ref and U ref as reference length and velocity, respectively, and nondimensional variables u * , p * , τ * and obtain the following formulation: in which where M is the nondimensional analogue of m.   Modelling and Simulation in Engineering several relative positions. The dimension of the computational domain is ½0, 2:2m × ½0, 0:41m. The circular obstacle C 1 is fixed at ð0:2, 0:2Þ, and C 2 is placed at different locations to alter spacing between the cylinders G p . The up and down walls of channel are set with no slip conditions. An inflow parabolic profile with u max = 0:3m/s is exposed to the inlet of the channel, and a do-nothing boundary condition at the outlet is chosen (see Figure 1).

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 LBBstable FEM pair P 2 − P 1 is chosen to approximate the dis-crete 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, n represents the iteration number, and ψ denotes a component of solution. Figure 2 shows the computational grids at refinement level 1.
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 G p = 0:1, 0:2, 0:3, it is concluded that for the high refinement levels, the degree of freedom is 190551 at G p = 0:1, 190941 at G p = 0:2, and 218349 at G p = 0:3.
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 C d and lift coefficient C l are readily available as postprocessing by Here, U mean is the average velocity, and D is the diameter of the obstacle.

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 G p = 0:1, 0:2, and 0:3. For the fixed value of Reynolds number Re = 20, the flow of the Bingham fluid is observed. Starting from Bn = 0 (the Newtonian case) up to Bn = 50, 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 (G p = 0:1, 0:2, and 0:3) and for different values of Bn. For all G p values, it is observed that velocity decreases with the increase in Bn (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. 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 Bn = 0 (Figure 6(a)). However, it increases with the increase in Bn (Figure 6(c)), and it is also observed that the pressure increases on downstream cylinder with the increase of gap spacing (Figure 8).
The impact of pertinent parameters on viscosity is plotted in Figures 9, 10,       with the increase in Bn. Higher values of viscosity observed even between the two cylinders as Bn increases. Moreover, small islands of high viscosity appear around the cylinders for all gap spacing and their size grows with an increase in Bn.
The steady flow in two-dimensional space is considered for analysis. Cut lines are marked exactly at the middle of two cylinders (x = 0:3, 0:35, and 0:4) to observe the velocity impact in the center of the two cylinders, and u and v components are considered for minute analysis. Plots show the dependence of stream lines on Re, Bn, andGp. For Gp = 0:1, the visual impact of velocity is observed for the u component in Figure 12(a). For Bn = 0 (non-Newtonian case), the fluidity is much higher than the rest of the Bn. The gradual increase in nonlinearity due to increase in Bn 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).
The impact of spacing on vertical velocity is negligible; however, for Newtonian case Bn = 0, it gets some peak as described in Figure 13.
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 Bn. 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 Bn whereas for right cylinder, the lift coefficient drops for after a certain threshold of Bn.
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 Bn 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 0 ≤ Bn ≤ 50 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 G p = 0:1, the drag and lifts are almost independent. But as the gap increases, the correlation becomes strong.

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 Bn 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 τ y : Yield stress μ p : Plastic viscosity _ γ: Rate of strain tensor _ γ: Shear rate τ: Stress tensor G p : Gap ratio u: Dimensional velocity vector u * : Dimensionless velocity vector U mean : Reference velocity p: Dimensional pressure p * : Dimensionless pressure m: Dimensional stress growth parameter M : Dimensionless stress growth parameter Re: Reynolds number Bn: Bingham number #EL: Number of elements DOF: Number of degrees of freedom C d : Drag coefficient C l : 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.