Response of Viscoelastic Turbulent Pipeflow Past Square Bar Roughness: The Effect on Mean Flow

Abstract

The influence of viscoelastic polymer additives on response and recovery of turbulent pipeflow over square bar roughness elements was examined using Direct Numerical Simulations at a Reynolds number of $5\times {10}^3$. Two different bar heights for the square bar roughness elements were examined, $h/D=0.05$ and $0.1$. A Finitely Extensible Non-linear Elastic-Peterlin (FENE-P) rheological model was employed for modeling viscoelastic fluid features. The rheological parameters for the simulation corresponded to a high concentration polymer of 160 ppm. Recirculation regions formed behind the bar elements by the viscoelastic fluid were shorter than those associated with Newtonian fluid, which was attributed to mixed effects of viscous and elastic forces due to the added polymers. The recovery of the mean viscoelastic flow was faster. The pressure losses on the surface of the roughness were larger compared to the Newtonian fluid, and the overall contribution to local drag was reduced due to viscoelastic effects.

1. Introduction

Response and recovery of turbulent pipeflow to abrupt surface variations, i.e., roughness elements, is a fundamental study of non-equilibrium flow behaviour with various industrial and engineering applications. For example, extraction and transport of various heavy crude and petroleum products are routinely carried out using pipelines [1]. The flow in most crude transport pipelines experiences small eddies and near-wall turbulence, which increases frictional drag. This incurs heavy energy losses in the transport of highly viscous fluids like crude or heavy-oil bitumen [2]. Thus, techniques for frictional drag reduction are of significant practical interest. It has been established that the addition of high molecular-weight non-Newtonian polymers in turbulent pipeflow produced upward of $80\%$ drag reduction, and it significantly reduced turbulent frictional drag on the pipe walls [3,4]. This phenomenon was first reported by Toms (1948) [5], and such polymer additives were named Drag-Reducing Agents (DRAs). The addition of such polymers led to the fluid experiencing elastic viscosity as well as Newtonian viscosity; hence, they are referred to as viscoelastic fluids [6]. Polymer additives for purposes of drag reduction have been used successfully in many applications including the transport of crude in pipelines [7], energy saving transport medium in residential heating devices [8], as well as sewer [9], firefighting systems [10], and hydrodynamic systems [11].

Transport of energy-related fluids through pipes and channels involves a high Reynolds number, and in most practical cases, includes turbulent flows over rough surfaces. Abrupt surface variations at high Reynolds numbers constitute a class of non-equilibrium or perturbed flows, which have been the focus of extensive research using Newtonian fluids. The behaviour of such flows is typically complex because the perturbations or step changes cause contraction in the flow, particularly at a high Reynolds number. There are several studies that investigate different aspects of Newtonian fluid flow past wall-mounted obstacles or perturbations, and non-equilibrium flow conditions. Smits et al. (1979) [12] studied the effects of abrupt surface variations on flow recovery and observed a second-order response, which led to long-lasting changes in the turbulent structures in the flow. This was followed by the study of Durst and Wang (1989) [13], who investigated, experimentally and numerically, the flow response to an axisymmetric perturbation. They found an overshoot of flow response in the vicinity of the perturbation due to the mixed effects of sudden contraction and expansion, as opposed to sudden expansion in the flow over a backward facing step. Dimaczek et al. (1989) [14] observed similar effects for wall-mounted axisymmetric perturbations over a wide range of Reynolds numbers. Studies of Jiménez (2004) [15] and Smits et al. (2019) [16] focused on turbulent flow over small step conditions or roughness elements, defined by the ratio of the lateral width of the object (h) with respect to the boundary layer width (${\delta }_o$), given as $h/{\delta }_o\ll 1$. These small step conditions were unique in terms of their effects on the overall flow response and downstream recovery. Smits et al. (2019) [16] studied the flow response and recovery over two roughness heights of $h/D=0.05$ and $0.1$ at Reynolds number of $1.56\times {10}^5$ using Particle Image Velocotimetry (PIV). This study identified that the sudden contraction and expansion of the incoming flow resulted in longer reattachment lengths and amplification of the Reynolds stresses near the vicinity of the roughness element. For larger roughness elements, the downstream collapse of Reynolds stresses was sudden, whereas it extended far downstream for smaller roughness elements.

