RANS-Based Modelling of Turbulent Flow in Submarine Pipe Bends: Effect of Computational Mesh and Turbulence Modelling

Abstract

Pipe bend is a critical integral component, widely used in slurry pipeline systems involving various engineering applications, including natural gas hydrate production. The aim of this study is to assess the capability of RANS-based CFD models to capture the main features of the turbulent single-phase flow in pipe bends, in view of the future investigation of the hydrate slurry flow in the same geometry. This is different from the available literature in which only a few accounted for the effects of a combination of computational mesh, turbulence model, and near-wall treatment approach. In this study, three types of mesh configuration were adopted to carry out the computations, namely unstructured mesh and two structured meshes with a uniform and nonuniform inflation layer, respectively. To explore the influence of the turbulence model, standard k-ε, low-Reynolds k-ε, and nonlinear eddy viscosity turbulence model were selected to close RANS equations. Pressure coefficient, mean axial velocity, turbulence intensity, secondary flow velocity, and magnitude of secondary flow were regarded as the critical variables to make a comprehensive sensitivity analysis. Predicted results suggest that turbulent kinetic energy is the most sensitive variable to the computational mesh while others tend to stabilize. The largest difference of turbulence kinetic energy was around 26% between unstructured mesh and structured mesh with a nonuniform inflation layer. Additionally, a fully resolved boundary layer can reduce the sensitivity of mesh on turbulent kinetic energy, especially for a nonlinear turbulence model. However, the large gradient and peak value of turbulence intensity near the inner wall of the bend was not captured by the case with a fully resolved boundary layer, compared with that of the wall function used. Furthermore, it has been confirmed that the same rule was detected also for different curvature ratios, Reynolds numbers, and dimensionless wall distance y+.

1. Introduction

A pipeline system is of paramount value for the transportation of natural gas hydrate slurry (NGHS). A pipe bend is one of the commonly used fitting components in slurry pipeline systems and hydrodynamic devices due to the flexibility of redirecting the flow. As known, even in the condition of dilute flow, excessive wear might happen at the outer wall due to the impingement of moving particles, which in turn results in a risk of leakage. Furthermore, pipe bends also induce excessive head loss related to the solid concentration. Disposing of tools for the prediction of the characteristics of the slurry flow inside the bend, as well as for the prediction of the induced loss of material, allows for the helping of the design and optimization of the slurry system. Computational fluid dynamics (CFD), an analysis of systems involving fluid flow, heat transfer, and associated phenomena such as chemical reactions by means of a computer-based simulation, provides a great opportunity in this sense, mostly since this approach has few limitations with respect to capital, scale, and measurement issues in laboratory tests. Simultaneously, numerical and modelling issues associated with CFD modelling have to be of concern.

To reduce the computational effort in most industrial applications, some averaging operator is usually applied to local instantaneous Navier-Stokes equations yielding the locally averaged quantities of the slurry flow, such as mean velocity and pressure fields, which raises nontrivial doubts regarding the proper type of averaging operator to apply and the modelling of the correlations arising from the averaging process. These aspects are still the subject of discussions among slurry flow modellers, however, no definite answer has been given even for the simplest case of single-phase incompressible flow. This is the consequence, on one hand, of the complexity of the turbulent flow inside a pipe bend and, on the other, of the strong assumptions underlying the Reynolds-Averaged Navier-Stokes (RANS) approach.

However, addressing the numerical simulation of the slurry flow in a certain geometry requires a clear understanding of how to model the same type of flow in the absence of particles. This is mainly due to two reasons. First of all, the single-phase flow can be interpreted as the limiting case of a slurry flow when the solid concentration tends to zero, and, therefore, slurry CFD models must be consistent with the single-phase ones. Second, the modelling of some of the correlations in slurry CFD models is usually performed via two-phase extensions of the commonly used turbulence models for single-phase flows. As a result, the goal of this study is to shed light on the capability of RANS-based CFD models to capture the physical processes of the turbulent single-phase flow in pipe bends, in light of the future investigation of the slurry flow in the same geometry. It should be mentioned that several studies adopted the large eddy simulation (LES) for simulating the turbulent flow in pipe bends with much more computational effort and more accurate results [1,2]. However, the LES approach is beyond the scope of this study since the goal is to provide the theoretical basis for the two-phase flow modeling where the RANS-based modelling is still applied.

Undoubtedly, more or less than half of the time spent on the CFD simulations is devoted to the preprocessing procedure, in particular, geometry building and mesh generation since the number and even shape of the grid can influence the convergence and stability of numerical simulation, finally resulting in an inaccurate solution. Generally, structured mesh composed of a quadrilateral in 2D and a hexahedron in 3D, respectively, is first used to discrete the numerical domain for the simple geometries. Gradually, the flow in complex geometry must be investigated and the unstructured mesh created, including a triangle in 2D, and a tetrahedron, prism, and pyramid in 3D. Despite the flexibility of unstructured mesh applied to a wide range of geometries, the structured mesh is still popular nowadays for the reason that a structured grid can be simply addressed by using the indices (i, j, and k) and the connectivity of a neighbour cell is very clear. Apart from this, governing equations can be solved by means of a straightforward solution algorithm without requiring the implementation of extra terms, provided that the grid is nearly orthogonal and without skewness. Conversely, once the grid is highly nonorthogonal, structured mesh might lead to an unphysical result, as the additional algorithms concerning the nonorthogonal issue have to be added which, in turn, results in numerical instability and divergence. In terms of an unstructured grid, a higher number of cells is needed to achieve the same accuracy, implying that a larger computational effort is required, compared with the structured mesh. In light of a case such as turbulent flow in a pipe bend, the boundary layer near the pipe wall might have a certain effect on the whole domain due to the large variety of variables here. For this reason, boundary layer resolution should also be considered with the application of a finer mesh near the wall to satisfy the log law of the wall or a fully resolved boundary layer up to the viscous sublayer.

