Smoothed Particle Hydrodynamics Simulations of Water Flow in a 90° Pipe Bend

Abstract

The flow through pipe bends and elbows occurs in a wide range of applications. While many experimental data are available for such flows in the literature, their numerical simulation is less abundant. Here, we present highly-resolved simulations of laminar and turbulent water flow in a ${90}^{°}$ pipe bend using Smoothed Particle Hydrodynamics (SPH) methods coupled to a Large-Eddy Simulation (LES) model for turbulence. Direct comparison with available experimental data is provided in terms of streamwise velocity profiles, turbulence intensity profiles and cross-sectional velocity maps at different stations upstream, inside and downstream of the pipe bend. The numerical results are in good agreement with the experimental data. In particular, maximum root-mean-square deviations from the experimental velocity profiles are always less than ∼1.4%. Convergence to the experimental measurements of the turbulent fluctuations is achieved by quadrupling the resolution necessary to guarantee convergence of the velocity profiles. At such resolution, the deviations from the experimental data are ∼0.8%. In addition, the cross-sectional velocity maps inside and downstream of the bend shows that the experimentally observed details of the secondary flow are also very well predicted by the numerical simulations.

1. Introduction

Pipe bends and elbows are frequently used to fit pipeline systems. They are important in many engineering applications such as the transportation of oil and gas in pipelines as well as other fluids, including hot water and steam for shorter distances, water for drinking or irrigation over long distances and hydrogen from the point of delivery to the point of demand, just to mention a few applications. The flow in bends is complex and inherently unsteady due to the generation of secondary flows in the form of vortices as a result of centrifugal forces produced by the bend curvature [1,2,3]. This secondary flow can produce strong vibrations and mechanical fatigue of the pipe system. The strength of the secondary flow will depend on the radius of curvature of the pipe bend, $R_{\mathrm{c}}$, and the Reynolds number, Re $=v_{\mathrm{B}}D/\nu$, where $v_{\mathrm{B}}$ is the bulk velocity, D is the pipe diameter and $\nu$ is the kinematic viscosity. Understanding its characteristics is important to determine the pressure losses and the heat/mass transfer in such bends, which are of considerable engineering importance. For example, visualization of turbulent motion in pipe flow has determined that flow separation occurs behind the bend inlet when $R_{\mathrm{c}}/D<1.5$, giving rise to increased pressure losses, while for $R_{\mathrm{c}}/D>1.5$ a pair of counter-rotating vortices forms within the bend that force the cross-sectional contours of the streamwise flow velocity to become C-shaped and its profiles to substantially distort from parabolic flow as they are shifted away from the center of curvature of the bend [4].

Experimental investigations of laminar and turbulent flow in ${90}^{°}$ pipe bends of circular and square section have been reported by a number of authors in the last four decades [1,2,3,5,6,7,8,9,10]. An extensive review on theoretical, experimental and numerical work on flow in curved pipes spanning from 1839 to 2016 was reported by Kalpakli Vester et al. [11]. Some of these studies include numerical simulations of the experiments using commercial fluid flow softwares such as ANSYS Fluent and OpenFOAM. A numerical assessment of the characteristics of laminar flow in a ${90}^{°}$ elbow was undertaken by van de Vosse et al. [12] by adopting finite element methods and, more recently, by Pantokratoras [13] using ANSYS Fluent, who found that the velocity profiles at the bend inlet are shifted towards the inner pipe wall when the bend curvature is high, while they remain almost symmetric at low curvature. In addition, at high curvature and Reynolds numbers, the velocity profiles are shifted towards the outer pipe wall at the bend outlet and towards the inner wall at low Reynolds numbers, while when the curvature is low the velocity profiles are shifted towards the outer wall at the bend exit independently of the Reynolds number. Numerical simulations of turbulent flow in ${90}^{°}$ bends aimed at comparing the Large-Eddy Simulation (LES) approach with various Reynolds-averaged Navier–Stokes (RANS) turbulence models were reported by Röhrig et al. [14], while flow separation in such bends was numerically studied by Dutta et al. [15] using the RANS approach coupled to a $k-ϵ$ model. In terms of the mean velocity profiles and turbulence intensities downstream of the bend outlet, the former authors found that LES provided much better agreement with the experimental data than the RANS approaches. Direct numerical simulations of turbulent flow in a ${90}^{°}$ pipe bend at Re $=5300$ were performed by Wang et al. [16] to investigate flow oscillation frequencies downstream of the bend and by Hufganel et al. [17] for flow at Re = 11,700 to study the nature of the three-dimensional wave-like structure that is responsible for the low-frequency switching affecting the Dean vortices downstream of the pipe bend and causing the fatigue of the pipe system. More recently, Du et al. [18] employed a LES turbulence model together with the proper orthogonal decomposition method to investigate the effects of varying the bend angle on pressure drop and mean flow characteristics in a small-diameter corrugated pipe. They found that increasing the bend angle from $0^{°}$ to ${90}^{°}$ produces stronger Dean vortices and increased pressure drops downstream of the pipe bend.