Extensive studies have looked at the effects of multiple, tandem wall-mounted structures on flow response and overall flow dynamics. Leonardi et al. (2003) [17] studied the flow response over tandem wall-mounted square bars, over a range of periodic separation ratios. It was determined that for large separation between two obstacles, the recirculation length was altered by the adverse pressure gradient imposed by the upstream wall of the subsequent object. Such flow behaviour was analogous to the study of wake characteristics behind flat plates by Hemmati et al. (2019) [18]. This was followed by the study of Goswami and Hemmati (2020) [19], who investigated the effects of turbulent flow over multiple roughness elements and different separation patters on the flow response and recovery. They considered two roughness element heights of $h/D=0.05$ and $0.1$, as well as two separation patterns: periodic and staggered. It was determined that the flow recovery was prolonged by the use of smaller roughness elements, while the separation patterns had negligible effects on the overall response and recovery as long as the elements are positioned outside the impinging zone of the upstream element. Further, Goswami and Hemmati (2021) [20] investigated the Reynolds number effects and scaling on the response and recovery of the flow over the roughness elements. They considered a square bar roughness element of $h/D=0.05$ and Reynolds numbers between $5\times {10}^3$ and $1.56\times {10}^5$. They observed an asymptotic trend in the reattachment lengths with increasing the Reynolds numbers, while the recovery trends followed a power-law behaviour of diffusion towards the centerline. Furthermore, the collapse of stresses towards the wall appeared earlier at low Reynolds numbers. While the studies on Newtonian flow past wall-mounted obstacles are abundant, viscoelastic fluid flow over such conditions has received limited attention despite its extensive industrial applications.

The numerical studies on viscoelastic flows mainly used Direct Numerical Simulations (DNS) [21,22,23,24,25] and Reynolds-Averaged-Navier–Stokes (RANS) turbulence models [26,27]. Experimentally, Particle Image Velocimetry (PIV) [28,29] and Laser-Doppler-Velocimetry (LDV) [30,31] were used to characterize the flow. Experimental studies mainly focused on investigating turbulence statistics, structure of polymeric flow, and the study of elastic instabilities over a wide range of sub-critical Reynolds numbers. Numerical studies used DNS with one or more non-linear differential models such as Oldroyd-B model [32], Giesekus model [33], Phan-Thein-Tanner (PTT) model [34], and Finitely Extensible Non-linear Elastic-Peterlin (FENE-P) model [35]. These studies were limited to lower Reynolds numbers because of the complexities in terms of modeling various rheological behaviour of the fluid and resolving various non-linear viscous and elastic effects [36]. Earlier numerical studies using DNS and several viscoelastic models were limited to smooth pipes, plane channel, and boundary layers. Azziez et al. (1996) [24] studied Giesekus, PTT and FENE-P models in simulating flow through 4:1 planar contraction. By comparing the performance of different rheological models with the experimental results, it was determined that the FENE-P model performed the best overall in terms of predicting the pressure drop, shear and normal stresses and drag reduction. This was followed by the work of Sureshkumar et al. (1997) [37], who demonstrated the first DNS channel flow using a viscoelastic fluid and FENE-P constitutive model. This study was performed at lower Reynolds numbers and a range of Weissenberg numbers ($Wi$, a dimensionless number characterizing the flow elasticity) to benchmark a set of criteria for onset of drag reduction in plane channel flows. The Finitely Extensible Non-linear Elastic (FENE) model was first proposed by Warner and Harold (1972) [38], as a solution to the gaps in the linear spring model that had no restriction on the deformation of the polymer chain. This was further extended by Bird et al. (1987) [35], who proposed the FENE-P model with Peterlin closure approximation. Here, the finite extension of the polymer chains was restricted by the parameter $L^2$, which is the dimensionless extensibility of the polymer chain. Since many such studies have complimented its accuracy and compatibility with numerical simulations, the FENE-P constitutive model is currently the most appropriate model for complex flow simulations.

Rothstein and McKinley (2001) [39] investigated the flow of polystyrene-added fluid through contraction–expansion geometry of varying expansion ratios using LDA. They studied the effects of contraction-expansion geometries on the pressure drop due to flow contraction and wake characteristics. They determined that while the addition of polymer-additive in the fluid does reduce the overall pressure and skin-friction drag, the intensity of such reduction vastly depends on the concentration of the added polymer. This was followed by Poole and Escudier (2003) [40], who investigated the influence of different polymer concentrations on the turbulent flow through sudden planar expansion using transverse LDA with a linear expansion ratio of $1.43$. They observed that the reattachment lengths by low polymer-concentration fluids increased, while for higher concentrations, it decreased significantly. The lower reattachment length was accompanied by a large reduction in turbulent intensity occurring at the point of flow separation for higher concentration solutions. Oliveira (2003) [41] further numerically investigated the effects of axisymmetric abrupt expansion geometries using FENE rheological models for low Reynolds numbers and a range of $Wi$, polymer extensibility and concentrations. The results showed similarities with the findings of Poole and Escudier (2003) [40], combined with a reduced pressure and skin-friction distribution along the walls. Further, Poole et al. (2007) [42] performed experiments on viscoelastic flow over a backward-facing step with varying polyacrylamide-concentration solutions at Reynolds numbers of 10–100. It was determined that as the polymer concentration increases, the combined effects of shear-thinning and viscoelasticity lead to reduced reattachment lengths behind the step, while an overshoot of velocity was observed compared to Newtonian fluid.