In general, fluid particles pass through a curved surface, such as a bend, which is mainly characterized by secondary flow, axial flow separation, and unsteadiness. Here a secondary flow occurs in both the laminar and turbulence states resulting from the imbalance of pressure gradient at the radial direction of the bend cross-section and the centrifugal force produced by the curved geometry, referring to a pair of counter-rotating vortices determined by the Dean number, which is also called the Dean vortex [3]. With this type of pressure-driven secondary flow superimposing on the primary flow, the primary streamwise flow is distorted to a large extent. The faster fluid is transferred from the inner part of the bend through the centerline of the cross-section to the outer side of the bend while the fluid with a lower momentum to the inner section is along the top and bottom wall. In the case of a fully developed pipe bend flow, such as that occurring, for instance, in long coiled pipes, the velocity field will be the same for every cross-section. Conversely, when a bend of finite length is installed within two straight pipes, the characteristics of the primary and secondary flows will change along the pipe centerline. In order to describe such a system, similar to Sudo et al. [4], the most straightforward coordinate system seems to be a combined cylindrical-toroidal one, as demonstrated in Figure 1, in which the axial positions of the cross-sections of the upstream pipe, the bend, and the downstream pipe are identified, respectively, by the distance from the bend inlet section, $z^{\prime }$, the angle from the bend inlet section, $\phi$, and the distance from the bend outlet section, $z$. Then, the points of each cross-section are identified by an azimuthal coordinate, $\theta$, and a radial coordinate r. Therefore, based on the above considerations, the primary and secondary flow patterns are expected to vary with $z^{\prime }$, $\phi$, and $z$. In some situations, the primary flow separates from the intrados of the pipe bend, giving rise to a streamwise vortex, provided that the curvature ratio defined as the proportion of curvature radius and pipe diameter R_c/d is less than 1.5 [5,6]. Clearly, the complex fluid dynamic behavior of pipe bend flows, especially in the turbulent regime, makes obtaining accurate, whole-field experimental data a very challenging task, especially in terms of velocity measurements.

To date, the turbulent flow through a pipe with a 90° or 180° bend has been widely explored by experimental and numerical approaches since it is of paramount importance to industrial practice. A comprehensive literature review on the studies of a curved pipe or channel from a historical point of view can be found in reference [7]. To the authors’ knowledge, the first research on the curved pipe can be dated back to 1855. Weisbach [8] carried out the experiment to measure the pressure loss in a sharp pipe bend. Followed by Weisbach, Thomson [9,10] pointed out that the difference between the flow in a curved pipe and wingding of the river and the secondary flow is induced by the imbalance of the pressure gradient and the centrifugal force. Williams et al. [11] first illustrated that the maximum mean velocity at the bend exit is displaced to the outer part of the pipe. Eustice [12,13] demonstrated that the pressure drop in a pipe bend is larger than that of a straight pipe. It is accepted that in one of the most notable works made by Dean [3,14], one pair of counter-rotating vortices were analytically calculated for laminar flow in a long pipe bend, which is determined by a parameter, namely the Dean number. Due to the limitation of the apparatus, the above-mentioned investigations are almost carried out in the laminar condition. With the development of measurement instruments, laser doppler velocimetry (LDV) and hot-wire anemometry have been applied to measure the mean quantities and turbulence statistics for three-dimensional components. Adler [15] first showed the mean velocity profiles and contours at a range of Reynolds numbers from 2000 to 12,000. Weske [5] conducted a series of experiments, aimed at exploring the velocity distribution along pipe bends for different types of cross-sections, such as circular, elliptical, square, and rectangular shapes. It was suggested that the shape and aspect ratio of the cross-section of the pipe have little effect on the mean velocity compared with the influence of the curvature ratio. Rowe [16] used a Pitot tube and yawmeter to measure the total pressure and yaw for a flow in the U bend and S bend. It was found that the secondary flow in a U bend started from 30° and decreased gradually until 90°. Then it became steady until the bend exit. In comparison to the U bend, the fully developed flow recovered faster downstream of the S bend. In order to further fully understand the profile of the flow, the longitudinal, circumferential, and radial components of mean and fluctuating velocities, and the Reynolds stresses in the pipe cross-section at several longitudinal stations were measured by using the technique of hot wire [4,17], extending the works of Enayet et al. [18] and Azzola et al. [19], who limited their investigations to measuring only the longitudinal and circumferential velocity components.

