Keywords

1 Introduction

The phenomenon that adding a little amount of additives into the turbulent flow would induce a significant reduction of turbulent skin frictional drag was called turbulent DR technology [1]. It owns the merits of remarkable DR effect, relative low cost and easy operation, and has been demonstrated to be of great potential in energy saving within the long-distance liquid transportation and circulation systems. To better apply this technique, the turbulent DR mechanism should be investigated intensively. In recent two decades, the rapid development of computer technology brought great changes to the studies on turbulent flow, where numerical simulation has become an indispensable approach. Among the commonly-used numerical approaches for simulation of turbulent flow, the computational workload of LES is smaller than that of direct numerical simulation (DNS) while the obtained information is much more comprehensive and detailed than that of Reynolds-averaged Navier-Stokes (RANS) simulation. It makes a more detailed but time-saving investigation on the turbulent flow with relatively high Reynolds number possible. Therefore, LES has a brilliant prospective in the research of turbulent flow especially with viscoelastic fluids.

Different from the LES of Newtonian fluid, a constitutive model describing the relationship between elastic stress and deformation tensor should be established first for the LES of viscoelastic fluid. An accurate constitutive model that matches the physical meaning is a critical guarantee to the reliability of LES. Up to now most of the constitutive models for viscoelastic fluid, such as generalized Maxwell model [2], Oldroyd-B model [3, 4], Giesekus model [5] and FENE-P model [6], were all built on polymer solutions because the long-chain polymer is the most extensively used DR additive in engineering application. Compared with Maxwell and Oldroyd-B models, the Giesekus model as well as FENE-P model can characterize the shear thinning behavior of polymer solutions much better and thus they are adopted by most researchers. However, it is reported in some studies that the apparent viscosity calculated with Giesekus model and FENE-P model still deviated from experimental data to some extent, indicating an unsatisfactory accuracy [7]. Inspired by the fact that constitutive model with multiple relaxation times can better describe the relaxation-deformation of microstructures formed in viscoelastic fluid, an N-parallel FENE-P model based on multiple relaxation times was proposed in our recent work [8]. The comparison with experimental data shows the N-parallel FENE-P model can further improve the computational accuracy of rheological characteristics, such as apparent viscosity and first normal stress difference.

For the LES of viscoelastic fluid, the effective SGS model is still absent. To the best of our knowledge, only a handful of researches have studied this problem. In 2010, Thais et al. [9] first adopted temporal approximate deconvolution model (TADM) to perform LES on the turbulent channel flow of polymer solutions, in which the filtered governing equation of LES was derived. Wang et al. [10] made a comparison between approximate deconvolution model (ADM) and TADM. It was found that the TADM was more suitable to LES study on viscoelastic turbulent channel flow. In 2015, comprehensively considering the characteristics of both momentum equations and constitutive equations of viscoelastic fluid, Li et al. [11] put forward a mixed SGS model called MCT based on coherent-structure Smagorinsky model (CSM) [12] and TADM. The forced isotropic turbulent flow of polymer solution and turbulent channel flow of surfactant solution were both simulated by using MCT. The calculation results such as turbulent energy spectrum, two-point spanwise correlations, and vortex tube structures agreed well with the DNS database. Although MCT made an important step to the mixed SGS model coupling spatial filtering and temporal filtering altogether, it was still limited by its spatial SGS model—CSM, with which an excessive energy dissipation is observed near the channel wall. In our recent work [13] we improved the CSM and proposed an improved mixed SGS model named MICT based on ICSM and TADM. Simulation results, for instance, DR rate, streamwise mean velocity, two-point spanwise correlations, were compared to demonstrate a better accuracy of MICT than conventional SGS models.

Regarding the turbulent flow itself, due to the intrinsic complexity and limitation of numerical approaches, together with the complicated non-Newtonian behaviors of viscoelastic fluid, the understanding of turbulent DR mechanism is generally limited within flows with relatively low Reynolds number. For the high Reynolds viscoelastic turbulent flows which are commonly seen in industry and engineering, the turbulent DR mechanism is still far from fully understood, which makes the quantitative guidance on engineering application of DR technology inadequate. Therefore, deeper studies need to be performed on the high Reynolds turbulent drag-reducing flows.