In recent years, there have been several studies that investigate the effects of viscoelastic polymer additives on flow past bluff bodies and wall-mounted obstacles, as well as on the implications of viscoelasticity on large-scale and very-large scale motions. Xiong et al. (2013) [43] carried out numerical simulations of two-dimensional flow past a circular cylinder using multiple rheological models and varying Reynolds number, Weissenberg number and polymer concentrations. An overshoot of velocity gradients and stresses were observed in the vicinity of the cylinder at high polymer conversations, while the intensities of the overshoot were very low compared to Newtonian fluid. Further, as the Reynolds number increased, vortex shedding was observed for the Newtonian fluid, which reduced gradually with increasing polymer concentrations and $Wi$ numbers. Tsukahara et al. (2014) [44] performed DNS simulations of turbulent viscoelastic flows in a channel with wall-mounted plates and studied the influence of viscoelasticity on turbulent structures and large-scale motions behind the wall-mounted plates. They observed that for the Newtonian fluid, three pairs of large-scale vortices occur behind the plates, which was significantly weakened in the viscoelastic fluid. This suppression of vortices was influenced by the elastic and viscous forces, which resulted in a significant reduction of drag behind the plate and a decreased skin-friction on the surface of the plate. They also observed a slight increase in reattachment lengths behind the plate with lower polymer concentrations and a reduction at higher concentrations. These findings were in agreement with those of Tsukahara et al. (2013) [45] and Poole and Escudier (2003) [40], who provided evidence of the mixed effects of elastic and viscous forces.

While investigations such as those of Dubief et al. (2013) [25] and Shaban et al. (2018) [29] extensively study the mechanisms of turbulent viscoelastic plane channel flows and smooth pipes, limited studies exist that investigate these effects with abrupt surface variations at high Reynolds numbers. Furthermore, viscoelastic flows past roughness elements and small-step ($h/{\delta }_o\ll 1$) conditions have not gained sufficient attention. Here, we examine the effects of viscoelastic fluid on flow over square bar roughness elements using a DNS and FENE-P rheological constitutive model at Reynolds number of $5\times {10}^3$. Particularly, we investigate the influence of higher-polymer concentration viscoelastic fluid on the flow response and recovery over perturbed flow conditions. For clarity, the results are compared with Newtonian flow over roughness element at the same Reynolds number. This paper is prepared in such a way that the problem description in Section 2 is followed by a discussion of results in Section 3 and conclusions in Section 4.

2. Problem Description

We examined the response and recovery of Newtonian and non-Newtonian turbulent fluid flow over square bar roughness elements of two different heights using DNS at Reynolds number of $5\times {10}^3$. First, the simulations were carried out using water, a Newtonian fluid, as a benchmark simulation. This was followed by examining the effects of viscoelastic polymer additives using a FENE-P rheological model incorporated into DNS. Finally, the results of viscoelastic simulations were compared with those of Newtonian DNS simulations to study the implications of viscoelastic fluid on the overall flow dynamics, response and recovery.

The governing equations for incompressible fluids are the continuity equation and the momentum equation, given by

$$\begin{array}{c}\hfill \nabla ·u=0,\end{array}$$

(1)

$$\begin{array}{c}\hfill \frac{\partial \rho u}{\partial t}+\nabla ·\left(\rho uu\right)=-\nabla ·p+\nabla ·\mathbf{\tau },\end{array}$$

(2)

where u is the velocity vector, $\rho$ is the fluid density, p is the isotropic pressure and $\mathbf{\tau }$ is the stress tensor. For Newtonian DNS simulations, the continuity (Equation (1)) and momentum (Equation (2)) equations were solved directly. For viscoelastic simulations, however, the stress tensor was split into a solvent or Newtonian stress component (${\mathbf{\tau }}_s$) and a non-Newtonian or polymeric stress component (${\mathbf{\tau }}_p$), such that

$$\mathbf{\tau }={\mathbf{\tau }}_s+{\mathbf{\tau }}_p.$$

The Newtonian stress component is defined as ${\mathbf{\tau }}_s=2{\mu }_s\mathit{D}$, where ${\mu }_s$ is the solvent viscosity and $\mathit{D}$ is the rate-of-strain tensor given by

$$\mathit{D}=\frac{1}{2}\left(\nabla ·u+\mid \nabla ·u{\mid }^T\right).$$

Substituting this into Equation (2) gives a modified momentum equation:

$$\frac{\partial \rho u}{\partial t}+\nabla ·\left(\rho uu\right)-{\mu }_s\nabla ·(\nabla u)=-\nabla ·p+\nabla ·{\mathbf{\tau }}_{\mathbf{p}},$$