At the same time, the turbulent flow in pipe bends has been the subject of several numerical investigations, helped by the growing capabilities of computers and multipurpose CFD codes. However, a distinction shall be made according to the turbulence modelling approach and the purpose of the study. In the initial stage, based on the finite difference method, Patankar et al. [20] solved RANS equations with the closure of the standard k-ε turbulence model to predict the turbulent flow in a U bend, referring to the measurements made by Rowe [21]. They claimed that the calculations were not in good agreement with experimental data due to the isotropic property of the eddy viscosity model. Al-Rafai [22] investigated the effect of the curvature ratio on the streamwise velocity and RSM velocity via experimental and numerical techniques. In their simulation, the polar mesh was utilized to discretize the domain based on the body-fitted coordinates employed onto PHOENICS. The numerical results qualitatively agree with the measured data while having a significant quantitative deviation when it comes to the mean axial velocity. Both methods proved that the secondary flow in a sharper bend is more intense. Lai et al. [23] took into account the influence of boundary layer resolution for the Reynolds stress model (RSM) on a simulation of turbulent flow in a U bend. Two high-Reynolds RSM turbulence models, different from the pressure-strain model, with the application of wall function, have worse performance in predicting the third cell, normal stress, and downstream recovery than is the case with the low-Reynolds RSM turbulence model, which indicates that a fully resolved viscous sublayer results in a more accurate solution. Note that the grid-independent analysis for two polar meshes was also performed in their study, however, only the mean velocities profile and the in-plane velocity of 67.5° were considered. Also, Wang and Shirazi [24] utilized polar mesh to predict the characteristics of flow in a 90° bend with high and low-Reynolds k-ε turbulence models. In contrast to Lai et al. [23], the low-Reynolds k-ε model with a fully resolved boundary has a slight improvement in predicting the axial velocity at the outer wall of the bend. In terms of grid influence, the mean axial velocity at 30° for two polar meshes was compared. Note that they only doubled the cell number of the radial direction and claimed the difference between the results from refined mesh and experimental data was less than 1~2%, which does not seem reliable to make grid independence. Pruvost et al. [25] further carried out various simulations to investigate the influence of boundary resolution and turbulence models on the flow profile in a 90° and a 180° bend, respectively, using a structured hexahedron grid in circumferential, and an unstructured hexahedron grid in the core region. The grid convergence was achieved by solely comparing the mean axial velocity and turbulence kinetic energy at the bend exit with different mesh sizes. Similar to Wang and Shirazi [24], it was concluded that there was little improvement for the 90° bend after the fully resolved boundary layer (dimensionless wall distance was around one), regarding the mean axial velocity and static wall pressure. However, the performance of the model with a fully resolved boundary layer has obvious improvement for 180° bend, compared with those cases where the wall function was applied. A structured hexahedron mesh with four corners having bad orthogonality has been used to simulate double bend flow, incorporating standard k-ε and renormalization group (RNG) turbulence models [26]. Both models gave satisfactory results for mean axial velocity and the existence of a vortex, whereas the RNG model performed rather better in predicting tangential velocity than the standard k-ε model. Nevertheless, the case with the RNG model occupied four times the computer resources than that of the standard k-ε model. Kim et al. [27] selected an automatically generated unstructured hexahedron grid to computationally reproduce the experiments of Sudo et al. [4] in order to find an optimum turbulence model and further verify their data. Wang et al. [28] simulated a turbulent flow of hydraulic oil through a 90° bend by using an unstructured mesh combined with a prismatic core and hexahedron cell near the wall. An O-type structured hexahedron grid was applied to simulate the flow in a bend with 90° [29,30]. Afterward, O-type mesh was also used by Hossain [31], aiming at reducing the skewness near the curved surface induced by H-typed mesh. They used an O-type mesh to compare the performance of two turbulence models in predicting the development of flow in a 180° circular channel. The results implied that the realizable k-ε turbulence model was slightly better than the RSM model. Alternatively, an unstructured mesh with a tetrahedral core and prism near the wall was utilized to discrete the bend domain [32,33,34].

Previous studies have made some great contributions to CFD calculation for turbulent flow in the curved pipe, mostly considering the turbulence model and the boundary layer resolution. However, to the authors’ best knowledge, the previous analyses either focused on some features of the flow neglecting others or, finally, provided only little clarification on the effect of numerically related parameters such as the computational mesh. As mentioned before, the computational mesh has a strong influence on the numerical stability and convergence, thereby affecting the accuracy of the result, which allows for determining whether the simulation fails or succeeds to some degree. As a result, the discussion of the mesh effects is of great value for the application of CFD modelling.

Based on these considerations, the effort is being paid to assess the capability of RANS-based CFD models to capture the main features of the turbulent single-phase flow in pipe bends, in view of future investigation of the hydrate slurry flow in the same geometry. Essentially, the goal is to provide a clear overview of the effect of the computational mesh on the predicted features of the flow that either affects the dynamics of solid particles or has a greater impact on the applications. This paper is a preliminary attempt to explore how the computational mesh could affect the characteristics of turbulent flow through a bend using OpenFOAM code in which three types of meshes (an unstructured mesh and two structured meshes), three turbulence models, as well as the issue of boundary layer resolution are all accounted for. The results obtained here can provide a theoretical basis for the solid–liquid flow in the pipe bend.

2. Numerical Methodology

Due to the complexity of turbulent flow through a curved pipe categorized by the secondary flow and separation flow, it is of vital interest to determine an optimum mesh configuration for engineering applications. On the grounds of an open-source code OpenFOAM via the CINECA platform, the cell-centred finite volume approach was applied to solve the three-dimensional steady-state RANS equations, with a pressure-velocity coupling achieved using a SIMPLE algorithm. The second-order central-differencing scheme was applied to discretize all terms in equations. There is no limitation on the time step since the simulation is typically time independent. Considering the convergence and stability of the modelling, the under-relaxation factors were accounted for. Except for the pressure set as 0.3, all of the variables were 0.7. The OpenFOAM solver was terminated after reaching the convergence criterion which was 1 × 10⁻⁵ and 1 × 10⁻⁶ for pressure and other variables, respectively.

