Reynolds Stress Model for Viscoelastic Drag-Reducing Flow Induced by Polymer Solution

Viscoelasticity drag-reducing flow by polymer solution can reduce pumping energy of pipe flow significantly. One of the simulation manners is direct numerical simulation (DNS). However, the computational time is too long to accept in engineering. Turbulent model is a powerful tool to solve engineering problems because of its fast computational ability. However, its precision is usually low. To solve this problem, we introduce DNS to provide accurate data to construct a high-precision turbulent model. A Reynolds stress model for viscoelastic polymer drag-reducing flow is established. The rheological behavior of the drag-reducing flow is described by the Giesekus constitutive Equation. Compared with the DNS data, mean velocity, mean conformation tensor, drag reduction, and stresses are predicted accurately in low Reynolds numbers and Weissenberg numbers but worsen as the two numbers increase. The computational time of the Reynolds stress model (RSM) is only 1/120,960 of DNS, showing the advantage of computational speed.


Introduction
Drag reduction (DR) phenomenon was first discovered by Toms [1]. He observed in his experiment that the addition of a long-chain polymer (polymethyl methacrylate) in monochlorobenzene dramatically reduced the turbulent skin friction by as high as 80%. The flow rate could be increased by the addition of the polymer at constant pressure gradient. Then, he reported these results at the First International Rheological Congress, so it is usually referred to as the "Toms Effect". The polymers that can reduce skin friction were later called drag-reducing agents (DRAs). The energy-saving effect of DRAs attracts many applications. The first famous application for polymer drag reduction was its use in the 48 inch diameter 800 mile length Alaska pipeline, carrying crude oil from the North slope in Alaska to Valdez in the south of Alaska [2]. After injecting a concentrated solution of a high-molecular-weight polymer downstream of pumping stations at homogeneous concentrations as low as 1 ppm [3], crude throughput was increased by up to 30%. Polymer DRAs were also successfully applied in other crude oil pipelines such as Iraq-Turkey, Bass Strait in Australia, Mumbai Offshore [4], and North Sea Offshore [5], and in finished hydrocarbon product lines [6].
The drag-reducing flow induced by polymer solution usually appears viscoelastic. Direct numerical simulation (DNS) can simulate this kind of viscoelastic turbulent flow in high precision [7][8][9][10][11]. More recent progresses are as follow. Dubief et al. [12] investigated the energetics of turbulence by correlating the work done by polymers on the flow with turbulent structures. Polymers are found to store and to release energy to the flow in a well-organized manner. Graham [13] proposed a tentative unified description of rheological drag reduction based on his numerical observations. Thais et al. [14] found that the spectra of its cross-flow component in viscoelastic flows exhibit a significantly higher energy level at a large scale. Pereira et al. [15] studied the polymer-turbulence interactions from an energetic standpoint for a range of Weissenberg numbers and found a cyclic mechanism of energy exchange between the polymers and turbulence that drives the flow through an oscillatory behavior. They also addressed the numerical simulation of thermo-fluid characteristics of triangular jets [16]. However, DNS needs numerous storages of the computer because very dense mesh is required to resolve small eddies in turbulent drag-reducing flow, such that computational time is too long to be accepted in engineering. Turbulent model computes the turbulent flow very quickly. Compared with DNS, turbulent model for viscoelastic drag-reducing flows develops slowly. There were zero-Equation models established by Edwards et al. [17] and Azouz et al. [18]. One-Equation models and two-Equation models were established by Durst et al. [19], and Hassid and Poreh [20][21][22]. Cruz et al. [23] and Pinho et al. [24] considered elongation thickening of drag-reducing fluid based on Newtonian fluid turbulent flow and derived a new low Reynolds number k-ε model for polymer drag-reducing flow. Elongation thickening is a very important factor of drag reduction. Thus, Pinho's work promoted the turbulent model greatly, but the precision is still not high enough due to the complexity of viscoelastic turbulent flow. Reynolds stress model (RSM) can simulate Newtonian turbulent flow in high precision, so that it has potential advantages to deal with the complex viscoelastic turbulent flow with polymer additives. If the RSM is established based on DNS data, the precision may be better.
In this paper, an RSM for viscoelastic drag-reducing flow is established based on DNS data. The goal is to find a new modeling way to solve drag-reducing flows in engineering with fast computation and good accuracy.