Flow simulations in curved pipes using Smoothed Particle Hydrodynamics (SPH) methods are much less abundant. The few examples found in the literature include the two-dimensional simulations by Hou et al. [19], who studied flow separation in right-angled bends for different turning angles; the three-dimensional SPH simulations by Alvarado-Rodríguez et al. [20], where the numerical results were compared with Sudo et al.’s [2] experimental measurements for turbulent flow in a square-sectioned ${90}^{°}$ pipe bend and Rup et al.’s [6] Fluent-based simulations of Sudo et al.’s experiment; and the calculations by Rosić et al. [21] aimed at testing the SPH capability to model supercritical flow across pipe bends. SPH is a fully Lagrangian, particle-based method that has become very popular in the last thirty years because of its many applications and ease of implementation [22,23,24]. In SPH, the state of the fluid is represented by a finite set of particles, where the physical properties carried by each particle are obtained from its close neighbors by means of an interpolation function, called the smoothing kernel. This way particles are allowed to interact with each other within the compact spherical support of the smoothing kernel. The resulting collective motion of the SPH particles mimics the flow of a fluid and so it can be modeled using the classical equations of hydrodynamics.

In this paper, we present the results of highly resolved SPH simulations of laser-Doppler velocimetry measurements of laminar and developing turbulent water flow in a ${90}^{°}$ pipe elbow reported by Enayet et al. [1]. These experiments were carried out for laminar flow at Re $=500$ and 1093 and for turbulent flow at a maximum obtainable Re of 43,000. The SPH calculations were performed using a weakly compressible (WCSPH) scheme [25] coupled to a Large-Eddy Simulation (LES) model for turbulence and non-reflecting outflow boundary conditions at the pipe outlet [20]. The numerical results are compared with the experimental data for both types of flow in terms of mean velocity profiles upstream of the bend inlet, cross-sectional mean velocity maps at different angles within the bend and downstream of the bend exit and turbulence intensity profiles upstream of the bend for Re = 43,000. At spatial resolutions higher than that for which convergence to the experimental velocity profiles is guaranteed, the results reproduce Enayet et al.’s [1] measurements with root-mean-square deviations $<1$%. In particular, the wall static pressure variation within the bend and the turbulence intensity profiles upstream of the bend inlet for Re = 43,000 are remarkably well reproduced by the SPH simulations when the resolution needed to guarantee convergence of the velocity profiles to the experimental data is quadrupled.

In what follows, the governing equations and a brief outline of the numerical framework are given in Section 2, while Section 3 describes the pipe flow model. Section 4 presents the numerical results by providing direct comparison with Enayet et al.’s [1] experiments and Section 5 summarizes the relevant conclusions.

2. Governing Equations and Numerical Framework

2.1. Exact Equations

Under prescribed inlet, outlet and solid wall boundary conditions, the differential equations describing the flow in a bent pipe are given by the Navier–Stokes equations

$$\begin{array}{ccc}\hfill \frac{d\rho }{dt}& =& -\rho \nabla ·\mathbf{v},\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \frac{d\mathbf{v}}{dt}& =& -\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{v}+\mathbf{g},\hfill \end{array}$$

(2)

where $\rho$ is the fluid density, $\mathbf{v}$ is the flow velocity vector, p is the pressure, $\nu =1.005\times {10}^{-6}$ m${}^2$ s${}^{-1}$ is the kinematic viscosity of water at 20 ${}^{°}$C and $\mathbf{g}$ is the gravitational acceleration. To complete the governing equations, the Murnaghan–Tait equation of state is employed where the pressure and density are related according to

$$p=p_0\left[{\left(\frac{\rho }{{\rho }_0}\right)}^{\gamma }-1\right],$$

(3)

where $\gamma =7$ for water, $p_0=c_0^2{\rho }_0/\gamma$, ${\rho }_0=998.29$ kg m${}^{-3}$ is the water density at 20 ${}^{°}$C and $c_0$ is a numerical sound speed taken to be at least 10 times higher than the maximum flow velocity in order to ensure weak compressibility and fluctuations of the density field $|\rho -{\rho }_0|/{\rho }_0\lesssim 0.01$ [25]. This enforces the incompressibility condition because the effects of any compression induced by such density fluctuations will be purely acoustic and superimposed to the main flow with almost no interaction.

2.2. LES Filtering

In the LES approach, the flow velocity must be separated into its mean component (or resolved scale) velocity, namely $\tilde{\mathbf{v}}$, and its fluctuating component (or sub-particle scale) velocity, ${\mathbf{v}}^{\prime }$, such that $\mathbf{v}=\tilde{\mathbf{v}}+{\mathbf{v}}^{\prime }$. The mean velocity is obtained by means of the density-weighted Favre-filtering

