Flow Characterization of Viscoelastic Fluids around Square Obstacle

This study focuses on the computational implementation of structured non-uniform finite volume method for the 2-D laminar flow of viscoelastic fluid past a square section of cylinder in a confined channel with a blockage ratio 1/4 for Re = 10-4, 5, 10 and 20. Oldroyd-B model (constant viscosity with elasticity) and the PTT model (shear-thinning with elasticity) are the constitutive models considered. In this study effects of the elasticity and inertia on the drag coefficients and stress fields around the square cylinder are obtained and discussed in detail. With an increase elasticity, drag coefficients get smaller due to stronger shear thinning effects for PTT fluid, however, the drag coefficients show slightly enhancement for the Oldroyd-B fluid.


Introduction
Numerical analysis of flow past the obstacles has received considerable attention for a long time [1][2][3][4][5]. Fluids around obstacles is encountered in many industrial applications such as shell and tube heat exchangers, coating processes, cooling towers, extruders and membrane processes. The bulk of the literature of fluid flow studies relate to flow over circular cylinders, followed by square, elliptic and rectangular cylinders or obstacles [6][7][8]. Indeed, even for the simplest shape of a circular cylinder which is free from geometrical singularities [7], the flow exhibits a rich variety of phenomena depending upon the nature of the mainstream flow (type of fluid Newtonian or non-Newtonian), blockage ratio of the cylinder (length to diameter ratio) and the characteristic Reynolds number of the flow. The behavior of a fluid flow past a square cylinder often has many complex phenomena such as flow separation, vortex shedding, recirculation length of wake flow (the flow path behind the square cylinder), distribution of the shear and normal profile around the solid surfaces of the square cylinder and, drag and lift force coefficients [9][10][11][12]. Majority of these studies in the literature deal with flow of a Newtonian fluid around a circular obstacle and to a smaller extend, around a square obstacle. For example Breuer et al. [9] examined laminar Newtonian flow around a square cylinder in a 2D channel using two different computational techniques, finite volume and lattice-Boltzmann automata. Their blockage ratio B, defined as the ratio between the obstacle dimension and channel height, was 1/8. They compared the results of the techniques in terms of velocity field, drag coefficient, Cd, recirculation length and Strouhal number. They observed an excellent agreement between the results of the techniques. A numerical study to investigate Newtonian flow past a square cylinder for Reynolds numbers Re ≤ 40 was conducted by Sen et al. [10] using a stabilized finite-element formulation with a non-uniform structured mesh. In order to mimic an unconfined flow, a small blockage ratio, B = 1/100, was employed. For comparison purposes they also used cylinders of elliptical and circular cross-sections. They investigated impact of Re on the flow separation angle and Cd. They found that, flow separation over the square cylinder occurs at a smaller Re compared to the other cylinders giving rise to the highest Cd among the investigated obstacles. The flow structures and wake flow characteristics, vortex shedding behavior of behind square cylinder at high Reynolds numbers were experimentally studied using particle image velocimetry (PIV) with different blockage ratios by Saha et al. [11]. They employed the two-dimensions flow past a stationary square cylinder at zero incidence for Reynolds number, Re  150 using a stabilized finite-element formulations. They also compared with a circular cylinder that the flow separated at a much lower Re from a square cylinder leading to the formation of a bigger wake. Consequently, at a given Re, the drag on a square cylinder was higher more than that on drag of a circular cylinder.
In spite of their important industrial implications, number of studies involving non-Newtonian fluids around enclosed obstacles is much smaller than those on Newtonian fluids in the literature. For example in their computational study Dhiman et al. [12] employed finite volume technique to investigate 2D flow of power-law fluids around a confined square cylinder. The fluids had index values between 0.5 ≤ n ≤ 2.0. Their results revealed that the effects of Re and B on the size of the recirculation zone and on Cd were stronger than that of the power-law index.
Rao et al. [13] extended the results on momentum and heat transfer characteristics to highly shear thinning fluids, especially n ≤ 0.5 at low Reynolds number. Fluid elements followed the contour of the square cylinder and flow remained attached to the surface. These works reveal that shear-thinning behavior increases both Cd and the rate of convective heat transfer from the square cylinder surface. Shear thinning behavior not only delays the formation of a visible wake but the resulting wake is also somewhat shorter than that of Newtonian fluid case. The shear thickening, on the other hand, has exactly the opposite influence on wake formation.
In the case of non-Newtonian fluid flow around confined obstacles, the literature is dominated by the studies with generalized Newtonian model to capture shear thinning or thickening effects. On the other hand, investigating the effects of viscoelasticity on the hydrodynamics of the flow around the obstacles has potentially crucial implications on many industrial applications. Therefore the authors believe that, the absence of viscoelastic flows around various obstacles in the literature merits a study on the hydrodynamics of such industrially important flows to reveal both microscopic and macroscopic flow characteristics. Hence, the objective of the present study is to investigate flow of a viscoelastic fluid, a linear PTT fluid and Oldroyd-B fluid, around a confined square obstacle computationally. Finite volume method is employed to solve coupled equations of continuity, motion and constitutive model along with appropriate boundary conditions. Effects of inertia in terms of Re, elasticity in terms of Weissenberg number, We, and constitutive equation parameters on the recirculation length, drag and on the flow field in terms of stress and velocity fields are examined and presented in detail.