With above background and based on our recent studies on the constitutive model and SGS model of viscoelastic fluid, we apply the LES to further explore the turbulent DR mechanism of high Reynolds viscoelastic turbulent flows in this study. The remainder of this paper is organized as follows. In Sect. 2, the N-parallel FENE-P model based on multiple relaxation times and the improved mixed SGS model MICT are briefly introduced. Then the governing equation of LES for viscoelastic turbulent drag-reducing channel flow is presented accordingly. Section 3 of this paper is devoted to describing the LES approach adopted in this paper. In Sect. 4, the flow characteristics and DR mechanism of high Reynolds turbulent drag-reducing channel flow are deeply studied and analyzed. In the final Section the most important findings of this study are summarized.

2 Governing Equation of LES for Viscoelastic Turbulent Drag-Reducing Flows

2.1 The N-Parallel FENE-P Constitutive Model

To improve the computational accuracy of constitutive model needed for LES of viscoelastic turbulent flows, an N-parallel FENE-P model based on multiple relaxation times was proposed in our recent work [8]. As shown in Fig. 1, the core idea of the proposed constitutive model is to connect N FENE-P models, which has single relaxation time, in parallel. With this parallel FENE-P model, stress-strain relationships of different microstructures formed in viscoelastic fluid are more truly modelled compared with the conventional FENE-P model, and it can better characterize the anisotropy of the relaxation-deformation in viscoelastic fluid. Comparative results indicate the N-parallel FENE-P constitutive model can better satisfy the relaxation-deformation process of microstructures, and characterize the rheological behaviors of viscoelastic fluid with much higher accuracy.

Fig. 1.
figure 1

Schematic of the N-parallel FENE-P model.

Full size image

2.2 The Improved Mixed SGS Model

In the LES of viscoelastic turbulent drag-reducing flow, different from that of Newtonian fluid, filtering of constitutive equation is needed. However, simulation results show the spatial filtering and spatial SGS model that applicable to Newtonian fluid are not so suitable to constitutive equation. Therefore, an improved mixed SGS model named MICT, proposed in our recent work [13] is adopted for the LES in this study. Figure 2 displays the core idea of MICT, that is, ICSM is employed to perform spatial filtering for the continuity and momentum equations within physical space, and TADM is applied to perform temporal filtering for the constitutive equation within time-domain. It has been proved the MICT has much higher computational accuracy in comparison with the MCT SGS model.

Fig. 2.
figure 2

Core idea of the MICT SGS model.

Full size image

2.3 LES Governing Equation of Viscoelastic Turbulent Channel Flow

The fully developed turbulent drag-reducing channel flow of viscoelastic fluid [15] is studied in this paper. Sketch map of the computational domain is shown in Fig. 3, where x, y, z represent the streamwise, wall-normal and spanwise directions respectively, the corresponding size of the plane channel is 10h × 5h × 2h, where h denotes the half-height of the plane channel.

Fig. 3.
figure 3

Sketch map of the computational domain for turbulent channel flow.

Full size image

Based on the N-parallel FENE-P model, the dimensionless governing equation of turbulent channel flow with viscoelastic fluids in Cartesian coordinate system reads