$$\tilde{\mathbf{v}}=\frac{1}{\overline{\rho }}\frac{1}{T}{\int }_t^{t+T}\rho (\mathbf{x},t)\mathbf{v}(\mathbf{x},t)dt,$$

(4)

where T is a sufficiently large time interval and $\overline{\rho }$ is the Reynolds-averaged density. After application of the Favre filtering, Equations (1) and (2) become

$$\begin{array}{ccc}\hfill \frac{d\overline{\rho }}{dt}& =& -\overline{\rho }\nabla ·\tilde{\mathbf{v}},\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill \frac{d\tilde{\mathbf{v}}}{dt}& =& -\frac{1}{\overline{\rho }}\nabla \overline{p}+\frac{\nu }{\overline{\rho }}\left[\nabla ·\left(\overline{\rho }\nabla \right)\right]\tilde{\mathbf{v}}+\frac{1}{\overline{\rho }}\nabla ·\mathbf{T}+\mathbf{g}.\hfill \end{array}$$

(6)

Here, $\mathbf{T}$ is the sub-particle stress tensor defined by

$$T_{ij}=\overline{\rho }{\nu }_t\left(2{\tilde{S}}_{ij}-\frac{2}{3}{\tilde{S}}_{kk}{\delta }_{ij}\right)-\frac{2}{3}\overline{\rho }{C}_I{\Delta }^2{\delta }_{ij}{\left|\tilde{S}\right|}^2,$$

(7)

where ${\tilde{S}}_{ij}$ is the Favre-filtered strain rate tensor given by

$${\tilde{S}}_{ij}=\frac{1}{2}\left(\frac{\partial {\tilde{v}}_i}{\partial x_j}+\frac{\partial {\tilde{v}}_j}{\partial x_i}\right),$$

(8)

$C_I=0.00066$, ${\nu }_t={(C_s\Delta )}^2\left|\tilde{S}\right|$, with $C_s=0.12$, is the Smagorinsky eddy viscosity, $|\tilde{S}|={\left(2{\tilde{S}}_{ij}{\tilde{S}}_{ij}\right)}^{1/2}$ is the local strain rate, ${\delta }_{ij}$ is the Kronecker delta and $\Delta$ is a measure of the finite particle size. For practical purposes, $\Delta$ is set equal to the local particle smoothing length h. A form of ${\nu }_t$ that improves the stress–strain relationship is given by ${\nu }_t=C_d{\nabla }^2\left|\tilde{S}\right|$, where $C_d$ is dynamically determined as a function of the smallest resolved scales [26]. The choice of the LES model is based on the results of the comparative assessment performed by Röhrig et al. [14] between the LES and various RANS models for turbulent flow in a ${90}^{°}$ pipe elbow. In particular, they found that LES reproduced the experimental velocity profiles and the turbulent intensities in the streamwise direction much better than the RANS schemes. On the other hand, it is well known that, for flow problems where experimental data are not available, LES is in general employed to generate reference data for the calibration and validation of turbulence closures.

2.3. The SPH Method

The Favre-filtered Equations (5) and (6) are solved using a variant of the DualSPHysics code [27]. Equation (5) is solved for the density fluctuations using the SPH representation

$$\frac{d{\rho }_a}{dt}=\sum _{b=1}^n{m}_b\left({\mathbf{v}}_a-{\mathbf{v}}_b\right)·{\nabla }_a{W}_{ab},$$

(9)

where ${\rho }_a$ is the particle-scale density, $W_{ab}=W\left(\right|{\mathbf{x}}_a-{\mathbf{x}}_b|,h_a)$ is the kernel function, $h_a$ is the smoothing length of particle a, $m_b$ is the mass of neighboring particle b and the summation is taken over all neighbors of particle a within the compact support of the kernel, i.e., for $|{\mathbf{x}}_a-{\mathbf{x}}_b|\le h$. The subscript a in the nabla operator means that the derivatives of the kernel function must be taken with respect to the coordinates of particle a. The tilde operator over the ensemble average velocity vector and the bar operator over the mean density have been dropped for simplicity. In SPH form, the momentum Equation (6) is given by

$$\begin{array}{ccc}\hfill \frac{d{\mathbf{v}}_a}{dt}=& -& \frac{1}{{\rho }_a}\sum _{b=1}^n\frac{m_b}{{\rho }_b}\left(p_a+p_b\right){\nabla }_a{W}_{ab}+4\nu \sum _{b=1}^n{m}_b\frac{{\mathbf{v}}_a-{\mathbf{v}}_b}{{\rho }_a+{\rho }_b}\frac{{\mathbf{x}}_{ab}·{\nabla }_a{W}_{ab}}{|{\mathbf{x}}_{ab}{|}^2+{ϵ}^2}\hfill \\ & +& \sum _{b=1}^n{m}_b\left(\frac{{\mathbf{T}}_a}{{\rho }_a^2}+\frac{{\mathbf{T}}_b}{{\rho }_b^2}\right)·{\nabla }_a{W}_{ab}+\mathbf{g},\hfill \end{array}$$