Numerical Method
In this study isothermal flow of a viscoelastic fluid over a 2D confined square cylinder is considered. The flow system is schematically depicted in Fig. 1. The ratio between heights of the square and the channel, referred to as the blockage ratio, is 1/4 (b/H = 1/4). In the computations the upstream region length is set as 1/6 of the total channel length, L to ensure fully developed flow region. The ratio between channel length and height was set as L/H = 30.
The formulation of the flow begins by considering basic equations of fluid flow, i.e. continuity and momentum equations given below.
where u is the velocity, p is the pressure, η 0 is the total viscosity and τ represents the polymer or non-Newtonian contribution to the deviatoric stress tensor. The constant β is the ratio between the solvent viscosity and the total viscosity (β = η s / η 0 ). A viscoelastic constitutive model provides the additional relation needed to solve the conservation equations. In this study, Linear PTT model [14,16,17] and Oldroyd-B model [14,15,17] given by Eqs.  extensional parameter, (ε). The extensional parameter is considered ε ≥ 0, as a constant. The extensional parameter imposes an upper bound on the extensional viscosity which is inversely proportional to ε [16,17]. It is worth mentioning that ε = 0.25 in the PTT model correspond to the flow behavior of extremities for concentrated polymer melts polymer solutions [17]. The linear PTT model [14,16]: .
The Oldroyd-B model [15,16]: where ƞ p is the zero shear rate viscosity for viscoelastic flow contribution and τ is the extra stress tensor, λ 1 , is the relaxation time of flow. Total viscosity, ƞ 0 , ƞ 0 = ƞ P + ƞ S , is the sum of the viscoelastic flow contribution viscosity and solvent or Newtonian flow contribution viscosity parts. The Reynolds number (Re) and Weissenberg number (We) which is defined as the ratio of characteristic fluid relaxation time to characteristic time scale in the flow are given through: Where H and U are the characteristic length and velocity in the flow, respectively.
A finite volume method [18][19][20] with non-uniform staggered grid is used to obtain discrete form of the flow equations. Second order central difference scheme is used for the approximation of the diffusion terms in the momentum equations. Convective terms in the constitutive equations are approximated by at least second-order accurate, bounded and non-uniform version of CUBISTA [21] scheme. This scheme is preferred due to its documented advantages on the higher order schemes when viscoelastic fluids are considered [22].
The implementation of the CUBISTA scheme is carried out via deferred correction method that was proposed by Rubin and Khosla [23], Hayase et al. [24]. The method is widely used to ensure stability of the higher-order schemes for the evaluation of the variables at the faces of the control volumes. Mathematically this technique may be stated as; where ϕ f is flow property (u, v, p, τ xx , τ yy , τ xy ), the first term in Eq. (6) is the result from the low order (LO) scheme, and is used to evaluate the coefficient of the discretized equation. The other term is obtained at the previous iteration. Upwind Differencing Scheme (UDS) is used to handle the first term in Eq. (6). High order (HO) results are also obtained using Eq. (7). The purpose of the convection scheme use is then to specify the values of ϕ f at the face of the control volume, based on existing values at the neighbouring cell centres, ϕ c (ϕ P ). The detailed numerical computations were given in the study of Tezel (2016).
The final form of the two dimensional discretized governing equations over a control volume can be expressed symbolically as follows using finite volume method: Where ϕ i,j reprensents i th and j th component of ϕ in control volume. The coefficients are also approximated by the upwind differencing scheme. Gradients of velocities in the source terms of the constitutive equations are computed by central differences at interior points of the domain. However, at the flow boundaries they are calculated by approximating velocity profiles in the form of polynomial functions [24].
The SIMPLE [26] method is employed to handle the coupled system of the continuity, momentum and constitutive equations. The set of linearized algebraic equations are then solved by using the Thomas algorithm or the tridiagonal matrix algorithm (TDMA). The solution process is reiterated until the maximum relative change of flow variables, Ø i,j (u, v, p, τ xx , τ yy , τ xy ) is less than a prescribed tolerance or residual as: Use of correct boundary conditions especially for stress components is crucial to capture the physics of the flow. Boundary values of stress components are obtained by solving constitutive equations with the known velocity boundary conditions. The following inlet conditions are imposed for x and y-components of the velocity, u and v respectively. Corresponding values of stresses are given in Eq. (11).
No slip boundary conditions are assumed at the solid-fluid interfaces: Stress boundary conditions for solid channel walls can then be expressed by the following equations: At the channel exit Neumann boundary conditions are imposed for the flow variables.