$$ \frac{{\partial u_{i}^{ + } }}{{\partial x_{i}^{*} }} = 0 $$
(1)
$$ \frac{{\partial u_{i}^{ + } }}{{\partial t^{*} }} + u_{j}^{ + } \frac{{\partial u_{i}^{ + } }}{{\partial x_{j}^{*} }} = \delta_{1i} - \frac{{\partial p^{ + '} }}{{\partial x_{i}^{*} }} + \frac{1}{{Re_{\tau } }}\frac{\partial }{{\partial x_{j}^{*} }}\left( {\frac{{\partial u_{i}^{ + } }}{{\partial x_{j}^{*} }}} \right) + \sum\limits_{m = 1}^{N} {\frac{{\beta_{m} }}{{We_{\tau ,m} }}\frac{{\partial \left[ {f\left( {r_{m} } \right)c_{ij,m}^{ + } } \right]}}{{\partial x_{j}^{*} }}} $$
(2)
$$ \frac{{\partial c_{ij,m}^{ + } }}{{\partial t^{*} }} + u_{k}^{ + } \frac{{\partial c_{ij,m}^{ + } }}{{\partial x_{k}^{*} }} - \frac{{\partial u_{i}^{ + } }}{{\partial x_{k}^{*} }}c_{kj,m}^{ + } - \frac{{\partial u_{j}^{ + } }}{{\partial x_{k}^{*} }}c_{ik,m}^{ + } { = }\frac{{Re_{\tau } }}{{We_{\tau ,m} }}\left[ {\delta_{ij,m} - f\left( {r_{m} } \right)c_{ij,m}^{ + } } \right] $$
(3)

where superscript ‘+’ represents the nondimensionalization, \( x_{i}^{*} = {{x_{i} } \mathord{\left/ {\vphantom {{x_{i} } h}} \right. \kern-0pt} h} \), \( t^{*} = {t \mathord{\left/ {\vphantom {t {\left( {h/u_{\tau } } \right)}}} \right. \kern-0pt} {\left( {h/u_{\tau } } \right)}} \), \( u_{i}^{ + } = {{u_{i} } \mathord{\left/ {\vphantom {{u_{i} } {u_{\tau } }}} \right. \kern-0pt} {u_{\tau } }} \), \( p^{ + } = {p \mathord{\left/ {\vphantom {p {\left( {\rho u_{\tau }^{2} } \right)}}} \right. \kern-0pt} {\left( {\rho u_{\tau }^{2} } \right)}} = \overline{{p^{ + } }} + p^{{ +^{\prime } }} \), \( {{\partial \overline{p} } \mathord{\left/ {\vphantom {{\partial \overline{p} } {\partial x_{i} }}} \right. \kern-0pt} {\partial x_{i} }} = - {{\rho u_{\tau }^{2} \delta_{1i} } \mathord{\left/ {\vphantom {{\rho u_{\tau }^{2} \delta_{1i} } h}} \right. \kern-0pt} h} \), \( u_{\tau } { = }\sqrt {{{\tau_{w} } \mathord{\left/ {\vphantom {{\tau_{w} } \rho }} \right. \kern-0pt} \rho }} \); \( u_{i}^{ + } \) represents the velocity component; \( p^{{{ + }^{\prime } }} \) denotes the pressure fluctuation;\( \beta_{m} \) is the contribution of the mth branching FENE-P model to zero-shear viscosity of viscoelastic solution, \( \beta_{m} { = }{{\eta_{{{\text{Vm}}}} } \mathord{\left/ {\vphantom {{\eta_{{{\text{Vm}}}} } {\eta_{\text{N}} }}} \right. \kern-0pt} {\eta_{\text{N}} }} \), \( \eta_{{{\text{V}}m}} \) and \( \eta_{\text{N}} \) respectively refer to the dynamic viscosity of the solute and solvent of viscoelastic solution; \( Re_{\tau } \) is the frictional Reynolds number, \( Re_{\tau } = {{\rho u_{\tau } h} \mathord{\left/ {\vphantom {{\rho u_{\tau } h} {\eta_{\text{N}} }}} \right. \kern-0pt} {\eta_{\text{N}} }} \); \( We_{\tau ,m} \) is the Weissenberg number, \( We_{\tau ,m} = {{\lambda_{m} \rho u_{\tau }^{2} } \mathord{\left/ {\vphantom {{\lambda_{m} \rho u_{\tau }^{2} } {\eta_{\text{N}} }}} \right. \kern-0pt} {\eta_{\text{N}} }} \), \( \lambda_{m} \) is the relaxation time; \( f\left( {r_{m} } \right) \) denotes nonlinear stretching factor, \( f\left( {\bar{r}_{m} } \right) = {{\left( {L^{2} - 3} \right)} \mathord{\left/ {\vphantom {{\left( {L^{2} - 3} \right)} {\left( {L^{2} - {\text{trace}}\left( {{\mathbf{c}}_{m}^{ + } } \right)} \right)}}} \right. \kern-0pt} {\left( {L^{2} - {\text{trace}}\left( {{\mathbf{c}}_{m}^{ + } } \right)} \right)}} \); \( c_{ij,m}^{ + } \) denotes the component of conformation tensor; the subscript ‘m’ represents the mth branching FENE-P model.

In this study, Eqs. (1)–(3) are filtered by the MICT SGS model. The filtered dimensionless LES governing equation of viscoelastic turbulent drag-reducing channel flow reads

$$ \frac{{\partial \bar{u}_{i}^{ + } }}{{\partial x_{i}^{*} }} = 0 $$
(4)
$$ \begin{aligned} \frac{{\partial \bar{u}_{i}^{ + } }}{{\partial t^{*} }} + \bar{u}_{j}^{ + } \frac{{\partial \bar{u}_{i}^{ + } }}{{\partial x_{j}^{*} }} = & \delta_{1i} - \frac{{\partial \bar{p}^{{ +^{\prime } }} }}{{\partial x_{i}^{*} }} + \frac{1}{{Re_{\tau } }}\frac{\partial }{{\partial x_{j}^{*} }}\left( {\frac{{\partial \bar{u}_{i}^{ + } }}{{\partial x_{j}^{*} }}} \right) + \sum\limits_{m = 1}^{N} {\frac{{\beta_{m} }}{{We_{\tau ,m} }}\frac{{\partial \left[ {f\left( {\bar{r}_{m} } \right)\bar{c}_{ij,m}^{ + } } \right]}}{{\partial x_{j}^{*} }}} \\ & + \,\sum\limits_{m = 1}^{N} {\frac{{\beta_{m} }}{{We_{\tau ,m} }}\frac{{\partial R_{ij,m} }}{{\partial x_{j}^{*} }}} + \frac{{\partial \tau_{ij} }}{{\partial x_{j}^{*} }} \\ \end{aligned} $$
(5)
$$ \begin{aligned} \frac{{\partial \bar{c}_{ij,m}^{ + } }}{{\partial t^{*} }} + \bar{u}_{k}^{ + } \frac{{\partial \bar{c}_{ij,m}^{ + } }}{{\partial x_{k}^{*} }} - \frac{{\partial \bar{u}_{i}^{ + } }}{{\partial x_{k}^{*} }}\bar{c}_{kj,m}^{ + } - \frac{{\partial \bar{u}_{j}^{ + } }}{{\partial x_{k}^{*} }}\bar{c}_{ik,m}^{ + } { = } & \frac{{Re_{\tau } }}{{We_{\tau ,m} }}\left[ {\delta_{ij,m} - f\left( {\bar{r}_{m} } \right)\bar{c}_{ij,m}^{ + } } \right]{ + }P_{ij,m} + Q_{ij,m} \\ & - \frac{{Re_{\tau } }}{{We_{\tau ,m} }}R_{ij,m} + \chi_{\text{c}} \left( {\bar{\gamma }_{ij,m} - \bar{c}_{ij,m}^{ + } } \right) \\ \end{aligned} $$
(6)

where the overbar “—” represents filtering; \( \tau_{ij} \) is the SGS shear stress in ICSM; \( R_{ij,m} \) represents the subfilter term related to nonlinear restoring force; \( P_{ij,m} \) and \( Q_{ij,m} \) are the subfilter terms induced by stretching of microstructures formed in viscoelastic solution; \( \chi_{\text{c}} \left( {\bar{\gamma }_{ij,m} - \bar{c}_{ij,m}^{ + } } \right) \) is a second-order regularization term; \( \chi_{\text{c}} \) denotes dissipative coefficient and is set as 1.0 in this work.

To calculate \( P_{ij,m} \), \( Q_{ij,m} \), \( R_{ij,m} \) and \( \chi_{\text{c}} \left( {\bar{\gamma }_{ij,m} - \bar{c}_{ij,m}^{ + } } \right) \) in Eqs. (5)–(6), the deconvolution velocity \( u_{i}^{*} \) for the approximation of unsolved velocity \( u_{i}^{ + } \), the deconvolution conformation tensor \( \phi_{ij,m} \) and \( \gamma_{ij,m} \) for the unfiltered conformation tensor \( c_{ij,m}^{ + } \) are established respectively as follows

$$ u_{i}^{*} = \sum\limits_{r = 0}^{p} {C_{r} \bar{u}_{i}^{ + (r + 1)} } ,\,\phi_{ij,m} = \sum\limits_{r = 0}^{p} {C_{r} \bar{c}_{ij,m}^{ + (r + 1)} } ,\,\gamma_{ij,m} = \sum\limits_{r = 0}^{q} {D_{r} \bar{c}_{ij,m}^{ + (r + 1)} } $$
(7)

where p and q refer to the deconvolution degrees, taken p = 3 and q = 2 in this study; \( C_{r} \) and \( D_{r} \) represent the optimal deconvolution coefficients corresponding to p and q, they can be calculated according to the binomial theorem, \( \left[ {\begin{array}{*{20}c} {C_{0} ,} & {C_{1} ,} & {C_{2} ,} & {C_{3} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}c} {0,} & {\sqrt 6 ,} & {\sqrt {4 + 2\sqrt 6 } - 2\sqrt 6 ,} & {1 - \sqrt {4 + 2\sqrt 6 } + \sqrt 6 } \\ \end{array} } \right],\left[ {\begin{array}{*{20}c} {D_{0} ,} & {D_{1} ,} & {D_{2} } \\ \end{array} } \right] = \left[ {\begin{array}{*{20}l} {15/8,} \hfill & { - 15/8,} \hfill & {1/4} \hfill \\ \end{array} } \right] \).

Based on Eq. (7), the additional subfilter terms \( P_{ij,m} \) and \( Q_{ij,m} \) can be calculated as fellows

$$ P_{ij,m} = \overline{{\frac{{\partial u_{i}^{*} }}{{\partial x_{k}^{*} }}\phi_{kj,m} }} - \frac{{\partial \bar{u}_{i}^{*} }}{{\partial x_{k}^{*} }}\bar{\phi }_{kj,m} $$
(8)
$$ Q_{ij,m} = \overline{{\frac{{\partial u_{j}^{*} }}{{\partial x_{k}^{*} }}\phi_{ki,m} }} - \frac{{\partial \bar{u}_{j}^{*} }}{{\partial x_{k}^{*} }}\bar{\phi }_{ki,m} $$
(9)

Furthermore, \( R_{ij,m} \) can be computed by the N-parallel FENE-P model below

$$ R_{ij,m} = \overline{{\left( {f\left( {r_{m} } \right)\phi_{ij,m} - \delta_{ij} } \right)}} - \left( {f\left( {\bar{r}_{m} } \right)\bar{\phi }_{ij,m} - \delta_{ij} } \right) $$
(10)

It is worth noting that the filtering of Eqs. (7)–(10) is carried out on the time-domain. Readers can refer to [11] for the details about the temporal filtering.

3 Numerical Approaches

In this study, the finite difference method (FDM) is applied to discretize the dimensionless governing Eqs. (4)–(6). The diffusion terms are discretized using the second-order central difference scheme and the second-order Adams-Bashforth scheme is adopted for time marching. To obtain high resolution numerical solutions physically and maintain the symmetrically positive definiteness of the conformation tensor, a second-order bounded scheme—MINMOD is applied to discretize the convection term in constitutive Eq. (6).

The projection algorithm [13] is utilized to solve the coupled discrete equations. Within this algorithm, the pressure fluctuation Poisson equation, which is constructed on the continuity equation, is directly solved on staggered mesh. The momentum equation is solved in two steps: (1) Ignore the pressure fluctuation gradient and calculate the intermediate velocities; (2) Substitute the intermediate velocities into the pressure fluctuation Poisson equation to obtain the pressure fluctuation. The finial velocities can be obtained by the summation of intermediate velocities and pressure fluctuation gradient. Considering that it is time-consuming to solve the pressure fluctuation Poisson equation with implicit iterations, the geometric multigrid (GMG) method [14] is employed to speed up the computation. The entire calculation procedures for the projection algorithm are presented as follows

  • Step1: set the initial fields of velocities, pressure fluctuation and conformation tensor first, and then calculate the coefficients and constant terms in momentum equation and constitutive equation;

  • Step2: discretize the momentum equation and calculate the intermediate velocities;

  • Step3: substitute the discrete momentum equation into continuity equation to get the discrete pressure fluctuation Poisson equation;

  • Step4: adopt the GMG method to solve the discrete pressure fluctuation Poisson equation to obtain the converged pressure fluctuation on current time layer;

  • Step5: calculate the velocity fields on current time layer with the intermediate velocities and pressure fluctuation;

  • Step6: discretize the constitutive equation and calculate the conformation tensor field on current time layer;

  • Step7: advance the time layer and return to Step 2. Repeat the procedures until reach the prescribed calculation time.

4 Study on DR Mechanism of High Reynolds Turbulent Drag-Reducing Flow with Viscoelastic Fluids

4.1 Calculation Conditions

Limited by the computer hardware and numerical approaches, previous studies on turbulent DR mechanism of viscoelastic fluid mainly focused on low Reynolds numbers. In engineering practice, however, the Reynolds number is always up to the order of 104. To better provide beneficial guidance to the engineering application of turbulent DR technology, the flow characteristics and DR mechanism in high Reynolds viscoelastic turbulent drag-reducing flow are investigated in present work. To achieve an overall perspective, we design seven test cases in which cases V1–V4 have different Weissenberg numbers while cases V2, V5, V6 have different solution concentrations. The detailed calculation parameters of all test cases are presented in Table 1, where N represents the Newtonian fluid (take water as example) and V denotes the viscoelastic fluid (take surfactant solution as example). In the LES simulation, the double parallel FENE-P model is adopted.

Table 1. Calculation parameters of LES of turbulent drag-reducing channel flow.
Full size table

For the boundary conditions of the test cases, the periodic boundary is adopted along both the streamwise and spanwise directions, non-slip boundary is imposed on the channel walls. Furthermore, it is essential to give the initial conditions in TADM SGS model. In this paper, the initial conditions are set as the same in [13]. To balance the computational burden and numerical accuracy, a mesh with 32 × 64 × 32 grid points is adopted in LES. The grid independence test can be found in [9, 11, 12] and the mesh has been proved dense enough to sustain the turbulence. It is worth noting that the uniform grid is adopted in streamwise and spanwise directions and non-uniform grid is used in wall-normal direction with denser mesh near the wall due to the ‘wall effect’. In the LES the maximum stretching length of the microstructures formed in surfactant solution is set as L = 100. The spatial filter width is set as 0.0634 and the temporal filter widths for velocity and conformation tensor are both set as 10 \( \Delta t^{ * } \).

4.2 Results Analysis and Discussion

  1. (1)

    Drag-reduction rate

    Compared with Newtonian fluids, the noticeable characteristic of viscoelastic fluid is DR effect, and it can be described quantitatively with the DR rate. In Table 2, the calculation results including dimensionless mean velocity, Reynolds number, frictional coefficient and DR rate of cases V1–V6 are presented in detail. We can easily observe that when Reynolds number reaches 3 × 104, the turbulent DR effect is very remarkable. For instance, the DR rate can reach up to 51.07% only when \( We_{\tau 1} \), \( \beta_{1} \), \( We_{\tau 2} \), and \( \beta_{2} \) are respectively set as 10, 0.1, 10 and 0.1 (case V1).

    Table 2. Results of turbulent drag reduction
    Full size table

    Further analyzing the influence of calculation parameters on the DR rate, some extraordinary phenomena in high Reynolds turbulent drag-reducing flow can be found in comparison with low Reynolds turbulent flow. In high Reynolds number situations, DR rate grows with the increase of \( We_{\tau } \) under constant \( \beta \), similarly with low Reynolds number situations. However, the influence from \( We_{\tau } \) on DR rate in cases V1–V4 is far less than that in low Reynolds flow. The main reason for this phenomenon lies in that the drag-reducing performance of viscoelastic fluid not only depends on \( We_{\tau } \), but also has close relation with the ratio of elastic effect to viscous effect \( {{We_{\tau } } \mathord{\left/ {\vphantom {{We_{\tau } } {Re_{\tau } }}} \right. \kern-0pt} {Re_{\tau } }} \). In high Reynolds flow, the variation of \( {{We_{\tau } } \mathord{\left/ {\vphantom {{We_{\tau } } {Re_{\tau } }}} \right. \kern-0pt} {Re_{\tau } }} \) with \( We_{\tau } \) is not obvious because the denominator is large, in this way the influence from \( We_{\tau } \) on DR rate is weakened. Meanwhile, it is also indicated in cases V5 and V6 that DR rate is not necessarily to increase with the rise of \( \beta \) under constant \( We_{\tau } \). This is because with the increase of solution concentration, the number of microstructures formed in viscoelastic fluid also increases, the elastic resistance induced by the drag additive increases obviously, which can largely offset the DR effect.

  2. (2)

    Root-mean-square of velocity fluctuations

    To further clarify the influence from viscoelasticity on characteristics of turbulent drag-reducing flow, it is necessary to analyze the turbulent fluctuation intensity which can be quantified by the root-mean-square (RMS) of dimensionless velocity fluctuation. Figure 4 shows the distributions of velocity fluctuation versus y+ in cases V1–V4. From the figure it can be obviously seen the RMS of dimensionless streamwise velocity fluctuation is far higher than those of the wall-normal velocity fluctuation and spanwise velocity fluctuation, which indicates that the RMS of dimensionless streamwise velocity fluctuation is dominant. Compared with the Newtonian fluid, the peaks of streamwise, wall-normal and spanwise velocity fluctuations of viscoelastic fluid all move away from the channel wall. However, in high Reynolds flow, the maximum RMS of dimensionless streamwise velocity fluctuation of viscoelastic fluid is almost the same with that of Newtonian fluid, as shown in Fig. 5(a). This phenomenon was also observed in the experiments, it implies the increase of RMS of streamwise velocity fluctuation is not an intrinsic characteristic of turbulent DR. Figure 5(b)–(c) illustrate the dimensionless wall-normal and spanwise velocity fluctuations are suppressed largely by the adding of additive. This trend is more obvious for the dimensionless wall-normal velocity, whose RMS is only 30% of that of Newtonian fluid. The suppression of dimensionless wall-normal velocity fluctuation directly results in the weakening of momentum transfer between turbulence and channel wall, which is one main reason for the DR.

    Fig. 4.
    figure 4

    Profiles of root-mean-square of dimensionless velocity fluctuation with y+.

    Full size image
    Fig. 5.
    figure 5

    Comparison of coherent structures (vortex tube structures) in turbulence at the dimensionless time 80

    Full size image
  3. (3)

    Contribution of shear stresses on the frictional resistance coefficient

    The balance relationship of shear stresses of viscoelastic fluid is different from that of Newtonian fluid, and different shear stresses have different impacts on turbulent frictional resistance. The contribution of shear stresses on the frictional resistance coefficient of viscoelastic fluid can be derived as follows

    $$ C_{f} = \frac{12}{{Re_{\text{b}} }} + \frac{6}{{\left( {U_{\text{b}}^{ + } } \right)^{2} }}\int_{0}^{1} {\left( { - \bar{u}^{{ +^{\prime } }} \bar{v}^{{ +^{\prime } }} } \right)\left( {1 - y^{ *} } \right)dy^{*} } + \frac{6}{{\left( {U_{\text{b}}^{ + } } \right)^{2} }}\int_{0}^{1} {\sum\limits_{m = 1}^{N} {\frac{{\beta_{m} f\left( {\bar{r}_{m} } \right)\bar{c}_{xy,m}^{ + } }}{{We_{\tau ,m} }}\left( {1 - y^{ *} } \right)} dy^{*} } $$
    (11)

    where on the right hand side of the equation, the terms from left to right denote influences of viscous stress, Reynolds shear stress and elastic stress on the frictional resistance coefficient, which are named as viscous contribution (VC), turbulent contribution (TC) and elastic contribution (EC), respectively. Especially, elastic contribution is equal to zero in Newtonian fluid.

    The proportions of VC, EC and TC to the overall turbulent frictional resistance coefficient for the test cases are presented detailedly in Table 3. We can find the frictional resistance coefficient of Newtonian fluid is much larger than that of viscoelastic fluid under same calculation condition. The VC of viscoelastic fluid is higher than that of Newtonian fluid while the TC of the former is much lower than the later, although TC is dominant in both fluids. Different from the low Reynolds situations, the influence of Weissenberg number is not obvious in high Reynolds flows. This is similar to the influence from Weissenberg number on DR rate in Table 2. However, the EC is more sensitive to \( \beta \) and the larger \( \beta \) is, the larger EC would be.

    Table 3. Contributions of different shear stresses on frictional resistance coefficient.
    Full size table
  4. (4)

    Turbulent coherent structure

    Turbulent coherent structure is a series of motions triggered randomly in the turbulent flow, which is closely related to the turbulent burst events. In this work, the following Q method proposed by Hunt et al. [16] is adopted to extract one of coherent structures (vortex tube structure) in the turbulent channel flows

    $$ Q = \frac{1}{2}\left( {\overline{W}_{ij} \overline{W}_{ij} - \bar{S}_{ij} \bar{S}_{ij} } \right) > 0 $$
    (12)

    where \( \overline{W}_{ij} \) is the filtered vorticity tensor, \( \overline{W}_{ij} { = }\frac{1}{2}\left( {\frac{{\partial \bar{u}_{j} }}{{\partial x_{i} }} - \frac{{\partial \bar{u}_{i} }}{{\partial x_{j} }}} \right) \); \( \overline{S}_{ij} \) is the filtered velocity-strain tensor, \( \bar{S}_{ij} { = }\frac{1}{2}\left( {\frac{{\partial \bar{u}_{j} }}{{\partial x_{i} }} + \frac{{\partial \bar{u}_{i} }}{{\partial x_{j} }}} \right) \).

Figure 5 demonstrates the comparison of vortex tube structures at the dimensionless time 80 (the turbulent flow is fully developed) with different Q values in cases N1 and V2. From the figure, it can be seen that the number of vortex tube structures in viscoelastic fluid is less than that in Newtonian fluid. Meanwhile, more vortex tube structures in viscoelastic fluid are in large scale. This indicates that the adding of additive in high Reynolds turbulent channel flow can largely suppress the formation of turbulent coherent structures, which further suppresses the intermittency of turbulent flow. In this way, the frequency and intensity of turbulent burst events are largely weakened, this phenomenon is also similar to the low Reynolds situations.

To investigate the intermittency of turbulent flow quantitatively, the skewness factor and flatness factor of streamwise velocity fluctuation are calculated below

$$ S\left( u \right) = {{\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{3} } \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{3} } \right\rangle } {\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{2} } \right\rangle^{1.5} }}} \right. \kern-0pt} {\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{2} } \right\rangle^{1.5} }} $$
(13)
$$ F\left( u \right) = {{\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{4} } \right\rangle } \mathord{\left/ {\vphantom {{\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{4} } \right\rangle } {\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{2} } \right\rangle^{2} }}} \right. \kern-0pt} {\left\langle {\left( {\bar{u}^{{ +^{\prime } }} } \right)^{2} } \right\rangle^{2} }} $$
(14)