(10)

where the symmetric form suggested by Colagrossi and Landrini [28] is employed for the pressure gradient term and those used by Lo and Shao [29] are adopted for the laminar viscous term and the sub-particle stresses. In particular, the form of the pressure gradient with a positive sign is chosen on the basis that it is variationally consistent with the use of Equation (9) for the density [30]. Here, ${\mathbf{x}}_{ab}={\mathbf{x}}_a-{\mathbf{x}}_b$ and ${ϵ}^2=0.01h_a^2$. To prevent the growth of numerical errors due to anisotropies in the spatial distribution of particles, they are moved along the pipe section using the equation [31]

$$\frac{d{\mathbf{x}}_a}{dt}={\mathbf{v}}_a+\frac{\beta x_0{v}_{\max }}{M}\sum _{b=1}^N{m}_b\frac{{\mathbf{x}}_{ab}}{{\left({\mathbf{x}}_{ab}·{\mathbf{x}}_{ab}\right)}^{3/2}},$$

(11)

where $\beta =0.04$, $v_{\max }$ is the maximum estimated fluid velocity, M is the total mass of fluid within the pipe and $x_0$ is given by

$$x_0=\frac{1}{N}\sum _{b=1}^N{\left({\mathbf{x}}_{ab}·{\mathbf{x}}_{ab}\right)}^{1/2},$$

(12)

where the summations in Equations (11) and (12) are taken over all fluid particles in the computational domain.

To improve the convergence properties of SPH, the interpolation kernel used in this paper is a Wendland C${}^2$ function [32,33]

$$W(q,h)=\frac{21}{2\pi h^3}{\left(1-q\right)}^4\left(1+4q\right),$$

(13)

for $q\le 1$ and zero otherwise, where $q=|\mathbf{x}-{\mathbf{x}}^{\prime }|/h$.

2.4. Boundary Conditions

The fluid flow in a truncated pipe requires the treatment of boundary conditions at the pipe inlet, the pipe outlet and the solid wall. No-slip boundary conditions ($\mathbf{v}=\mathbf{0}$) are implemented at the solid wall using the method of dynamic boundary particles proposed by Crespo et al. [34,35]. In this work, a layer of particles defines the pipe wall and two layers of ghost particles are placed outside the computational domain around the pipe wall. These outer boundary particles are treated as normal fluid particles. Although they are updated using Equation (10), they are forced to remain stationary and preserve their initial positions. The density of the wall and outer boundary particles is determined by the density of all neighboring fluid particles lying within their support domains. To avoid particle penetration across the pipe wall, the wall particles will exert artificial repulsive forces to nearby fluid particles, which are derived from the source of the momentum Equation (10).

The pipe inlet and outlet are treated as open boundary conditions. Proper inlet boundary conditions are implemented by placing an inflow zone in front of the entrance plane, as shown in Figure 1. The inflow zone is filled with five columns of regularly distributed particles and has a length equal to $5\Delta x_0$, where $\Delta x_0$ is the initial fluid particle separation in the flow direction. Inflow particles are then allowed to flow in as needed with a prescribed density and velocity. This way boundary information propagates into the flow domain since fluid particles close to the inlet plane will use information from neighboring inflow particles lying within their support domains. The outlet boundary is treated using non-reflecting outflow boundary conditions [20]. This method requires defining an outflow zone in front of the outlet plane of the same length as the inflow zone (see Figure 1). Fluid particles entering the outflow zone are moved by solving the outgoing wave equation

$$\frac{\partial \mathbf{v}}{\partial t}+v_x\frac{\partial \mathbf{v}}{\partial x}-\nu \left(\frac{{\partial }^2\mathbf{v}}{\partial y^2}+\frac{{\partial }^2\mathbf{v}}{\partial z^2}\right)=\mathbf{0},$$

(14)

where $\mathbf{v}=(v_x,v_y,v_z)$. In this equation, the assumption is made that the streamwise flow is in the x-direction. Since the inflow rate may differ from the outflow rate, the number of particles entering and leaving the fluid domain may not be the same. Therefore, particles leaving the outflow zone are stored in a reservoir buffer zone with a zero velocity. This way, when an inflow particle crosses the inlet plane and enters the fluid domain, a particle is removed from the buffer and placed into the inflow zone. At the beginning of a simulation, the reservoir buffer is not necessarily empty but must contain an arbitrary number of stationary particles as needed.

A numerically stable SPH representation of Equation (14) is given by

$$\frac{\partial {\mathbf{v}}_o}{\partial t}=-v_{x,o}\sum _{b=1}^n\frac{m_b}{{\overline{\rho }}_{ob}}\left({\mathbf{v}}_b-{\mathbf{v}}_o\right)\frac{\partial W_{ob}}{\partial x_o}+2\nu \sum _{b=1}^n\frac{m_b}{{\rho }_b}\frac{\left({\mathbf{v}}_b-{\mathbf{v}}_o\right)}{|{\mathbf{x}}_{ob}{|}^2+{ϵ}^2}\left(y_{ob}\frac{\partial W_{ob}}{\partial y_o}+z_{ob}\frac{\partial W_{ob}}{\partial z_o}\right),$$