Numerical Results and Discussion
A Reynolds number range 10 -4 ≤ Re ≤ 20 was investigated numerically, where Re is based on the cylinder diameter b and the maximum flow velocity u max of the parabolic inflow profile as seen in Fig. 1. The following section starts with a description of the different flow patterns observed with respect to Re and We. The subsequent sections present a detailed comparison of the computed results based on velocity, pressure and stress profiles at several positions in the flow field for both of the viscoelastic model.
Furthermore, the computations are analyzed and compared in terms of drag coefficient.

Streamlines around the obstacle
We begin presenting the computational results of the flow system by the inertial and elasticity effects on the recirculation patterns behind the obstacle depicted in Figs. 2(a) -(d).
The values of the constitutive model parameters used to obtain these results are ε = 0.25 and β = 0.2 . Dimensionless recirculation length that is also known as the wake region is defined by Breuer et al. [9] as the distance between the obstacle surface and reattachment point of streamlines to form the encapsulated region behind the obstacle. It should be noted that at a given Re upper limit of We is determined by the stability of the computations. The higher the value of Re, the lower We that can be attained for stable computations. For example when Re is set as 20, the maximum value of the attainable We is 3 and 6 for PTT and Oldroyd-B fluid respectively. Fig. 2 illustrate that increasing fluid elasticity or inertia leads to larger recirculation lengths and eventually formation of symmetric vortexes as depicted in Fig. 2(c) similar to the results reported by Breuer et al. [9]. Larger recirculation lengths and observed vortexes can be attributed to Hoop stresses getting stronger at increased fluid inertia and elasticity. Hu and Joseph [27] used UCM model, which has the same behavior as Oldroyd-B, and reported similar trend in the flow around circular cylinder. They also observed larger vortices with increasing elasticity.
The flow quantities such as stress and velocity fields are determined through the complex interactions between inertia, elasticity and shear thinning that gets stronger at elevated Re and We. In Figs. 3(a) -(c), viscoelastic fluid effects on the recirculation length (Lr) are compared at different Re. Oldroyd-B fluid with constant viscosity flow has larger wake region size or recirculation length depicted in Table 1. Huang and Feng [28] obtained the same wake lengths for Oldroyd-B flow around circular cylinder at Re = 10 when We is increased zero to one. Also, Lr is increased with Re numbers as in Newtonian flow.
Viscoelastic wake behind square obstacle is longer than the Newtonian wake (We = 0). Because, the wake formation has strong dependence on the structure of the flow and on the presence of vortices in the region of highest stress that resides downstream of the rear stagnation point of obstacle surface. Therefore, lacking shear thinning property, Oldroyd-B fluid leads to higher stresses and larger wake field around the obstacle compared to PTT and Newtonian flows. Oldroyd-B flow vortex centers also shift upward and downward direction with respect to PTT fluid due to expanding wake region in Figs. 3(b) and 3(c). This can be clearly seen in Table 2. The vortex pair size increases in wake region as Re number increases as in Fig. 3(c), when vortex intensities increase for both model listed in Table 2.
However, at Re = 10, the magnitude of vorticity becomes lower. For low Re number, inertial force gets lower and impact of the elasticity on the flow is comparatively smaller than at high Re number. So, at Re = 20, the elasticity of polymer molecule is higher for Oldroyd-B than PTT flow in Fig. 3(c). So, vortex intensity values in wake region (1.3079 < ψ < 1.3179) gets higher at x locations away from the obstacle at a constant y = 2.2006 for Oldroyd-B flow at Re = 20 as depicted in Table 2. Also, this vortex enhancement may be attributed to large and constant elongational viscosity of Oldroyd-B flow. Because, large elongational viscosity delays the acceleration of the fluid which results in the increase of vortex size. Fig. 4 presents velocity distribution of component, u the flow field for Re = 20 with different We numbers along the center line. Axial velocity, u, dominates the flow field as compared to vertical velocity, v. At x = 19 that is the position prior to the cylinder region, the velocity profiles differ from Poiseuille flow. Streamwise velocity, u, has a local minimum points for each We numbers. As the We number increases, the local minimum velocity decreases as in Fig. 4  velocity increases at elevated We. In the wake region, x = 22, u-velocity profiles become highly modified due to the presence of the obstacle. For both model, flow is nearly independent of We number at Re = 20. These results associated with the velocity field can be attributed to shear thinning property of the fluid that becomes more pronounced at elevated We values as shown in Fig. 4(a) for PTT fluid. On the other hand, u velocity for Oldroyd-B has more deformation due to more elastic behavior and the maximum velocity region shifts (see Fig. 4(b)). In other words, the elasticity changes dramatically the flow at high We numbers. The streamwise velocity at high We is slower to recover the undisturbed bulk velocity than low We flows especially in the wake region as seen in Figs. 2(a) -(d). Fig. 4(c) shows influence of We on profiles of the streamwise velocity component, u, along the centerline (y = 2) for PTT fluid. We also compare these results with the corresponding numerically computed velocity profiles for a Newtonian fluid under at Re = 20 (see Fig. 4(d)). Upstream of the cylinder the flow is essentially independent of We along the cylinder, in agreement with Figs. 4(c) -(d). For PTT flow, required length to achieve the fully developed velocity, ∂u / ∂x = 0, is nearly the same for all We numbers plotted in Figs. 4(a) and 4(c) while fully developed velocity in magnitude is lower than Newtonian flow. Axial velocity is shifted downstream and thus decreases relative to Newtonian flow with increasing We. For Oldroyd-B flow elasticity of the fluid also leads to increase in the required length for velocity recovery as shown in Fig. 4(d). Fig. 5 shows the pressure variation in the flow. For viscoelastic flow case, pressure gradients are generated by mainly elastic effects in near front and rear stagnation points of the cylinder. Viscoelastic pressure drop values are also higher than Newtonian flow case due to increasing effect of longitudinal flow as seen in Figs. 5(a) and 5(b). Fig. 5(a) depicts pressure profiles for PTT flow. As We increases, pressure value increases as shown in Fig. 5(a). This result can be attributed to the smaller extensional viscosities associated with higher We for a PTT fluid. Pressure drop also gets amplified as in Fig. 5(b) due to the breaking of the fore-aft symmetry of flow around obstacle (see Fig. 4(d)) at higher We for Oldroyd-B.