(3)

where ${\mathbf{\tau }}_p$ is defined in the multimode from (${\sum }_{k=1}^n{\mathbf{\tau }}_{pk}=0$) as a symmetric tensor (${\mathbf{\tau }}_{pk}$) obtained by the sum of the contributions of individual relaxation modes. The majority of viscoelastic materials are composed of molecular structures of different sizes, and therefore they have different relaxation times. The multimode formulation enables the inclusion of a vast spectrum of relaxation times to obtain more realistic results. The expression of ${\mathbf{\tau }}_{pk}$ depends on the viscoelastic constitutive equations employed. For viscoelastic simulations here, the polymeric stress tensor was computed using the FENE-P model constitutive equation [35] as follows,

$$\left(1+\frac{\frac{3}{1-3/{L_k}^2}+\frac{{\lambda }_k}{{\mu }_{pk}}tr\left({\mathbf{\tau }}_{\mathit{pk}}\right)}{{L_k}^2}\right){\mathbf{\tau }}_{\mathit{pk}}+{\lambda }_k\stackrel{\nabla }{{\mathbf{\tau }}_{\mathit{pk}}}=2\left(\frac{1}{1-3/{L_k}^2}\right){\mu }_{pk}D,$$

(4)

where $\lambda$ is the polymer relaxation time, $L^2$ is the extensibility of the polymer molecule, tr(${\mathbf{\tau }}_p$) is the trace of the polymer stress tensor, and $\stackrel{\nabla }{{\mathbf{\tau }}_{\mathit{pk}}}$ is the upper convected derivative given as

$$\stackrel{\nabla }{{\mathbf{\tau }}_{\mathit{pk}}}=\frac{D{\mathbf{\tau }}_{\mathit{pk}}}{Dt}-[\nabla u·{\mathbf{\tau }}_{\mathit{pk}}]-{[\nabla u·{\mathbf{\tau }}_{\mathit{pk}}]}^T,$$

(5)

and $\frac{D{\mathbf{\tau }}_{\mathit{pk}}}{Dt}$ is the total derivative of the polymer-stress tensor, defined as

$$\frac{D{\mathbf{\tau }}_{\mathit{pk}}}{Dt}=\frac{\partial {\mathbf{\tau }}_{\mathit{pk}}}{\partial t}+u·\nabla {\mathbf{\tau }}_{\mathit{pk}}$$

(6)

The simulation setup in the current study mimicked those of Favero et al. (2010) [21] and Holmes et al. (2012) [23], both of which incorporated the FENE-P rheological model in DNS using OpenFOAM. They then proceeded with verifying and validating the accuracy of this model. The FENE-P model was further validated by Dubief et al. (2013) [25] in studying elasto-inertial turbulence. To this effect, we simulated the experimental channel flow of Shaban et al. (2018) [29] by incorporating the FENE-P model into DNS in OpenFOAM, which showed a good agreement in Figure 1. Hereon forward, the overbar symbol ($\overline{\varphi }$, where $\varphi$ is any flow parameter) identifies time-averaged results. The grid spacing for the simulation was $\Delta x^{+}=1.85$, the Weissenberg number ($W{i}_H$) was $1.82$, the polymer extensibility ($L^2$) was 200 and the viscosity ratio, $\beta$, was $0.22$, consistent with the experiments.

The schematics of the axisymmetric computational domain, as shown in Figure 2, was designed following the setup of Smits et al. (2019) [16] and Yamagata et al. (2014) [46], which was verified by Goswami and Hemmati (2020, 2021) [19,20] to properly capture the main features of the geometrically-symmetric pipe flow. The axisymmetric perturbation was introduced by placing a square cross section ring with heights of $h/D=0.05$ and $0.1$ in the pipe. Here, the square bar roughness element has a height of h, and the inner diameter of the pipe is D. All simulations were performed based on an axisymmetric assumption. The domain extended from $-30D$ to $+40D$ in the $x-$direction, and $R=D/2$ in the $y-$direction. A non-homogeneous spatial grid with a total of $1.2\times {10}^6$ hexahedral elements was used to discretize the domain. The grid distribution was such that the fine mesh was placed close to the wall boundaries, as shown in Figure 3. The minimum spatial grid requirement for a DNS simulation was satisfied by coinciding with the order of magnitude of element size, $O[\Delta x]$, with the order of magnitude of the Kolmogorov length scales, $O\left[\eta \right]$, at the pipe walls and in the vicinity of the roughness element. The Kolmogorov length scales are calculated using $\eta ={({{\nu }_s}^3/\epsilon )}^{1/4}$, based on Newtonian fluid viscosity, where $\epsilon$ is the mean dissipation rate calculated using the approximation of $\epsilon \approx {\nu }_s\left(\overline{\frac{\partial u_i}{\partial x_j}\frac{\partial u_j}{\partial x_i}}\right)$.