(15)

where the subscript o is employed to denote outflow particles, ${\mathbf{x}}_{ob}={\mathbf{x}}_o-{\mathbf{x}}_b$, $y_{ob}=y_o-y_b$, $z_{ob}=z_o-z_b$ and ${\overline{\rho }}_{ob}=({\rho }_o+{\rho }_b)/2$. The numerical stability of Equation (15) can be improved if the streamwise velocity of outflowing particles, $v_{x,o}$, is smoothed according to the summation

$$v_{x,o}=\sum _{b=1}^n\frac{m_b}{{\rho }_b}{v}_{x,b}{W}_{ob},$$

(16)

which is equivalent to averaging the convective velocity close to the outlet plane. Since for outflow particles nearby the outlet plane the summations in Equations (15) and (16) will include both types of particles (i.e., fluid and outflow particles) the information carried by fluid particles crossing the outlet plane is propagated outwards so that feedback noise from the outlet boundary is greatly reduced. As a final remark, the position of outflow particles is determined by solving the equation

$$\frac{d{\mathbf{x}}_o}{dt}={\mathbf{v}}_o,$$

(17)

which must be integrated using the same loop as the inner fluid particles.

2.5. Time Marching Scheme

The temporal integration of Equations (9)–(11) is performed using a second-order Verlet scheme, which is implemented according to the following difference formulae

$$\begin{array}{ccc}\hfill {\rho }_a^{n+1}& =& {\rho }_a^{n-1}+2\Delta t{\left(\frac{d{\rho }_a}{dt}\right)}^n,\hfill \\ \hfill {\mathbf{v}}_a^{n+1}& =& {\mathbf{v}}_a^{n-1}+2\Delta t{\left(\frac{d{\mathbf{v}}_a}{dt}\right)}^n,\hfill \\ \hfill {\mathbf{x}}_a^{n+1}& =& {\mathbf{x}}_a^n+\Delta t{\mathbf{v}}_a^n+0.5\Delta t^2{\left(\frac{d{\mathbf{v}}_a}{dt}\right)}^n,\hfill \end{array}$$

(18)

where $\Delta t=t^{n+1}-t^n$. To improve the numerical coupling among Equations (9)–(11) during the evolution, the above scheme is replaced every 40 timesteps by the alternative differences

$$\begin{array}{ccc}\hfill {\rho }_a^{n+1}& =& {\rho }_a^n+\Delta t{\left(\frac{d{\rho }_a}{dt}\right)}^n,\hfill \\ \hfill {\mathbf{v}}_a^{n+1}& =& {\mathbf{v}}_a^n+\Delta t{\left(\frac{d{\mathbf{v}}_a}{dt}\right)}^n,\hfill \\ \hfill {\mathbf{x}}_a^{n+1}& =& {\mathbf{x}}_a^n+\Delta t{\mathbf{v}}_a^n+0.5\Delta t^2{\left(\frac{d{\mathbf{v}}_a}{dt}\right)}^n.\hfill \end{array}$$

(19)

To ensure numerical stability, the time step is constrained by the following criteria

$$\begin{array}{ccc}\hfill \Delta t_f& =& 0.25\underset{a}{min}{\left(\frac{h}{|d{\mathbf{v}}_a/dt|}\right)}^{1/2},\hfill \\ \hfill \Delta t_{v,a}& =& \underset{b}{max}\left|\frac{h{\mathbf{x}}_{ab}·{\mathbf{v}}_{ab}}{{\mathbf{x}}_{ab}·{\mathbf{x}}_{ab}+{ϵ}^2}\right|,\hfill \\ \hfill \Delta t_{\nu }& =& 0.40\underset{a}{min}\left(\frac{h}{c_a+\Delta t_{v,a}+2.5\nu /h}\right),\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill \Delta t& =& min\left(\Delta t_f,\Delta t_{\nu }\right),\hfill \end{array}$$

(21)

where the maximum and minima are taken over all particles in the system, $c_a$ is the sound speed for particle a given by

$$c_a^2=c_0^2{\left(\frac{{\rho }_a}{{\rho }_0}\right)}^{\gamma -1},$$

(22)

and

$$\left|\frac{d{\mathbf{v}}_a}{dt}\right|={\left[{\left(\frac{d{v}_{x,a}}{dt}\right)}^2+{\left(\frac{d{v}_{y,a}}{dt}\right)}^2+{\left(\frac{d{v}_{z,a}}{dt}\right)}^2\right]}^{1/2},$$

(23)

is the magnitude of the fluid acceleration obtained by direct evaluation of the source terms in Equation (10).

3. Pipe Bend Model Description