Normal and Shear stress profiles around the obstacle
In order to further examine the behavior of viscoelastic flows around the obstacle, we plot shear and normal stress profiles as in Fig. (6). τ xx and τ xy variation around the cylinder are given along x and y direction of the flow. Both shear and normal stresses are zero at the front and the rear stagnant points since the velocity gradient (∂u / ∂x = 0) is zero at the centerline as in Figs. 6(a) and 6(b). In the case of PTT fluid, at Re = 0 and 10, absolute value of normal and shear stresses decreases as We increases due to stress decaying of flow field as shown in Figs. 6 and 8. Higher the We is, stronger the shear thinning effects that lead to the lower stress values. As the inertial effects get amplified, τ xx and τ xy along the x direction becomes smaller at y = 0. Also, shear stress peak is observed at the front stagnation point due to sudden change in the boundary condition (at singularity point) for all We numbers as shown in Fig. 8(a). For the Oldrody-B fluid, absolute magnitude of stresses increases compared to PTT in the flow field as in Figs. 7 and 9. In Fig. 7(a) overshoot in τ xy profile occurs at We = 15. It may be resulted from u velocity gradient with respect to y direction at x = 19. On the other hand, as We gets higher, there is an increase in all stress gradients as shown in Fig. 7(b). The saturation of normal stresses at  the obstacle surface was also observed in creeping flow of Oldroyd-B fluid around cylinder [29,30]. At Re = 10, for PTT in the wake region, shear stress relaxation occurs in an oscillation manner as in Fig. 8(a). However, Oldroyd-B have more absolute shear stresses, stress relaxation is more quickly and suddenly owing to more flexible behavior of polymer chain, when the stress source removed (away from the cylinder), it is returned into undeformed case of the flow as in Fig. 9(a). Another important flow feature is the response of normal stress shear viscosity decreases along with the normal stress. For Oldroyd-B fluid, all elastic or normal stresses are nearly recovered as suggested by the plots in Fig. 7(b). However, as Re increases, elastic stresses reach a maximum value (We = 5) at the channel wall as in Fig. 9(b). Oldroyd-B flow relaxes quickly elastic stresses at the center of the obstacle compared as PTT flow.

