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

Authors: Asif Mehmood; Rashid Mahmood; Afraz Hussain Majeed; Farah Jabeen Awan
Journal: Modelling and Simulation in Engineering
 (2021)

---

## 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-stableP2−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.


---

## Body

## 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 20th 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(1)∇.u=0,(2)ρu⋅∇u=−∇p+∇.τ,where the all terms have their specific meanings. A simplified rheological relationship is proposed by Bingham [4] fluid for viscoplastic materials as
(3)γ̇=0,τ≤τy,τ=τyγ̇+μpγ̇,τ>τy,where τy, τ, μp, and γ̇, represent the yield stress, stress tensor, plastic viscosity, and shear rate, respectively. The strain tensor is defined as
(4)γ̇≡∇u+∇uT.Here, velocity vector is denoted asu. The shear rate and stress magnitude are defined as
(5)γ̇≡12γ̇:γ̇12,τ≡12τ:τ12.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
(6)τ=τyγ̇1−exp−mγ̇+μpγ̇.Here, the parameterm is showing the stress growth. Owing to Equation (4), the viscosity can be written as
(7)η=τyγ̇1−exp−mγ̇+μp,which is applicable in the whole flow domain. In order to introduce the nondimensional numbers for the simulation purposes, we choose LrefandUref as reference length and velocity, respectively, and nondimensional variables u∗,p∗,τ∗ and obtain the following formulation:
(8)∇.u∗=0,Reu∗.∇u∗=−∇p∗+∇.τ∗,in which
(9)τ∗=Bnγ̇∗1−exp−Mγ̇∗+1γ̇∗,where Re=ρUrefLref/μp and Bn=τyLref/μpUref are the Reynolds and Bingham number. The stress growth parameter m is now given by M=mUref/Lref. The viscosity in dimensionless form is
(10)η∗=Bnγ̇∗1−exp−Mγ̇∗+1.M is the nondimensional analogue of m.


## 3. Problem Setup

### 3.1. Flow Configuration

Two equally sized cylinders with diameterD=0.1m are placed in a rectangular channel at several relative positions. The dimension of the computational domain is 0,2.2m×0,0.41m. The circular obstacle C1 is fixed at 0.2,0.2, and C2 is placed at different locations to alter spacing between the cylinders Gp. The up and down walls of channel are set with no slip conditions. An inflow parabolic profile with umax=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).Figure 1
Schematic diagram of the problem.


### 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 pairP2−P1 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:
(11)ψn+1−ψnψn+1≤10−6.Here,n represents the iteration number, and ψ denotes a component of solution. Figure 2 shows the computational grids at refinement level 1.Figure 2
The computational grids at refinement level 1. (a)GP=0.1m. (b) GP=0.2m. (c) GP=0.3m.
(a)(b)(c)Following this strategy of refinement, Table1 shows the number of elements and degrees of freedom at different refinement levels.Table 1
Mesh statistics at different refinement levels.
Refinement levelsGp=0.1Gp=0.2Gp=0.3#ELDOF#ELDOF#ELDOF168618066861806694181821230298812703048126030333197044341990446420244515439828169403282444034824755958116555972116766018117456108021958410570192361080019581730544530342516244961250044472486314310865463371108996689751174029117698190551117958190941136230218349Domain discretization of channel with pair of circular obstacles in tandem arrangement (several gap spacing) at different refinement levels is provided in Table1. From the evaluated data about the degrees of freedom at Gp=0.1,0.2,0.3, it is concluded that for the high refinement levels, the degree of freedom is 190551 at Gp=0.1, 190941 at Gp=0.2, and 218349 at Gp=0.3. 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:(12)Fd=∫Sρν∂ut∂nny−pnxdS,Fl=−∫Sρν∂ut∂nnx−pnydS.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 coefficientCd and lift coefficient Cl are readily available as postprocessing by
(13)Cd=2FdρU2meanD,(14)Cl=2FlρU2meanD.Here,Umean is the average velocity, and D is the diameter of the obstacle.


## 3.1. Flow Configuration

Two equally sized cylinders with diameterD=0.1m are placed in a rectangular channel at several relative positions. The dimension of the computational domain is 0,2.2m×0,0.41m. The circular obstacle C1 is fixed at 0.2,0.2, and C2 is placed at different locations to alter spacing between the cylinders Gp. The up and down walls of channel are set with no slip conditions. An inflow parabolic profile with umax=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).Figure 1
Schematic diagram of the problem.