Table 4 gives the average skewness factor and flatness factor of the seven test cases. The result of skewness factor indicates the turbulent flow field distributions for Newtonian fluid and viscoelastic fluid are both asymmetrical. However, the absolute value of skewness factor for viscoelastic fluid is larger than that of the Newtonian fluid, which indicates a higher asymmetry and larger deviation from Gaussian field in the viscoelastic fluid. It can also be found the flatness factor of Newtonian fluid is larger than that of viscoelastic fluid, indicating that the probability density function of streamwise velocity in viscoelastic fluid is flatter than that of Newtonian fluid. Therefore, the overall intermittency of turbulent channel flow is suppressed and the turbulent fluctuation and burst events are obviously weakened.

Table 4. Calculation results of average skewness factor and flatness factor.
Full size table

5 Conclusions

In this paper, we apply the LES approach to investigate turbulent DR mechanism in viscoelastic turbulent drag-reducing flow with high Reynolds number based on an N-parallel FENE-P model and MICT SGS model. The main conclusions of this work are as follows

  1. (1)

    There is a notable difference of flow characteristics between high Reynolds and low Reynolds turbulent drag-reducing flows. For instance, the RMS peak value of dimensionless streamwise velocity fluctuation in high Reynolds turbulent flow is not obviously strengthened as in low Reynolds turbulent flow.

  2. (2)

    The influence of calculation parameters on frictional resistance coefficients in high Reynolds turbulent drag-reducing flow also differs from that in low Reynolds situations. For example, the effect of Weissenberg number exerted to turbulent flow is not obvious while influence of \( \beta \) on elastic contribution is more remarkable.

  3. (3)

    The number of vortex tube structures in high Reynolds turbulent drag-reducing flow with viscoelastic fluid is much less than that in Newtonian turbulent flow. Compared with the Newtonian fluid, the turbulent flow field of viscoelastic fluid is featured by higher asymmetry and more large-scale motions with low frequency.