Viscoelastic Drag Phenomena over the square cylinder
The corresponding drag coefficient can be splitted into three parts. These are drag coefficient due to normal stress contribution, drag coefficient due to shear stress contribution, drag coefficient due to pressure drop contribution [31]. All these contributions are expected to be functions of Re and We numbers and also coupled each other.
In the case of Oldroyd-B fluid, investigation of drag coefficients becomes more involved due to occurrence of higher stress gradients than those of PTT fluid. The drag increases with We number as tabulated in Table 3.
This increase in drag is associated with a remarkably long recirculation region as shown in Table 1. Another reason is the formation of Hoop stresses at high We flows in the wake region. For PTT fluid, this region is more stabilized owing to decaying of stress field. When We increases, shear and normal stresses contributions to drag coefficient are small due to constant viscosity of Oldroyd-B fluid at a given Re. At Re = 5, up to We = 4, there is an effective balance between normal stress and shear stress drag contribution leading to nearly the same drag coefficients. For We = 5 and 6, pressure drop contribution to drag coefficient decreases because pressure magnitude suddenly drops as in Fig. 7(c), while the normal stress contribution to drag coefficients gets amplified. Furthermore, Cd increases slightly with increased We at high Re. This behavior of drag can be explained by the secondary flow formation resulting from vortex pairs in the wake region.
At Re = 20, large deformation is observed in normal stress field especially at We = 3. Therefore, normal stress contribution to drag coefficient increases from 85 to 108, along with stronger pressure contribution as shown in Table 3. Also, the rapid change of pressure and the formation of normal stress boundary layer generates high extensions and shears around the obstacle. For We < 3, negative normal stress field around the obstacle has also decreasing effect in magnitude of normal stress contribution to drag  coefficient. Table 4 shows drag coefficients variation with polymeric viscosity ratio for PTT fluid. In our governing equations, polymer concentration effect is characterized by the polymer viscosity ratio, wr. Polymeric concentration affects both the viscosity and the relaxation time of polymer molecules in viscoelastic medium. We carry out a few tests for various wr for We = 1, 2 and 3 at fixed Re = 20 tabulated in Table 4. Although shear and normal stresses contribution to drag increase, the variation of Cd is nearly same for the range of We. However, pressure drop contribution has slightly increased with wr. The higher is the value of wr, the more drag enhances for We = 3 case due to increased elasticity effect. As demonstrated in Table 5, at higher wr, drag coefficient increases at a constant We. At high We flows, further increase of polymer concentration makes no significant contribution to drag force around the obstacle as shown in Table 5.

Conclusions
In this study numerical studies are carried out for steady laminar flows of viscoelastic fluids around the confined square cylinder to reach physical mechanism of flow behavior around the square obstacle. Oldroyd-B model (constant viscosity with elasticity) and the PTT model (shear thinning with elasticity) are used to capture viscoelasticity. Flows are simulated at various Reynolds and Weissenberg numbers by utilizing the finite volume method. The results obtained in this study allow one to draw the following conclusions: • Increasing fluid elasticity or inertia leads to larger recirculation lengths and eventually formation of symmetric vortexes as in Fig. 2(d) and Fig. 3 Normal stress relaxation occurs between obstacle surface and the channel wall. PTT fluid delays normal stress momentum transfer from obstacle surface to the channel wall. As the deformation rate is raised, shear viscosity decreases approaching to the channel wall and normal stress loss is observed. As shown in Figs. 7(b) and 9(b), for Oldroyd-B model, all elastic or normal stresses are nearly recovered and absorbed by the channel wall. It relaxes quickly elastic stresses compared as PTT model.