## 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 pairP2−P1 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:
(11)ψn+1−ψnψn+1≤10−6.Here,n represents the iteration number, and ψ denotes a component of solution. Figure 2 shows the computational grids at refinement level 1.Figure 2
The computational grids at refinement level 1. (a)GP=0.1m. (b) GP=0.2m. (c) GP=0.3m.
(a)(b)(c)Following this strategy of refinement, Table1 shows the number of elements and degrees of freedom at different refinement levels.Table 1
Mesh statistics at different refinement levels.
Refinement levelsGp=0.1Gp=0.2Gp=0.3#ELDOF#ELDOF#ELDOF168618066861806694181821230298812703048126030333197044341990446420244515439828169403282444034824755958116555972116766018117456108021958410570192361080019581730544530342516244961250044472486314310865463371108996689751174029117698190551117958190941136230218349Domain discretization of channel with pair of circular obstacles in tandem arrangement (several gap spacing) at different refinement levels is provided in Table1. From the evaluated data about the degrees of freedom at Gp=0.1,0.2,0.3, it is concluded that for the high refinement levels, the degree of freedom is 190551 at Gp=0.1, 190941 at Gp=0.2, and 218349 at Gp=0.3. 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:(12)Fd=∫Sρν∂ut∂nny−pnxdS,Fl=−∫Sρν∂ut∂nnx−pnydS.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 coefficientCd and lift coefficient Cl are readily available as postprocessing by
(13)Cd=2FdρU2meanD,(14)Cl=2FlρU2meanD.Here,Umean 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 spacingGp=0.1,0.2,and0.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 (Gp=0.1,0.2,and0.3) and for different values of Bn. For all Gp 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.Table 2
Influence ofBnand separation distance on drag coefficientCd.
BnCd for left cylinderCd for right cylinderGp=0.1Gp=0.2Gp=0.3Gp=0.1Gp=0.2Gp=0.302.6802202.7207382.7569181.2399241.7041892.00348913.1320793.2238813.2783681.7319342.2986332.67233255.1372765.4542555.5966193.8059264.8192915.431213107.6781708.3100398.6145426.3679757.9450508.7469711510.2562811.2111711.672508.85809711.0167611.974152012.8276514.1396014.7465611.3123714.0507715.167582515.3740417.0689017.8216813.7557617.0554618.342793017.8877619.9920520.8977216.1981820.0365721.510333520.3752022.9018723.9691318.6327022.9990724.669394022.8408425.7940727.0344721.0672325.9457127.819744525.2889528.6716830.0946923.4977128.8798530.963365027.7236631.5386833.1486025.9278431.8064734.10000Table 3
Influence ofBnand separation distance on lift coefficientCl.
BnLeft cylinder liftRight cylinder liftGp=0.1Gp=0.2Gp=0.3Gp=0.1Gp=0.2Gp=0.300.0012010.0071930.0063010.0341270.0295420.01796610.0792630.0186590.0200470.0595780.0497590.0315850.2822850.0676470.0633120.1297550.0619620.027698100.3720300.1573880.1424200.1540890.007122-0.04540150.4437840.2404560.2151860.155495-0.05300-0.13735200.5059670.3137980.2844770.146449-0.10727-0.23050250.5673200.3759150.3498710.136279-0.15612-0.31986300.6263030.4354850.4120740.125025-0.20105-0.40430350.6895890.4877220.4711270.117755-0.24332-0.48369400.7472400.5329270.5262680.101797-0.28519-0.55947450.7979910.5767250.5797890.080313-0.32648-0.63356500.8472870.6185340.6337240.058068-0.36710-0.70511Figure 3
Velocity profiles for variousBnwith Gp=0.1 andRe=20. (a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)Figure 4
Velocity profiles for variousBn with Gp=0.2 andRe=20. (a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)Figure 5
Velocity profiles for variousBn with Gp=0.3 andRe=20. (a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)Figures6, 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).Figure 6
Pressure profiles for variousBnwith Gp=0.1 andRe=20.(a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)Figure 7
Pressure profiles for variousBn with Gp=0.2 andRe=20. (a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)Figure 8
Pressure profiles for variousBn with Gp=0.3 andRe=20. (a) Bn=0. (b) Bn=10. (c) Bn=50.
(a)(b)(c)The impact of pertinent parameters on viscosity is plotted in Figures9, 10, and 11. It is observed that viscosity increases 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.Figure 9
Viscosity profiles for variousBn with Gp=0.1 and Re=20. (a) Bn=10. (b) Bn=50.
(a)(b)Figure 10
Viscosity profiles for variousBn with Gp=0.2 andRe=20. (a) Bn=10. (b) Bn=50.
(a)(b)Figure 11
Viscosity profiles for variousBn with Gp=0.3 andRe=20. (a) Bn=10. (b) Bn=50.
(a)(b)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,and0.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).Figure 12
Impact onu velocity for different Bn. (a) Gp=0.1m,x=0.3m. (b) Gp=0.2m,x=0.35m. (c) Gp=0.3m,x=0.4m.
(a)(b)(c)The impact of spacing on vertical velocity is negligible; however, for Newtonian caseBn=0, it gets some peak as described in Figure 13.Figure 13
Impact onv velocity for different Bn. (a) Gp=0.1m,x=0.3m. (b) Gp=0.2m,x=0.35m. (c) Gp=0.3m,x=0.4m.
(a)(b)(c)For a better analysis, the data presented in Tables2 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.Figure 14
Influence ofBnand separation distance on drag and lift coefficient. (a) Drag for left cylinder. (b) Lift for left cylinder. (c) Drag for right cylinder. (d) Lift for right cylinder.
(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 ofBn ranging from 0 to 50. The software package SPSS 23.0 is used for this purpose.In Table4, 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 Gp=0.1, the drag and lifts are almost independent. But as the gap increases, the correlation becomes strong.Table 4
Correlation between drag and lift.
Left cylinder (upstream)Right cylinder (downstream)Gp=0.10.969-0.054Gp=0.20.993-0.992Gp=0.30.998-0.995


## 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 numberBn 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.



---

*Source: 5583110-2021-04-28.xml*