2.1. Governing Equations

Similar to previous works in which the RANS equations were solved for predicting the turbulent flow in pipe bends [35,36,37], the steady-state RANS equations for incompressible fluid flow are as below. Note that the physical meaning of each variable residing in all following equations is omitted for the sake of space.

$$\frac{\partial U_i}{\partial x_i}=0$$

(1)

$$\frac{\partial (U_i{U}_j)}{\partial x_j}=-\frac{1}{\mathsf{\rho }}\frac{\partial P}{\partial x_i}+\frac{\partial }{\partial x_j}(\mathsf{\nu }\frac{\partial U_i}{\partial x_j}-\overline{u_i^{’}{u}_j^{’}})+g$$

(2)

2.2. Turbulence Model

As for the straight pipe flow, the standard k-ε model has been proven to be accurate enough to model the turbulence encountered in industrial applications. However, for the pipe bend associated with secondary flow, the nonlinear cubic k-ε model is preferable due to the pressure-strain and curvature effects considered in the Reynolds stress formulation ($\overline{u_i^{’}{u}_j^{’}}$), according to the suggestions of previous explorations. Here, three turbulence models with different accuracy levels were initially used for exploring the effect of the turbulence model, namely the standard k-ε model, low-Reynolds k-ε model (LaunderSharmaKE), and nonlinear eddy viscosity model (LienCubicKE). The transport equations of each turbulence model are as following:

2.2.1. Standard k-ε Turbulence Model

The standard k-ε turbulence model of Launder and Spalding [38] is one of the widely used models in various engineering fields due to its simplicity and robust property. It is worth noting that this model is not suitable for a flow with curved streamlines and a strong pressure gradient.

$$\frac{\mathrm{D}k}{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}]+P_k-\mathsf{\epsilon }$$

(3)

$$\frac{\mathrm{D}\epsilon }{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }})\frac{\partial \epsilon }{\partial x_j}]+C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}$$

(4)

where the related closure terms and coefficients of Equations (3) and (4) are summarized in

Table 1. Related terms and coefficient in standard k-ε turbulence model.

$\overline{{\mathit{u}}_{\mathit{i}}^{’}{\mathit{u}}_{\mathit{j}}^{’}}/\mathit{k}$	${\mathit{\nu }}_{\mathit{t}}$	${\mathit{P}}_{\mathit{k}}$	S	${\mathit{S}}_{\mathit{i}\mathit{j}}$
$\frac{2}{3}{\delta }_{ij}-\frac{{\nu }_t}{k}{S}_{ij}$	$C_{\mu }\frac{k^2}{\epsilon }$	${\nu }_t{S}^2$	$\sqrt{2S_{ij}{S}_{ij}}$	$\frac{1}{2}(U_{i,j}+U_{j,i})$
$C_{\epsilon 1}$	$C_{\epsilon 2}$	$C_{\mu }$	${\sigma }_k$	${\sigma }_{\epsilon }$
1.44	1.92	0.09	1	1.3

.

2.2.2. Low-Reynolds k-ε Turbulence Model

One of the low-Reynolds k-ε turbulence models implemented in OpenFOAM was developed by Launder and Sharma [39], which was chosen in this work.

$$\frac{\mathrm{D}k}{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}]+P_k-\tilde{\epsilon }-D$$

(5)

$$\frac{\mathrm{D}\tilde{\epsilon }}{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_{\tilde{\epsilon }}})\frac{\partial \tilde{\epsilon }}{\partial x_j}]+C_{\epsilon 1}{f}_1\frac{\tilde{\epsilon }}{k}{P}_k-C_{\epsilon 2}{f}_2\frac{{\tilde{\epsilon }}^2}{k}+E$$

(6)

where the related closure terms and coefficients of Equations (5) and (6) are given in

Table 2. Related terms and coefficient in low-Reynolds k-ε turbulence model.

${\mathit{\nu }}_{\mathit{t}}$	$\tilde{\mathit{\epsilon }}$	$\mathit{D}$	$\mathit{E}$	${\mathit{f}}_{\mathit{\mu }}$	${\mathit{f}}_1$
$C_{\mu }{f}_{\mu }\frac{k^2}{\tilde{\epsilon }}$	$\epsilon -D$	$2\nu {(\frac{\partial \sqrt{k}}{\partial x_i})}^2$	$2\nu {\nu }_t{(\frac{{\partial }^2{U}_i}{\partial x_i\partial x_k})}^2$	$e^{\frac{-3.4}{{(1+R_T/50)}^2}}$	1
$f_2$	$R_T$	$C_{\epsilon 1}$	$C_{\epsilon 2}$	$C_{\mu }$	${\sigma }_k$	${\sigma }_{\epsilon }$
$1-0.3e^{-R_T^2}$	$\frac{k^2}{\nu \tilde{\epsilon }}$	1.44	0.92	0.09	1	1.22

.

2.2.3. Nonlinear Viscosity Turbulence Model

The LienCubicKE model proposed by Lien et al. [40] in OpenFOAM can be applied as both the high-Reynolds k-ε model with the usage of wall function and the low-Reynolds k-ε model with the fully resolved boundary layer. The model includes third-order terms to account for the streamline curvature effects.

$$\frac{\mathrm{D}k}{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_k})\frac{\partial k}{\partial x_j}]+P_k-\epsilon$$

(7)