Table 1. Grid resolution, $\Delta x/\eta$, for the current DNS study.

Case	$\mathit{x}/\mathit{D}\mathbf{=}$	1	2	5	8
Current Case		$2.06$	$2.27$	$2.90$	$3.54$
Hemmati et al. (2018) [47]		$1.85$	$2.59$	$3.24$	$3.97$
Ashrafian et al. (2004) [48]		≈3.2	−	−	−
Dubief et al. (2013) [25]		≈2.3	−	−	−

shows the grid resolution of the current grid as a variation of $\Delta x/\eta$ in the axial direction up to $x/D=8$. The $\Delta x/\eta$ was less than 4 at least up to $x/D=8$, while near the vicinity of the roughness element, they were below ∼2. These ratios were comparable to other DNS studies on wakes [47] and viscoelastic pipeflows [25,48].

The parameters used in the simulations are given in

Table 2. Parameter space for the current study.

Parameter	Value
Extensibility of the polymer, $L^2$	200
Ratio of solvent to zero-shear viscosity, $\beta$	$0.22$
Polymer relations time, $\lambda$ (s)	$0.023$
Reynolds number, $R{e}_D$	5000
Frictional Reynolds number, $R{e}_{\tau }$	384

. For all viscoelastic simulations, the non-dimensional extensibility ($L^2$) and parameter $\beta$ were fixed at 200 and 0.22, where $\beta$ is defined as the ratio of solvent viscosity to the zero-shear viscosity of the polymer solution, given by $\beta ={\mu }_s/({\mu }_s+{\mu }_p)$. The parameter $\beta$ was selected such that we obtain a shear-thinning solution behaviour. Further, the corresponding polymer concentration of 160 ppm availed lowering of the solution viscosity, which enables the shear-thinning fluid behaviour analyzed in this study, as well as allowing for a higher drag reduction [29]. $\lambda$ is the polymer relaxation time, which was obtained experimentally by Shaban et al. (2018) [29] for a solution of polyacrylamide polymer at 160 ppm. The simulations were performed at $R{e}_D=5\times {10}^3$, defined by pipe diameter (D) and bulk velocity ($U_b$) across the pipe. Using the friction factor correlation given by Mckeon et al. (2004) [49], the upstream friction velocity corresponded to the friction Reynolds number of $R{e}_{\tau }=384$. The Reynolds number based on the square bar height and solvent viscosity was 250 and 500, and the linear expansion ratios were $ER=D/(D-2h)=1.11$ for $h/D=0.05$ and $ER=1.25$ for $h/D=0.1$.

Identical boundary conditions were used for all the simulations. The inlet boundary condition was a constant uniform velocity, $u=U_b$. The separation distance between the inlet and the roughness element was kept sufficiently long to satisfy the hydrodynamic entrance length condition at this Reynolds number. A Neumann-type outflow condition based on $\partial \phi /\partial n=0$, where $\phi$ is any flow variable, was imposed at the outlet, and a no-slip type boundary condition was imposed at the walls and the roughness element. The timestep, $\Delta t$, for all simulations was set so that the maximum Courant–Friedrichs–Lewy (CFL) number remained below 0.8 for Newtonian simulations and 0.2 for viscoelastic simulations. The smaller timestep requirement for viscoelastic simulations guarantees the boundedness of ${\mathbf{\tau }}_{\mathit{pk}}$ [25]. All simulations were performed using OpenFOAM, an open-source finite-volume-method platform. All discretized equations were solved using PimpleFoam, a transient solver for incompressible, turbulent flow. PimpleFoam solver incorporates the PIMPLE algorithm [50], which is the combination of PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms. The second-order accurate and bounded numerical discretization schemes for spatial discretization and the backward Euler scheme for temporal discretization were employed for the simulations. The convergence criterion of ${10}^{-6}$ was set for the root-mean-square of momentum residuals, for each timestep. The convergence was guaranteed by running all simulations for ten through-times, where a through-time is defined as the time fluid takes to travel from inlet to outlet without any disturbance. All simulations were completed using 48 Intel Platinum 8160 F Skylake 2.1 GHz cores and 192 GB of memory, using $3.2\times {10}^4$ core hours on Compute Canada clusters.

3. Results and Discussion

We begin by looking at the viscoelastic fluid flow response over the square bar roughness element. The Newtonian fluid flow simulation is our base flow to establish the differences due to non-Newtonian fluid effects. This is followed by the investigation of flow recovery and the distributions of pressure and wall shear stresses by the viscoelastic fluid.

3.1. Flow Response