Instantaneous Equations
Viscoelastic drag-reducing flow can be described by the following governing Equations [25].
(1) Continuity Equation: ∂u (2) Momentum Equation: (3) Giesekus constitutive Equation: where ρ and λ are the density and the relaxation time of drag-reducing fluid respectively. α is the mobility factor, which determines the extensional viscosity. p is pressure. u i (i = x, y, z) are the velocity components in the x, y, z directions. c ij is the conformation tensor. µ is the zero-shear-rate viscosity of solvent. η is the zero-shear-rate viscosity of drag-reducing polymer. We choose the following dimensionless transformation: and use the definition of Weissenberg number We τ = ρλu 2 τ /(µ + η), the definition of frictional Reynolds number Re τ = ρu τ h/(µ + η), and the ratio β = µ/(µ + η). u τ is the frictional velocity (u τ = τ w /ρ, τ w is the wall shear stress). h is the half height of the channel, as shown in Figure 1. Using these definitions, Equations (1)-(3) can be transformed to be the following dimensionless Equations: The above Reynolds stress model was used to simulate the fully developed viscoelastic dragreducing channel flow. The computational domain is shown in Figure 1. Periodic boundary conditions were imposed in both the streamwise (x-) and spanwise (z-) directions, while nonslip boundary conditions were adopted for the top and bottom walls. Computational parameters were:  The numerical method of DNS is a fractional step method. Adams-Bashforth scheme was used to ensure the second-order accuracy of velocity. An implicit scheme was used for the pressure term. Staggered grid was applied to avoid unphysical oscillations of pressure. A second-order finite difference scheme was used for spatial discretization. Uniform mesh was used in the x and z directions due to the periodic boundary condition. To capture small eddies near the walls, nonuniform mesh was used in the y direction. Grid number was 64 × 64 × 64. Drag reduction is defined as: where Cf is the calculated friction factor, CfDean is evaluated by Dean's correlation [28], DR% is the drag reduction.
Mean streamwise velocity, drag reduction, Reynolds stress, and fluctuation intensity were obtained and compared with the results of DNS to validate the model. Bulk mean variables are listed in Table 1. The results obtained by the turbulent model (RSM) agree well with the DNS results. The relative deviations of mean streamwise velocity, Reynolds number, frictional factor, and drag reduction were 0.57%, 0.57%, 1.27%, and 9.3%. Thus, the prediction of the bulk mean variables by RSM is accurate.   (4) Dimensionless continuity Equation: (6) Dimensionless Giesekus constitutive Equation:

Time-Average Equations
Reynolds stress model focuses on the time-average effect of turbulent flow. Thus, instantaneous variables in the above Equations can be considered as the summation of time-average variables and fluctuation variables, that is, ϕ + = ϕ + + ϕ + (ϕ + represents instantaneous variables (u + i , p + , c + ij etc.); superscripts "-" and "'" represent time-average variables and fluctuation variables respectively). After this decomposition, time-average operations can be made for the instantaneous Equations to obtain the following time-average Equations: (1) Time-average momentum Equations: (2) Time-average constitutive Equations: Equations (2)-(7) so that we can obtain the following fluctuation Equations. (3) Fluctuation Equations: Equation (9) can also be rewrote as: Equation (9) × u + j + Equation (10) × u + i and do the time average, we can obtain the following Reynolds stress transport Equations.
(4) Reynolds stress transport Equations: All these Equations contain high-order moments A ij , B ij , C ij , φ ij , T ij , ε ij , E ij , F ij and u + i u + j . They are new unknowns. Apparently, the number of unknowns exceeds the number of Equations, so that all the high-order unknowns need to be modeled as functions of time-average variables.

Modeling of High-Order Moments of Fluctuations
Different from Newtonian turbulent flow, the above high-order moments of drag-reducing flow are all related to viscoelasticity. The terms directly related to viscoelasticity (A ij , B ij , C ij , E ij , F ij ) do not appear in the turbulent models of Newtonian flow, such that they are modeled for the first time. The other terms (φ ij , T ij , ε ij ) have implicit and complex relations with viscoelasticity, such that the modeling manner of these terms uses the same way of Newtonian flow. The modeling process is as follows.