$$\frac{\mathrm{D}\epsilon }{\mathrm{D}t}=\frac{\partial }{\partial x_j}[(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }})\frac{\partial \epsilon }{\partial x_j}]+C_{\epsilon 1}{f}_1\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}{f}_2\frac{{\epsilon }^2}{k}+E$$

(8)

Where the related closure terms and coefficients of Equations (7) and (8) are listed in

Table 3. Related terms and coefficient in LienCubicKE turbulence model.

${\mathit{\nu }}_{\mathit{t}}$	${\mathit{f}}_{\mathit{\mu }}$	${\mathit{f}}_1$		${\mathit{f}}_2$	${\mathit{y}}^{*}$	${\mathit{R}}_{\mathit{T}}$	${\mathit{P}}_{\mathit{k}}^{’}$
$C_{\mu }{f}_{\mu }\frac{k^2}{\epsilon }$	$(1-e^{-0.198y^{*}})$($1+\frac{5.29}{y^{*}}$)	$1+\frac{P_k^{’}}{P_k}$		$1-0.3e^{-R_T^2}$	$y\sqrt{k}/\nu$	$\frac{k^2}{\nu \epsilon }$	$1.33(1-0.3e^{-R_T^2})(1+2\nu \frac{k}{y^2})e^{-0.00375y^{*2}}$
$\overline{u_i^{’}{u}_j^{’}}/k$
$\frac{2}{3}{\delta }_{ij}-\frac{{\nu }_t}{k}{S}_{ij}$+$C_1\frac{{\nu }_t}{\epsilon }(S_{ik}{S}_{kj}-\frac{1}{3}{\delta }_{ij}{S}_{kl}{S}_{kl})+C_2\frac{{\nu }_t}{\epsilon }({\Omega }_{ik}{S}_{kj}+{\Omega }_{jk}{S}_{ki})+C_3\frac{{\nu }_t}{\epsilon }({\Omega }_{ik}{\Omega }_{jk}-\frac{1}{3}{\delta }_{ij}{\Omega }_{lk}{\Omega }_{lk})+C_4\frac{{\nu }_{t}k}{{\epsilon }^2}(S_{ki}{\Omega }_{lj}+S_{kj}{\Omega }_{li})S_{kl}+C_5\frac{{\nu }_{t}k}{{\epsilon }^2}(S_{kl}{S}_{kl}-{\Omega }_{kl}{\Omega }_{kl})S_{ij}$
$C_{\mu }$	$C_1$	$C_2$	$C_3$	$C_4$	$C_5$	$S$	$\Omega$	${\Omega }_{ij}$
$\frac{2/3}{1.25+S+0.9\Omega }$	$\frac{3/4}{1000+S^3}$	$\frac{15/4}{1000+S^3}$	$\frac{19/4}{1000+S^3}$	$-10C_{\mu }^2$	$-2C_{\mu }^2$	$\frac{k}{\epsilon }\sqrt{\frac{1}{2}{S}_{ij}{S}_{ij}}$	$\frac{k}{\epsilon }\sqrt{\frac{1}{2}{\Omega }_{ij}{\Omega }_{ij}}$	$\frac{1}{2}[U_{i,j}-U_{j,i}]$

.

2.3. Computational Domain and Mesh

In this study, the experimental setup (Figure 1) for turbulent airflow through a 90° bend made by Sudo et al. [4] was used in simulations. The bend had an inner diameter d = 2R = 104 mm and a curvature radius of R_c = 208 mm. In this study, the computational domain was long enough to achieve a fully-developed flow condition, including a straight pipe in the upstream (30d), bend and a straight pipe in the downstream (20d). Two different meshes produced by ANSYS Workbench and one created by blockMesh were utilized to discretize the whole domain. Figure 2a shows the schematic of the straight pipe with a 90° bend. In the upstream and downstream mesh named in part one, a larger cell size than that of the bend (part 2) was used to decrease the total number of elements which in turn reduced the computational effort, as shown in Figure 2. Considering the influence of mesh and inflation layer on the flow field, the detailed sensitivity analysis of meshes was presented hereafter.

2.4. Boundary Conditions

Typically, three types of boundary conditions were required to be applied in pipe flow simulations. First, the uniform velocity, the turbulent kinetic energy, and the turbulent dissipation rate were prescribed at the pipe inlet (k = k_in, ε = ε_in, U = 0, V = 0, and W = W_in). At the pipe wall, the nonslip condition was imposed (U = V = W = 0). It should be noted that the log-law wall function had been applied for the high-Reynolds turbulence model while the boundary layer was fully resolved for capturing the viscous sublayer for the low-Reynolds turbulence model. Once the wall function was used, the dimensionless wall distance y+ from the first near-wall cell was strictly dropped in the range from 30 to about 130 in such a case, whereas a fully resolved boundary layer needs y+ of 1, according to the theory of channel flow. Finally, similar to references [35,36,37], the normal gradient of any variable has been considered as zero at the pipe outlet ($\partial \psi /\partial n=0$), except for a reference pressure value which was set as zero (relative pressure, P_r = 0).

3. Effect of Computational Mesh and Turbulence Modeling

3.1. Parameters Used in Benchmarking