The flow response was evaluated by analyzing the nature of mean wake features behind the roughness element. Figure 4 and Figure 5 show the mean streamline plots for Newtonian and viscoelastic fluid for two roughness heights. For all cases, the flow separation occurred at the upstream edge of the roughness element, followed by the formation of a recirculation bubble before the flow reattaching to the bulk flow in the downstream region of the roughness element. From the first glance, it was evident that the longitudinal width (length) of the separation bubble for viscoelastic fluid in Figure 4b and Figure 5b was significantly smaller than that of the Newtonian fluid in Figure 4a and Figure 5a, respectively. Furthermore, the lateral width (height) of the recirculation bubble exceeded that of the roughness element by ∼5% for $h/D=0.05$ and ∼15% for $h/D=0.1$ in the case of the Newtonian fluid. For viscoelastic fluid, however, the lateral width of the recirculation bubble did not exceed the height of the roughness element for both bar heights.

A quantitative comparison of the recirculation lengths obtained from Newtonian and viscoelastic flow past a roughness element is shown in

Table 3. The mean recirculation length ($\overline{L_r}/h$) for Newtonian and non-Newtonian flows over roughness elements.

Study	$\mathit{h}/\mathit{D}$	$\overline{{\mathit{L}}_{\mathit{r}}}/\mathit{h}$	$\mathbf{\Delta }\overline{{\mathit{L}}_{\mathit{r}}}(\%)$
Newtonian	$0.05$	$23.28$	−
Viscoelastic	$0.05$	$6.28$	73
Newtonian	$0.1$	$34.18$	−
Viscoelastic	$0.1$	$17.18$	$49.74$

. For the smaller roughness height, the mean recirculation length formed by the viscoelastic fluid exhibited a difference of ∼73% from that of the Newtonian fluid. For the larger roughness height, the observed difference was ∼50%. For Newtonian fluid, the reattachment length differed between the two bar heights by ∼32%. Since the larger bar element ($h/D=0.1$) created a larger contraction in the flow, it led to a larger pressure gradient compared to that of the smaller square element. Thus, both the recirculation bubble and the reattachment length were larger for Newtonian flow over roughness element of height $h/D=0.1$ compared to the smaller element. Viscoelastic fluid, on the contrary, suppressed the velocity gradients in the vicinity of the roughness element, due to which the effects of large pressure gradients were diminished. This led to a shorter recirculation length compared to that of the Newtonian fluid.

The profiles of mean axial velocity obtained from Newtonian and viscoelastic fluids are shown in Figure 6 and Figure 7 for $h/D=0.05$ and $0.1$, respectively. Profiles are plotted at three different axial locations downstream of the roughness elements: $x/h=60,80$ and 100. For Newtonian fluid (Figure 6a and Figure 7a), the profiles for different bar height cases were quite distinct. In the case of $h/D=0.05$, the velocity near the centre at $x/h=60$ was ∼30% below the fully developed profile. As the flow progressed, the variations at $x/h=100$ became negligible near the centre, while slight variations (∼0.2%) were observed closer to the wall. In contrast, velocity profiles for $h/D=0.1$ showed faster progression towards the fully developed profiles with ∼5% difference at the centre between $x/h=60$ and 100. Smits et al. (2019) [16] stated that the velocity for a smaller roughness element exceeded the fully developed profile, while that of a larger roughness element was over-flattened by Reynolds shear stresses. The trends of velocity profiles observed here contrasted the trends observed by Smits et al. (2019) [16] at $Re=1.56\times {10}^5$, which hinted at the effects of the Reynolds number on the flow dynamics and scaling. This finding agreed with experiments of Hultmark et al. (2012) [51], who recorded similar effects in pipeflow with velocity in the centre of the pipe scaling with increasing Reynolds number. The study of Goswami and Hemmati (2021) [20] showed that the flow response and recovery scaled with changing Reynolds number and observed a reduced mean centerline velocity with increasing the Reynolds number, which was consistent with the observation of the current study.

In the case of viscoelastic fluids (Figure 6b and Figure 7b), velocity profiles vastly differed from the trends observed for Newtonian cases. For both bar heights, velocities appeared to have recovered close to fully developed profiles as early as $x/h=60$. For $h/D=0.05$, the velocity at $y=R$ was slightly (∼2.5%) below that of the fully developed profile. As the flow progressed between $x/h=60$ and 100, the variation in velocity was negligible, showing a trend similar to that of Newtonian flow over the smaller roughness element ($h/D=0.05$). In the case of a larger square bar, velocity close to the center of the pipe ($y=R$) appeared to recover to a fully developed profile, while slight variations were observed in $0.2\ll y/R\ll 0.8$. The lower free-shear Reynolds number for the viscoelastic flow hinted at laminarization effects as a result of significantly reduced solution viscosity to favour shear-thinning effects. Thus, the profiles of velocity appeared to be laminar, while the peak velocities remained high.

3.2. Flow Recovery