High-Order Moments Directly Related to Viscoelasticity
The additional turbulent diffusion terms introduced by viscoelasticity (A ij , E ij ) and the correlation of conformation tensor B ij are much smaller than turbulent diffusion. Thus, they can be neglected. To conform the physical process and consider the easy use of the model, C ij is assumed to have relations with Reynolds stress, mean strain rate, and so on. According to the positive definiteness of Reynolds stress and nonpositive definiteness of mean strain rate, the dimensional form of C ij can be expressed as: ) and Kolmogorov velocity scale ((εν) 1/4 ). ν (= µ/ρ) is kinetic viscosity of fluid.
Equation (8) is nondimensionalized using ε + = ε/ u 3 τ /h , fluid density, viscosity, and frictional velocity and combined with the definitions of Re τ , We τ , and β. The model of C ij is as follows: where ϕ 1 = 0.05, ϕ 2 = −0.01. From observations, F ij and C ij have the following simple relation:

Results and Discussion
The above Reynolds stress model was used to simulate the fully developed viscoelastic drag-reducing channel flow. The computational domain is shown in Figure 1. Periodic boundary conditions were imposed in both the streamwise (x-) and spanwise (z-) directions, while nonslip boundary conditions were adopted for the top and bottom walls. Computational parameters were: Re τ = 150, We τ = 10, α = 0.001, β = 0.8.
The numerical method of DNS is a fractional step method. Adams-Bashforth scheme was used to ensure the second-order accuracy of velocity. An implicit scheme was used for the pressure term. Staggered grid was applied to avoid unphysical oscillations of pressure. A second-order finite difference scheme was used for spatial discretization. Uniform mesh was used in the x and z directions due to the periodic boundary condition. To capture small eddies near the walls, nonuniform mesh was used in the y direction. Grid number was 64 × 64 × 64. Drag reduction is defined as: where C f is the calculated friction factor, C fDean is evaluated by Dean's correlation [28], DR% is the drag reduction. Mean streamwise velocity, drag reduction, Reynolds stress, and fluctuation intensity were obtained and compared with the results of DNS to validate the model. Bulk mean variables are listed in Table 1. The results obtained by the turbulent model (RSM) agree well with the DNS results. The relative deviations of mean streamwise velocity, Reynolds number, frictional factor, and drag reduction were 0.57%, 0.57%, 1.27%, and 9.3%. Thus, the prediction of the bulk mean variables by RSM is accurate.  Figure 2 is the comparison of time-average streamwise velocity profiles between RSM and DNS. The DNS results for Newtonian turbulent flow agree well with the typical distribution of viscous sublayer (u + = y + ), buffer layer (u + = 5 ln y + − 3.05), and turbulent core region (u + = 2.5 ln y + + 5.5) [29], showing that the fluid has achieved fully developed turbulent flow. The time-average velocity of RSM coincides well with DNS in the viscous sublayer, is lower than DNS in the buffer layer, and is higher than DNS in the turbulent core region. The mean deviation is only 9.1%, which is much lower than previous turbulent models (as high as 30-50% or more). This indicates the established model can predict the time-average streamwise velocity profile accurately. Fluctuation intensities are compared in Figure 3. Predicted v + rms and w + rms are larger than DNS, while u + rms is smaller than DNS. The peak value positions of u + rms v + rms , and w + rms of RSM are basically the same as DNS.
the buffer layer, and is higher than DNS in the turbulent core region. The mean deviation is only 9.1%, which is much lower than previous turbulent models (as high as 30-50% or more). This indicates the established model can predict the time-average streamwise velocity profile accurately. Fluctuation intensities are compared in Figure 3. Predicted rms v + and rms w + are larger than DNS, while rms u + is smaller than DNS. The peak value positions of rms u + rms v + , and rms w + of RSM are basically the same as DNS. Stress balance of Reynolds stress, viscous stress, and viscoelastic stress is shown in Figure 4. Prediction results of the three stresses agree well with those of DNS. Predicted Reynolds stress is smaller than DNS. Predicted viscoelastic stress is almost the same as that of DNS. which is much lower than previous turbulent models (as high as 30-50% or more). This indicates the established model can predict the time-average streamwise velocity profile accurately. Fluctuation intensities are compared in Figure 3. Predicted rms v + and rms w + are larger than DNS, while rms u + is smaller than DNS. The peak value positions of rms u + rms v + , and rms w + of RSM are basically the same as DNS. Stress balance of Reynolds stress, viscous stress, and viscoelastic stress is shown in Figure 4. Prediction results of the three stresses agree well with those of DNS. Predicted Reynolds stress is smaller than DNS. Predicted viscoelastic stress is almost the same as that of DNS. Stress balance of Reynolds stress, viscous stress, and viscoelastic stress is shown in Figure 4. Prediction results of the three stresses agree well with those of DNS. Predicted Reynolds stress is smaller than DNS. Predicted viscoelastic stress is almost the same as that of DNS.  C mij represents the mean conformation tensor c + ij . The four main components of the conformation tensor, with subscripts 1, 2, 3 representing streamwise, wall-normal, and spanwise directions, are compared in Figure 5. It shows that the four components agree well with DNS. The peak value of C m11 is 15% lower than the value in DNS. C m12 of RSM coincides well with DNS. It determines the value of viscoelastic stress, such that the deviation of viscoelastic stress is small in Figure 4. The RSM for the above case is basically verified. We further verify the RSM using the DNS results in literature [30] at different Weissenberg numbers. Table 2 shows the bulk mean variables such as drag reduction. It is apparent that the RSM can also simulate the drag-reducing flow accurately at Weτ = 12.5 and 30, which correspond to low and medium drag reduction, respectively. Figure 6 and Figure 7 show that the critical features (e.g., average velocity profile and Reynolds stress distribution) can be simulated by the RSM in satisfied precision.  The RSM for the above case is basically verified. We further verify the RSM using the DNS results in literature [30] at different Weissenberg numbers. Table 2 shows the bulk mean variables such as drag reduction. It is apparent that the RSM can also simulate the drag-reducing flow accurately at We τ = 12.5 and 30, which correspond to low and medium drag reduction, respectively. Figures 6 and 7 show that the critical features (e.g., average velocity profile and Reynolds stress distribution) can be simulated by the RSM in satisfied precision.  RSM only concerns time-average variables so that it does not need to resolve the small eddies near the walls like DNS. The complex and time-consuming numerical methods, such as numerous iterations of pressure, can be avoided. Therefore, the computational time of RSM can be much smaller than that of DNS. Table 3 verifies that RSM can largely save the computational time because the acceleration ratio of computational time of DNS and RSM is as high as 120,960.