Here the preliminary study was obtained by numerically simulating the experiments of [4] since comprehensive measurements of turbulent flow through a bend with 90° have been made. The parameters used in their experiments were as follows: the bend had an inner diameter of d = 2R = 104 mm and a curvature radius of R_c = 208 mm, and the carrier fluid was air flowing (density ρ = 1.22 kg/m³ and kinematic viscosity ν = 1.46 × 10⁻⁵ m²/s) at a bulk-mean velocity of W_in = 8.7 m/s. Therefore, the bulk Reynolds number was 6 × 10⁴, and the curvature ratio, R_c/d, was equal to two. The hot-wire anemometer technique was applied to measure the mean and fluctuating velocity fields in three directions, namely the longitudinal, circumferential, and radial directions. Reynolds stresses in several longitudinal cross-sections were also obtained. Additionally, the pressure coefficient, turbulent kinetic energy, secondary flow intensity, and primary flow deviation along the pipe were presented.

3.2. Types of Computational Meshes

As shown in Figure 2, the dimensionless distance from pipe wall y+ = 30 and y+ = 1 were both generated for three types of meshes, which depends on the choice of turbulence model. The standard k-ε turbulence model is accurate only if the first grid resides in the log-law region, here corresponding to y+ = 30. However, the low-Reynolds turbulence model, such as LaunderSharmaKE in OpenFOAM, requires a fully resolved boundary layer with y+ of 1. Note that a nonlinear eddy viscosity turbulence model such as LienCubicKE implemented in the current version of OpenFOAM is combined with its low-Reynolds and high-Reynolds models. Therefore, both meshes with y+ = 30 and y+ = 1 were used to test the mesh sensitivity. The inflation layer was applied to determine how the boundary layer was resolved. In terms of grid-independence analysis, the refinement of mesh was made by progressively reducing the cell size of two parts marked in Figure 2a. In the ANYSYS workbench, the body size tool was applied to define the cell size of those two parts for unstructured mesh and structured mesh with the nonuniform inflation layer while the number of divisions was modified in blockMesh for structured mesh with the analogous uniform inflation layer. All of the cases run in this study with detailed information related to mesh types, y+, cell number, and mesh quality, are listed in

Table 4. The details of computational meshes used in this study.

Mesh Types	Abbreviation	The Cell Size of Part 1 (m)	The Cell Size of Part 2 (m)	Number of Inflation Layer		Growth Ratio of Inflation Layer		Total Cells Number		Average Nonorthogonality		Max Skewness
				y+ = 30	y+ = 1	y+ = 30	y+ = 1	y+ = 30	y+ = 1	y+ = 30	y+ = 1	y+ = 30	y+ = 1
Unstructured mesh	UM1	0.0095	0.007	5	14	1.05	1.325	713,135	1,260,911	12.46	9.54	0.745	0.756
	UM2	0.008	0.006					1,093,009	1,843,736	12.74	9.83	0.736	0.791
	UM3	0.0065	0.005					1,843,134	2,935,319	12.99	10.22	0.862	0.780
	UM4	0.0055	0.004					3,054,005	-	13.19	-	0.803	-
	UM5	0.004	0.003					5,066,295	-	13.59	-	0.843	-
Structured mesh with uniform inflation layer	SUM1	0.0095	0.007					163,584	347,616	12.91	9.43	0.529	0.577
	SUM2	0.008	0.006					241,808	497,840	13.38	9.67	0.594	0.64
	SUM3	0.0065	0.005					377,400	744,600	13.98	10.13	0.67	0.711
Structured mesh with nonuniform inflation layer	SNUM1	0.0095	0.007					246,144	553,824	5.81	4.42	1.175	0.601
	SNUM2	0.008	0.006					360,808	785,288	5.98	4.59	1.020	0.630
	SNUM3	0.0065	0.005					579,309	1,205,589	6.20	4.81	0.837	0.715
	SNUM4	0.0055	0.004					953,988	-	6.43	-	0.670	-
	SNUM5	0.004	0.003					1,761,984	-	6.69	-	0.605	-

. To keep the similar inflation layer among the three types of meshes, only the body size of the core region was changed, and the inflation layer was kept at an unchanged number of layers and expansion ratio. It should be mentioned that only three meshes were tested for structured mesh with a uniform inflation layer for the reason that this type of mesh converged quickly with respect to the effect of grid size. Additionally, a nonuniformed mesh size was used for part one of this mesh type, thereby the total cell number is relatively smaller than that of the other structured type. But it is ensured that the grid size in part two, particularly near the bend, is nearly the same as other mesh types.

3.3. Grid Sensitivity Analysis

Among all cases, the pressure coefficient, turbulence kinetic energy (Equation (9)), and secondary flow magnitude (Equation (10)) along the pipe from upstream to downstream of the bend, axial velocity, in-plane velocity, and turbulence intensity at the cross-section of the bend exit were regarded as the identifier to verify if the solution reached grid-independent since those variables were used to characterize such a flow, according to a critical literature review. All variables were normalized by a factor according to Sudo et al. [4]. For the sake of brevity, only some representative results are presented in this paper. It should be noted that although some variables are shown in a qualitative manner, the results can still provide valuable information.

The grid convergence tests were carried out by considering three types of typical meshes, as illustrated in Table 4. UM1 to UM5, SUM1 to SUM3, and SNUM1 to SNUM5

$$k_a=\frac{4}{\pi d^2{W}_{in}^2}{{\displaystyle \int }}_{-\pi /2}^{\pi /2}{{\displaystyle \int }}_0^R\frac{1}{2}(\overline{u^{’2}}+\overline{v^{’2}}+\overline{w^{’2}})r\mathrm{d}r\mathrm{d}\theta$$

(9)