The flow configuration investigated here corresponds to the experimental test bend analyzed by Enayet et al. [1]. This represents a pipe of diameter $D=48$ mm, which is turned around an angle of ${90}^{°}$ with a radius of curvature $R_{\mathrm{c}}=2.8D=134.4$ mm, as depicted in Figure 2. Consistent entrance and exit flow is established by attaching the bend to linear piping. Straight pipes of length $5D=240$ mm and $10D=480$ mm are attached to the bend upstream of its inlet plane and downstream of its outlet plane, respectively. The numerical simulations were aimed at reproducing the laser-Doppler data acquired by Enayet et al. [1]. To do so, the simulations were carried out for three different sets of inlet conditions, corresponding to flat inlet streamwise velocities: (a) $v_x=10.5$ mm s${}^{-1}$ (Re $=500$); (b) $v_x=23$ mm s${}^{-1}$ (Re $=1093$); and (c) $v_x=0.92$ m s${}^{-1}$ (Re = 43,000). For these models, the mass flow rate is calculated according to

$$\dot{m}=\frac{\pi }{4}\rho \nu D\mathrm{Re},$$

(24)

and the bulk velocity is given by

$$v_{\mathrm{B}}=\frac{4\dot{m}}{\pi \rho D^2}.$$

(25)

To enforce direct comparison with the experimental data, flow measurements are made at the pipe stations $S_1=1D$ and $S_2=0.58D$ upstream of the bend inlet, corresponding to distances from the pipe inlet of 192 and 212.16 mm, respectively, at angles $\theta ={15}^{°}$, ${30}^{°}$, ${60}^{°}$ and ${75}^{°}$ inside the bend, and at stations $S_3=1D$, $S_4=2D$ and $S_5=6D$ downstream of the bend outlet, corresponding to distances of 432, 384 and 192 mm from the pipe exit, respectively.

Important quantities that must be compared with the experimental data for Re = 43,000 are the turbulence intensity and the pressure coefficient given by

$$C_{\mathrm{p}}=\frac{2(p-p_{\mathrm{ref}})}{\rho v_{\mathrm{B}}^2},$$

(26)

where $p_{\mathrm{ref}}$ is the wall static pressure measured in station $S_1$. In particular, we are interested in the turbulence intensity in the streamwise direction. In turbulent flow, the streamwise velocity can be separated into its mean value, $\tilde{v}$, and its fluctuating part, $v^{\prime }$, such that $v\left(t\right)=\tilde{v}+v^{\prime }\left(t\right)$, where $\tilde{v}$ is calculated as a mean over time, i.e.,

$$\tilde{v}=\frac{1}{T}{\int }_t^{t+T}v\left(t\right)dt.$$

(27)

Here, T is a sufficiently long time interval which can be set equal to the time required to achieve a stationary flow. From knowledge of $\tilde{v}$, the turbulent fluctuation in the streamwise direction follows as

$$v^{\prime }\left(t\right)=v\left(t\right)-\tilde{v}.$$

(28)

The turbulence intensity is defined as the root-mean-square (rms) of the square of the fluctuating velocity, namely

$$\overline{v^{\prime }{\left(t\right)}^2}=\frac{1}{T}{\int }_t^{t+T}{\left[v\left(t\right)-\tilde{v}\right]}^{2}dt,$$

(29)

so that

$$v_{\mathrm{rms}}^{\prime }=\sqrt{\overline{v^{\prime }{\left(t\right)}^2}},$$

(30)

which gives the turbulence intensity in the streamwise direction.

Validation and Particle Independence Test

For validation purposes, the same geometrical configuration of Figure 2 was adopted in preliminary simulations, which were optimized via a convergence study by increasing the number of SPH fluid particles. Figure 3 shows the results of the convergence study for varying spatial resolution from $N=$ 1,080,924 to 4,042,185 particles, where the streamwise velocity profiles in the horizontal (y) and vertical (z) planes normalized to $v_{\mathrm{B}}$ are compared with Enayet et al.’s [1] experimental measurements at pipe station $S_2$ for Re = 43,000. Asymptotic convergence to the experimental data is observed for $N=$ 4,042,185 particles. At this resolution, the root-mean-square (RMSE) deviations between the numerical results and the experimental data are ≈1.3% for the horizontal profile and ≈1.4% for the vertical profile.

4. Results

The validation test was performed to assess convergence to the experimental streamwise velocity profiles at pipe station $S_2=0.58D$ upstream of the bend for Re = 43,000. Further comparison of the numerical results with Enayet et al.’s [1] measurements at different angles within the bend and stations downstream of the bend outlet was performed by further increasing the spatial resolution to $N=$ 8,064,442 and 16,056,513 particles. In particular, Enayet et al.’s experimental measurements were considered to be accurate enough to provide a useful benchmark test for checking the numerical accuracy of any flow prediction method.

4.1. Laminar Flow