The mean flow recovery was compared between the Newtonian and viscoelastic flows by locating the axial position where the profiles appear to approach an equilibrium state or the fully-developed state. The recovery location, $X_r/h$, was assumed if the consecutive variations in the flow gradients in the axial (streamwise) direction were less than $5\%$. The results in

Table 4. The recovery locations compared between Newtonian and viscoelastic flows, for both roughness heights.

Study	$\mathit{h}/\mathit{D}$	${\mathit{X}}_{\mathit{r}}/\mathit{h}$	$\mathbf{\Delta }{\mathit{X}}_{\mathit{r}}(\%)$
Newtonian	$0.05$	350	−
Viscoelastic	$0.05$	100	71
Newtonian	$0.1$	200	−
Viscoelastic	$0.1$	60	70

show the recovery locations, $X_r/h$, for both Newtonian and viscoelastic flow and both roughness heights. The mean Newtonian flow recovered within $X_r=350h$ for $h/D=0.05$ and $X_r=200h$ for $h/D=0.1$. The mean flow recovery of the smaller roughness element was prolonged, while that of the larger element was accelerated. These observations are consistent with the study of Goswami and Hemmati (2020) [19]. Similarly, mean viscoelastic flow recovered within $X_r=100h$ for $h/D=0.05$ and $X_r=60h$ for $h/D=0.1$. Comparable to the Newtonian flow counterpart, the mean viscoelastic flow showed faster recovery with $h/D=0.1$, while a prolonged recovery by $h/D=0.05$.

For $h/D=0.05$, the mean axial flow showed an accelerated recovery, ∼71% faster for the viscoelastic flow compared to the Newtonian flow. Similarly to $h/D=0.1$, the viscoelastic flow recovered ∼70% faster compared to the Newtonian flow. This hinted at an accelerated recovery for viscoelastic flows compared to Newtonian flow, which may be attributed to the effect of reduced fluid viscosity due to the addition of the polymeric additive, resulting in the shear-thinning effects and faster mean flow recovery.

3.3. Distribution of Pressure and Wall Shear Stresses

Effects of viscoelasticity on the flow characteristics and the roughness element were more clearly illustrated in the distributions of pressure and skin friction. First, the contribution of skin friction and pressure drag were examined on the overall local drag acting on the perturbation surfaces. The distribution of pressure and wall shear stresses were defined by pressure coefficient ($C_p$) and skin friction coefficients ($C_f$), respectively. The coefficient of pressure is given as

$$C_p=(p-p_{\infty })/\frac{1}{2}\rho {U_b}^2,$$

where $p_{\infty }$ is the static pressure of the incoming flow. Similarly, the skin friction coefficient is given as

$$C_f={\tau }_w/\frac{1}{2}\rho {U_b}^2,$$

where ${\tau }_w$ is the wall shear stress.

The distribution of mean pressure coefficient ($\overline{C_p}$) was examined in Figure 8 for both roughness heights. The profiles were drawn along the wall for $x/h=0-50$ downstream of the roughness elements. For the Newtonian fluid, the pressure rose abruptly after $x/h\sim 5$ in the case of $h/D=0.05$, while the sudden increase in pressure for $h/D=0.1$ occurred at $x/h\sim 20$. These locations roughly coincided with the centre of the recirculation bubble behind the roughness element. This was in good agreement with Leonardi et al. (2003) [17] and Hemmati et al. (2019) [18], who stated that the flow starts experiencing adverse pressure gradients downstream of the centre of the recirculation bubble. In the case of viscoelastic fluid, however, the initial increase in the pressure was noted near the stagnation point. Farther downstream, pressure reduced for the smaller bar height ($h/D=0.05$), while it stabilized for $h/D=0.1$. For the latter, pressure dropped as the flow progressed downstream. This suggests that high pressure gradients formed due to the increased near-wall flow interactions for Newtonian flow diffused quickly in the case of viscoelastic flow. Thus, it facilitated the flow recovery towards an equilibrium state faster.

Profiles of mean pressure coefficient on the upstream and downstream faces of the roughness element are shown in Figure 9. Here, it was apparent that larger roughness elements produced greater pressure loss with Newtonian fluid. This observation was consistent with the findings of Liu et al. (2019) [52] on Newtonian flow over an array of periodic roughness elements. Similarities were observed with viscoelastic flow, where pressure losses were higher for a larger element compared to a smaller element. Furthermore, pressure on surfaces of the roughness element was significantly lower for the Viscoelastic fluid compared to the Newtonian fluid for both bar heights. This hinted at a reduction in generated drag on the perturbation surfaces due to Viscoelastic effects.