$$I_s=\frac{4}{\pi d^2{W}_{in}^2}{{\displaystyle \int }}_{-\pi /2}^{\pi /2}{{\displaystyle \int }}_0^R\frac{1}{2}(U^2+V^2)r\mathrm{d}r\mathrm{d}\theta$$

(10)

denote the different mesh sizes’ composition, ranging from coarse to fine for unstructured, structured with uniform inflation layer, and structured mesh with nonuniform mesh, respectively. During the grid refinement process, the total cell number was increased from coarse to fine grid as most researchers did. As anticipated, the turbulence models can lead to some differences among the results. Specifically, the converging rate changed with turbulence models for different mesh types. For instance, as shown in the top graph of Figure 3, all of the variables converged quickly, however, the turbulence kinetic energy attained convergence only in very fine meshes (UM5 and SNUM5), referring to the application of the nonlinear turbulence model with y+ of 30. Clearly, the secondary flow appears at the bend inlet, reaches the maximum value at the vicinity of the bend exit, and its magnitude starts to reduce rapidly. Similarly, the turbulence intensity begins to increase, arrives at its highest point in the vicinity of the bend exit, and then decreases. It is worth noting that there is no need to compare the numerical results with experimental data in the current stage (sensitivity analysis), aiming at exploring the effect of numerical and modeling factors on the numerical solutions while the validation of the model definitely requires it. The possible reason for the differences in the numerical results and experimental measurements could be the application of RANS-based turbulence modeling, which is incapable of accurately capturing the turbulence field in such a complex phenomenon, however, the predicted results should satisfy most of the engineering applications. Also, a similar observation was detected for the contour plot of turbulence intensity (Figure 4) when the standard k-ε turbulence model was used, however, this phenomenon only appeared in the unstructured mesh. Clearly, the largest local fluctuation (12%) appeared in the region near the inner wall due to the steep velocity gradient (Figure 5) along with the depression of the mean velocity, whereas a relatively weak fluctuation happened to the outer wall since the velocity gradient there was gentle. However, as demonstrated in Figure 5, the mean axial velocity at the bend exit was slightly sensitive to the mesh sizes, whereas the changes in the other variables were negligible. It should be mentioned that the left-hand side was the inner wall of the bend while the right represented the outer wall in all of the contour and vector plots in this study. Obviously, the tongue-shaped velocity distribution appeared due to the remarkable distortion of primary fluid induced by the secondary flow near the bend exit. It could be also detected from the three-dimensional flow structure by streamlining, as shown in Figure 6. It was clear that the convergence rate was not significantly affected by the mesh sizes once the boundary layer was fully resolved, corresponding to y+ of 1. UM3, SUM3, and SNUM3 were fine enough to get a grid-independent solution for all turbulence models considered in this work, which was probably a consequence of a more accurate solution obtained from the resolved boundary layer than that with the usage of the wall function.

In terms of other variables, such as pressure coefficients (C_p = (P_a−P_ref)/0.5ρW_in²), indicating the difference of wall static pressure in various locations and a reference point in the upstream (z = −17.6d) divided by 0.5ρW_in² and the velocity vector of the secondary flow at bend exit are demonstrated in Figure 7 and Figure 8, respectively. For the sake of space, the results for those two variables calculated by the nonlinear turbulence model were taken as an example. As for the pressure coefficients (Figure 7), overall, no difference was observed, neither from the different mesh sizes nor the mesh types. However, the pressure drop predicted by two structured meshes was slightly larger than that of unstructured mesh. Additionally, all pressure drop calculated by the low-Reynolds turbulence model was larger than that of the high-Reynolds models. This is also true for the nonlinear turbulence model with y+ of 1. Due to the insensitivity to mesh sizes, a secondary velocity vector was only presented for the finest mesh. Obviously, all of the mesh types could capture well the Dean vortex (Figure 8) which should be a pair of counter-rotating vortices caused by the symmetric pipe. Usually, the secondary flow induced by the centrifugal force occurs in the cross section of the bend starting from φ = 30° and the Dean vortex circulates outwards in the center part of the pipe and inwards near the upper and lower walls [20]. Not surprisingly, the Dean vortex predicted by SNUM5 was relatively weaker especially in the four corners near the inflation layer due to the certain skewness in this region. This difference can be neglected.

The types of meshes still have some influence on the turbulence intensity when a nonlinear turbulence model with a wall function (y+ = 30) was applied. Figure 9 demonstrates the turbulence intensity at the bend exit for the three types of meshes. Again, the variable shown here is grid independent. It was obvious that a large deviation happens to the inner wall of the bend among the different mesh types. A relatively larger gradient of turbulence intensity was captured by UM5 and SNUM5, and a wide range of larger isolines was found in the simulation SNUM5. This phenomenon was consistent with the average turbulence kinetic energy illustrated in Figure 3, in the sense that the average turbulence kinetic energy of the cross-section for SNUM5 was larger than that of the other two meshes. It was also found that the variable was almost the same between SUM3 and UM5 due to the wide range of larger isolines in SUM3, although there was no strong gradient of turbulence kinetic energy detected. This was also true for the effect of mesh types on the axial velocity at the bend exit, as demonstrated in Figure 10. As can be seen from Figure 11, the shape of the isolines near the inner wall predicted by SNUM5 was different from the other two to some degree. There were also some spurious types of mesh effects with the standard k-ε model, however, this is more moderate.