For laminar flow at Re$=500$, Figure 4 shows the horizontal (left) and vertical (right) streamwise velocity profiles as compared with the experimental measurements at $S_2=0.58D$ upstream of the bend inlet. The same plots for Re $=1093$ are displayed in Figure 5. The profiles for $N=$ 4,042,185, 8,064,442 and 16,056,513 particles closely overlap each other and match the experimental data with RMSE deviations ≲0.8%. When the Reynolds number is doubled from 500 to 1093, the thickness of the inlet boundary layer decreases from ≈0.26D to ≈0.23D, resulting in a larger central region (potential core) of uniform velocity. These features are also in very good agreement with the experimental observations.

Figure 6 shows maps of the streamwise velocity for Re $=500$ (left) and 1093 (right) in three cross-stream planes at $\theta ={30}^{°}$, ${60}^{°}$ and ${75}^{°}$ in the bend and at $S_3=1D$ downstream of the bend outlet for the $N=$ 16,056,513 run. Almost indistinguishable plots to Figure 6 were also obtained for $N=$ 8,064,442 particles. The velocity maps reproduce the shape of Enayet et al.’s [1] experimental contours at the same stations and the values of $v/v_{\mathrm{B}}$ associated with the colors in the maps fit very well their contour level numbers (see their Figure 4 and Figure 6). As a result of doubling the Reynolds number, the transverse pressure gradient within the bend increases and, as a consequence, the secondary flow becomes more intense. In both cases, the secondary flow develops gradually and displaces the core region of maximum velocity to the outer arc of the bend. At the $\theta ={30}^{°}$ plane, a thickening of the shear layer towards the inner arc of the bend is already visible. This effect is greater for Re $=1093$ than for Re $=500$. When the flow crosses the $\theta ={60}^{°}$ plane, the secondary flow is almost fully developed in both cases, leading to a C-shaped distortion of the core and a displacement of the peak velocity from the plane of symmetry. Owing to the increased centrifugal forces and transverse pressure gradients, the secondary flow is stronger for the Re $=1093$ case. At $\theta ={75}^{°}$, the flow distortion is such that the region of maximum velocity takes a horseshoe shape, which is more clearly defined for Re $=1093$ than for Re $=500$ because of the more intense swirling secondary flow in the former case. At the exit of the bend, the secondary motion persists and does not show a significant decay at station $S_3=1D$ downstream of the bend outlet, as shown in Figure 6.

4.2. Turbulent Flow

At Re = 43,000, the inlet boundary layers are turbulent. Figure 7 displays the horizontal (left) and vertical (right) streamwise velocity profiles for $N=$ 4,042,185, 8,064,442 and 16,056,513 particles as compared to the experimental measurements (empty circles) at $S_2=0.58D$ upstream of the bend inlet. As in Figure 4 and Figure 5 for the laminar flows, at these high resolutions, the profiles closely overlap each other, meaning that the simulations are well inside the convergence region. It is evident from the results in Figure 7 that, when passing from Re $=1093$ to Re = 43,000, the inlet boundary layers become much thinner. The simulations predict a thickness of ≈0.09D, which is in good agreement with the experimental observations for the same turbulent flow. As a consequence of the slimming of the inlet boundary layers, the region of maximum (uniform) velocity at the center of the pipe section becomes much larger, extending over most of the pipe diameter. The RMSE deviation of the numerical profiles from the experimental data is always less than about 0.7%. The horizontal and vertical profiles of the streamwise velocity rms fluctuations (turbulence intensity) are displayed in Figure 8. While asymptotic convergence for the streamwise velocity profiles in both the laminar and turbulent approach flows 0.58 pipe diameters upstream of the bend inlet is always achieved with $N=$ 4,042,185 particles, Figure 8 shows that the same level of convergence for the rms fluctuations needs to quadruple the number of particles. This clearly means that the turbulent fluctuations are by far more sensitive to resolution than the streamwise flow. The turbulence intensity is minimum in the potential core region and increases in the inlet boundary layers where the flow is turbulent. The numerical simulations match the experimental data for the turbulence intensity with RMSE deviations which are always ≲0.8% when $N=$ 16,056,513.

Cross-sectional maps of the streamwise velocity are shown in Figure 9 for $\theta ={30}^{°}$, ${60}^{°}$ and ${75}^{°}$ in the bend and in planes $S_3=1D$ and $S_5=6D$ downstream of the bend outlet. By the $\theta ={30}^{°}$ cross-sectional plane, the region of maximum velocity (potential core) was displaced towards the inside of the bend as a result of the adjustment of the flow to the cross-stream pressure gradient imposed by the bend curvature. This feature is in perfect agreement with the experimental contours of Enayet et al. (their Figure 9). At $\theta ={60}^{°}$ and ${75}^{°}$, the potential core looks highly distorted due to the presence of a strong swirling secondary flow, which is largely confined near the inner wall of the pipe where the velocity gradients are steeper and grows with distance through the bend. Close to the bend exit, the secondary flow extends over the entire cross-sectional area. However, downstream of the bend outlet, the transverse pressure gradients disappear, causing the velocity peak to move towards the outer wall as shown at pipe station $S_3=1D$. In contrast to laminar flow, where a C-shaped core forms within the bend, in turbulent flow, the core becomes C-shaped at the exit of the pipe bend. The secondary flow persists down the straight pipe section and at station $S_5$ the potential core embraces most of the pipe circular cross-section and deforms into a well-pronounced horseshoe shape.