DNS
RSM Acceleration Ratio 604,800s 5s 120,960 To further examine the precision of the RSM at higher Reynolds numbers, it is compared with the typical DNS results made by Housiadas et al. [31] in Table 4. The drag reduction of the RSM agrees well with the DNS data in low (Reτ = 125) and medium (Reτ = 180) Reynolds numbers. The deviation becomes larger at higher Reynolds number (Reτ = 395). The model was also used to predict the effect of Weissenberg number at the higher Reynolds number (Table 5). Drag reduction increases with increasing Weτ and tends to a stable value about more than 40%. RSM only concerns time-average variables so that it does not need to resolve the small eddies near the walls like DNS. The complex and time-consuming numerical methods, such as numerous iterations of pressure, can be avoided. Therefore, the computational time of RSM can be much smaller than that of DNS. Table 3 verifies that RSM can largely save the computational time because the acceleration ratio of computational time of DNS and RSM is as high as 120,960. To further examine the precision of the RSM at higher Reynolds numbers, it is compared with the typical DNS results made by Housiadas et al. [31] in Table 4. The drag reduction of the RSM agrees well with the DNS data in low (Re τ = 125) and medium (Re τ = 180) Reynolds numbers. The deviation becomes larger at higher Reynolds number (Re τ = 395). The model was also used to predict the effect of Weissenberg number at the higher Reynolds number (Table 5). Drag reduction increases with increasing We τ and tends to a stable value about more than 40%. Table 4. Drag reduction (DR%) comparison with Housiadas' results [31] at different Re τ .

Conclusions
Through the discussions, it is apparent that the newly established Reynolds stress turbulent model can predict bulk mean velocity and drag reduction accurately. The predictions of stresses and mean conformation tensor are also good. The time-average streamwise velocity outside the viscous sublayer is quite different from that in DNS. Streamwise fluctuation intensity is much smaller than that in DNS. This is probably because the modeling of the redistribution term φ ij is the same as that for Newtonian fluid. This term may affect viscoelasticity and pressure in a complex way, which needs further study in future to find a more accurate modeling manner. The acceleration ratio is as high as 120,960, showing the extremely high speed of the new turbulent model of polymeric drag-reducing flow. The Reynolds stress model appears to have good accuracy at lower Reynolds numbers compared with typical DNS results. The model predictions significantly worsen as the Reynolds number and Weissenberg number increase, limiting the applicability of the model up to modest levels of drag reduction. This suggests (in association with the other limitations already stated before) that there is room for model improvement in future work.