Surprisingly, as can be seen from Figure 9, the performance of a nonlinear turbulence model with a fully resolved boundary layer (y+ = 1) was not improved as expected, regarding the turbulence intensity. Its peak value was lower than that obtained from the model with the application of the wall function and no large gradient occurred near the inner wall. So far, the reason for this observation is still not yet clear.

It is accepted that the adoption of a log-law wall function ranged from y+ of 30 to 300, more or less. Given that, y+ of 40 and 50 were also considered in this study to confirm whether the observations detected above will be changed with a different y+. Here, only the predictions of turbulence intensity made by structured mesh with a uniform inflation layer and a nonlinear turbulence model with y+ of 40 and 50 are presented (Figure 11). The results demonstrated that the same rule was also suitable for the larger y+ than 30.

3.4. Considerations of Computational Effort

The mesh had a significant influence on the stability and convergence of CFD modelling, thereby probably affecting the accuracy of the solutions. As described in Table 4, two peculiarities were used to evaluate the quality of meshes, namely average nonorthogonality and maximum skewness. It seems that SUM and UM meshes are nonorthogonal to a large extent and they have a lower skewness, compared with SNUM. The value was nearly double for average nonorthogonality. However, the overall mesh quality of all of the meshes used here was fairly satisfactory without confronting any convergence issues during calculation. Consequently, the effect of mesh quality can be neglected in this work.

In terms of computational effort, of course, many more cells were required for unstructured mesh, in comparison with structured mesh. Obviously, the total number of unstructured mesh was more or less five times that of structured mesh. Therefore, the calculation time of the unstructured mesh needed was rather longer than the structured mesh. Even though the half geometry of the pipe bend allows it to be used, the computational effort was still huge, aiming at getting grid-independent solutions.

3.5. Influence of Curvature Ratio and Reynolds Number

The curvature ratio and Reynolds number are the crucial parameters for turbulent pipe bend flow. As is known, a longer bend and a larger Reynolds number can alleviate the effect of bend in general. Furthermore, this influence induced by the curvature ratio was much stronger than the Reynolds number. In order to get rid of the effects of the Reynolds number and the curvature ratio, another two curvature ratios (R_c/d = 1 and R_c/d = 3) and the Reynolds number (Re = 35,616 and Re = 88,328) were chosen for conducting the mesh sensitivity analysis. As expected, it was found that the rule is also valid for these situations no matter what the values of curvature ratio and Reynolds number were. Again, only the turbulence intensity contour at the bend exit for the nonlinear turbulence model and the unstructured mesh with a uniform inflation layer is presented. As illustrated in Figure 12, the large gradient can be also observed near the inner wall when y+ is equal to one due to the large adverse pressure gradient and flow separation. However, the peak value for y+ of 30 was also larger than that of y+ of 1, which was consistent with the previous predictions. The turbulence intensity for a longer bend regarding the different boundary layer treatments is depicted in Figure 13. Clearly, the turbulence intensity decreased significantly with the increases in curvature ratio since the bend effect was relatively weak. As shown in Figure 14 and Figure 15, the large peak value and gradient of turbulence intensity in the vicinity of the inner wall were only captured by the model with the wall function used. Of course, a much larger turbulence intensity appeared in the condition of higher Re than that of the lower one.

4. Conclusions

In this study, RANS-based simulations of the turbulent flow in pipe bends were performed using the OpenFOAM code. The experiment of Sudo et al. [4] was simulated to investigate how the CFD solution was affected by computational mesh type and near-wall mesh resolution, as well as the turbulence model. Three kinds of meshes and turbulence models were taken into account in CFD modelling. Based on the current investigations, the following preliminary conclusions could be drawn:

(1)

Computational mesh has a certain impact on those critical variables concerned in this study. The most sensitive variable regarding turbulent pipe bend flow was turbulent kinetic energy. The sensitivity of the turbulence intensity to mesh sizes also relies on the application of the turbulence model (Figure 3 and Figure 4). The sensitivity of mesh sizes on the turbulence kinetic energy can be reduced, provided that the boundary layer is fully resolved, especially for the configurations where the turbulence model based on the nonlinear eddy viscosity hypothesis is utilized, combined with unstructured mesh and structured mesh with a nonuniform inflation layer.

(2)

In contrast, other variables, such as pressure coefficient (Figure 7), secondary flow velocity, and its intensity (Figure 3) tend to be stable with the refinement of computational mesh, regardless of a slight change of axial velocity with mesh sizes (Figure 5).

(3)

The influence of mesh types on turbulence intensity is remarkable (Figure 9), as the smaller peak value and the gradient close to the near wall were captured by structured mesh with a uniform inflation layer when it comes to y+ of 30, compared with another two meshes. However, the fully resolved boundary layer seemed not to improve the results since the gradient and peak value near the inner wall had not been detected. The reason for this will be investigated in future work.

(4)

Contrarily, the pressure coefficient (Figure 7), secondary flow velocity vector (Figure 8), and its intensity (Figure 3) remained unchanged no matter which mesh type was applied, regardless of visible change in axial velocity with mesh sizes (Figure 10).

(5)

The same rule was still true for the different y+ (Figure 11), curvature ratios (Figure 12 and Figure 13), and Reynolds numbers (Figure 14 and Figure 15). Even though the flow separation occurred in a sharp bend, a similar phenomenon was observed.

Altogether, a recommendation for the computational mesh of turbulent flow through pipe bends could be a structured mesh associated with uniform inflation layers. These findings pave the way for further exploring the hydrate slurry flow in a pipe bend.