The static pressure at the pipe wall can be calculated using Equation (26) as a function of position along the pipe and azimuthal angle $\varphi$ on the circular pipe wall. At Re = 43,000, the variation of the pressure coefficient is sufficiently large to allow for accurate experimental measurements of the wall static pressure. To provide a comparison with the experimentally obtained results of Enayet et al. (their Figure 7), we calculate the pressure coefficient at azimuthal angles $\varphi =0^{°}$, ${45}^{°}$,${90}^{°}$, ${135}^{°}$ and ${180}^{°}$ on the circular wall, as shown by the inset of Figure 10, at stations $S_2=0.58D$ upstream of the bend inlet, at cross-sections $\theta ={15}^{°}$, ${30}^{°}$, ${60}^{°}$ and ${75}^{°}$ in the bend and at stations $S_3=1D$ and $S_4=2D$ in the straight pipe section downstream of the bend outlet.

It is clear from the results in Figure 10 that the variation of the static wall pressure declines rapidly at about one diameter from the bend exit. The numerical results (filled circles) match fairly well the experimental data (open circles) with RMSE deviations between ∼1.8 ($\varphi =0^{°}$) and ∼0.8% ($\varphi ={180}^{°}$). An interesting experimental observation which is reproduced by the numerical simulations is that the negative pressure variation (at $0^{°}$ on the inner arc of the bend) is approximately twice the positive pressure gradient on the opposite side, i.e., at ${180}^{°}$ on the outer arc of the bend. The form of the pressure coefficient across the bend reveals the sequences of favorable to adverse pressure gradient along the inner and outer arcs of the bend. A comparison of both sequences shows that the pressure variation is stronger along the inner arc, which is a common feature for this type of flow.

5. Conclusions

The laminar and turbulent flows of water through a ${90}^{°}$ pipe bend were simulated numerically using Smoothed Particle Hydrodynamics (SPH) techniques coupled to a large-eddy simulation (LES) approach and non-reflective outflow boundary conditions at the pipe exit. The simulations were addressed to embody a comparative assessment of how well the proposed SPH method can predict the laser-Doppler measurements of water flow in a ${90}^{°}$ pipe elbow reported by Enayet et al. [1].

For laminar flow at Re $=500$ and 1093, the SPH simulations were able to reproduce the experimental streamwise velocity profiles upstream of the bend inlet with sufficiently good accuracy. For a total number of particles ≥4,042,185, the numerically obtained profiles were seen to deviate from the experimental data with root-mean-square errors (RMSE) that were always less than 1%. Similar accurate results were also obtained for the velocity profiles upstream of the bend inlet at Re = 43,000. At this Re-value, the simulations predict that the inlet boundary layer is turbulent and much thinner (≈0.09D, where D is the pipe diameter) than that for the laminar case at Re $=1093$ (≈0.23D), which in both cases matched the experimentally measured values. Although asymptotic convergence to the experimental streamwise velocity profiles was already achieved for $N=$ 4,042,185 particles for both laminar and turbulent flows, it was necessary to increase the resolution to $N=$ 16,056,513 particles to observe the same level of convergence for the turbulence intensity profiles in the streamwise direction close to the bend inlet. At such high resolution, the numerically obtained profiles were seen to deviate from the experimental data with RMSE errors of ∼0.8%. Moreover, the details of the flow structure through the bend, concerning the strength of the swirling secondary flow and the deformation of the potential core, were also very well predicted by the present SPH simulations. Intensity velocity maps at different cross-sectional planes within the bend were seen to accurately match the experimentally obtained contour levels, showing that the occurring secondary flow can be precisely captured by the SPH simulations. Another important quantity characterizing the flow in curved pipes is the variation of the wall static pressure along the inner and outer arc of the bend. This provides information about the pressure gradients at different stations on the bend wall. The numerical simulations predicted the experimentally measured wall static pressures with RMSE deviations between ∼0.8 and ∼1.8% and reproduced very well the dependence of the pressure coefficient with position across the bend. The numerical solutions indicate that global convergence to the results of the experimental test bend analyzed by Enayet et al. [1] with the present SPH scheme requires working with a minimum of 16 million particles to correctly reproduce the behavior of resolution-sensitive quantities, such as the turbulence intensity and the variation of the wall static pressure, when coupling the scheme to a LES model. The good agreement demonstrated between the numerical solution and the experimental data for varying Re shows that the present SPH scheme can handle turbulent flows with pronounced streamline curvatures and potential separation regions with an acceptably good degree of prediction.