The distribution of pressure and skin friction on the top surface of the perturbation is shown in Figure 10 and Figure 11. The trend for the pressure profiles followed an initial drop for all cases, which was followed by a gradual increase. These were associated with the flow separation on the top surface of the roughness element. Location of the initial collapse in pressure suggested a strong average pressure gradient forming by flow separation on the leading edge of the roughness element. This suggested a local flow acceleration, which was evident by the sharp initial rise in skin friction from Figure 11. As the flow progressed downstream, there was a sudden increase in pressure near the wall (see Figure 8) along with a drop in skin friction. Thus, the effect of viscosity became dominant in the response and recovery on mean flow. Furthermore, pressure losses were higher for the larger roughness element in both fluids. For the viscoelastic fluid, these were significantly lower than those observed for the Newtonian fluid. In the case of the latter, the pressure coefficient increased again up to $x/h=1$ for both element heights. The skin friction (Figure 11) was initially large, but then it experienced a rapid reduction to an equilibrium constant value for both bar heights, reflecting the acceleration of the flow near the roughness leading edge. Negative skin friction occurred for $h/D=0.1$, indicating the existence of a recirculating region on top of the roughness element.

For viscoelastic fluid, the initial drop in pressure was followed by a gradual increase towards relatively constant values near the downstream edge of the roughness element. The pressure losses for the larger bar were higher compared to the smaller bar, following the trend of the Newtonian fluid. The skin friction for $h/D=0.05$ dropped gradually up to $x/h=0.9$ before an abrupt drop was observed near the downstream face of the roughness. For the larger bar element ($h/D=0.1$), it increased to achieve a constant value suggesting that the flow was accelerating near the downstream edge. This may be due to the fact that, as mentioned before, the larger bar element extended outside the boundary layer, whereas the smaller bar height was submerged in the boundary layer. The magnitude of the skin friction was higher in the Viscoelastic flow, which hinted at a higher contribution of viscous stresses to the local drag induced on the roughness. This was while the contribution of pressure drag due to the larger pressure losses on the front and back surfaces of the roughness element was lowered. Thus, the overall induced drag on the roughness element was reduced with viscoelastic fluid, leading to a smaller recirculation region and a faster recovery.

4. Conclusions

The implications of viscoelastic fluid on the response and recovery of mean flow over two roughness element heights ($h/D=0.05$ and $0.1$) were examined using Direct Numerical Simulations at $Re=5\times {10}^3$. For viscoelastic simulations, a Finitely Extensible Non-linear Elastic-Peterlin (FENE-P) rheological model was considered with the rheological parameters corresponding to a high-molecular weight polymer additive of 160 ppm concentration. The response behaviour was evaluated by studying the mean wake features behind the roughness element. The streamlines for viscoelastic fluid showed a shorter recirculating region for both bar heights compared to that of the Newtonian fluid. The lateral width of the recirculating bubble did not extend above the height of the roughness element contrary to the Newtonian flow. The reattachment length was significantly reduced, e.g., ∼73% for $h/D=0.05$ and ∼50% for $h/D=0.1$. The larger roughness element created a larger contraction in the flow, which was noticeable from the Newtonian fluid. Due to the high suppression of velocity gradients in the vicinity of the roughness element, compared to Newtonian fluid, the effects of pressure gradients were not significant in viscoelastic fluid, even for the larger roughness element. Thus, it led to a shorter reattachment length. Looking at the velocity field, the progression of velocity was slower for the smaller bar height compared to the larger element in Newtonian flow. For viscoelastic fluid, the velocity near the centre of the pipe appeared to have recovered close to the fully developed profile for the smaller roughness element, while a full recovery was observed for the larger element with only slight variations. The mean flow recovery for the Newtonian fluid was slower than that of viscoelastic flow. This was consistent for both bar heights with ∼71% faster recovery of viscoelastic flow for $x/h=0.05$ compared to ∼70% for $x/D=0.1$. This accelerated recovery of mean flow was attributed to the shear-thinning effects due to the reduced solution viscosity by the addition of polymeric additives.

The distribution of pressure on the pipe wall downstream of the roughness element coincided with the centre of the separation region, after which the flow started experiencing adverse pressure gradients. The pressure for viscoelastic fluid appeared to subside to an equilibrium state faster than the Newtonian flow. The distribution of pressure on the front and back surfaces of the roughness element indicated a larger pressure loss by the viscoelastic fluid compared to the Newtonian fluid. Furthermore, pressure losses by a larger element were greater than that of the smaller element with viscoelastic flow, an observation similar to that of Newtonian flows. This suggested that the roughness height is a prominent parameter in perturbed flows independent of Newtonian or viscoelastic characteristics of the flow. The contribution of viscous stresses to the overall local drag in viscoelastic flow was higher on the top surface of the roughness element, while that of pressure drag was lowered due to large pressure losses. Thus, the overall localized drag on the roughness element was greatly reduced for viscoelastic fluid, which led to a faster recovery.