Proposal of Damping Function for Low-Reynolds-Number 𝑘-𝜀 Model Applicable in Prediction of Turbulent Viscoelastic-Fluid Flow

Abstract

A low-Reynolds-number k-𝜀 model applicable for viscoelastic fluid was proposed to predict the frictional-drag reduction and the turbulence modification in a wall-bounded turbulent flow. In this model, an additional damping function was introduced into the model of eddy viscosity, while the treatment of the turbulent kinetic energy (k) and its dissipation rate (𝜀) is an extension of the model for Newtonian fluids. For constructing the damping function, we considered the influence of viscoelasticity on the turbulent eddy motion and its dissipative scale and investigated the frequency response for the constitutive equation based on the Giesekus fluid model. Assessment of the proposed model’s performance in several rheological conditions for drag-reduced turbulent channel flows demonstrated good agreement with DNS (direct numerical simulation) data.

1. Introduction

It is known that minute amount of additives, such as polymer and surfactant, can reduce the frictional drag in wall-bounded turbulent flows at high Reynolds numbers. This effect, often referred to as “Toms effect” or “turbulent drag reduction (DR),” has been focused as an efficient energy-saving technology for a large variety of applications (e.g., oil pipelines and district heating/cooling systems). The possibility of substantially decreasing the energy consumption for the transport of liquids stimulated intensive studies of this effect. Many investigators have considered the Toms effect to be a phenomenon that is closely related to a viscoelasticity of the dilute polymer (or surfactant) solution and to a modification of the turbulent structures in the flow. Predictions of the level of DR as well as the modulated turbulent flow of viscoelastic liquids are nontrivial issues to apply the Toms effect in practical situations.

Since the direct numerical simulation (DNS) is one of the important tools to investigate turbulence phenomena qualitatively and quantitatively, DNSs of viscoelastic flows have been increasingly performed by a number of researchers after the pioneering works in 1997-1998 [1, 2]: see reviews [3, 4]. They used several kinds of viscoelastic models for a rheological constitutive equation, such as the FENE-P (finitely extensible nonlinear elastic-Peterlin) model, the Oldroyd-B model, and the Giesekus model, to describe the behaviour of viscoelastic fluids. Most DNS studies aimed at understanding the mechanism of DR and actually revealed some important characteristics common to turbulent channel and pipe flows subject to DR. For instance, Kim and Sureshkumar [5] recently proposed a plausible mechanism associated with the effect of dynamic interactions between turbulent vortical structures and viscoelastic stress on DR. One of their interesting conclusions is the fact that the interplay between the turbulent-eddy (turnover) time scale and the relaxation time of viscoelastic fluid produces the differences in the flow characteristics between low DR regime and high DR regime, which are, respectively, observed at low and high Weissenberg numbers. In viscoelastic flows with a relaxation time that is small but larger than the vortex time scale in the vicinity of a wall, near-wall vortices can be affected by the viscoelastic stress, whereas weaker and larger vortices in the outer layer remain unaffected. With increasing relaxation time, the eddies with longer time scales would interact with the viscoelasticity and the outer-layer modification would be more pronounced. A reason of the existence of maximum DR limit was also explained by Kim and Sureshkumar: the viscoelastic effects eventually encompass all dynamically relevant vortex time scales with very high Weissenberg numbers and, hence, the DR saturates. This result exemplifies the Lumley’s time criteria [6].

As for high Reynolds number simulations, the available information on the viscoelastic turbulent flow is even more scarce than for the Newtonian turbulence, despite of the development of computer resources in the recent decades. A noticeable exception is the DNS performed by Thais et al. [7, 8], where high DR flows at the friction Reynolds number Re_𝜏=1000 was achieved and the Reynolds number similarity was investigated systematically. Viscoelastic flows through complex geometries are also of practically importance and interesting subjects. In flows around bluff bodies, the fluid viscoelasticity would significantly distort the streamlines of the mean flow, regardless of either laminar or turbulent state. Tsukahara et al. [9, 10] demonstrated the flow modification as well as the turbulence modulation in an orifice flow of the viscoelastic fluid, by means of DNS. However, both the Reynolds number and the Weissenberg number examined by the existing DNS studies are still lower than typical level in industrial applications. Therefore, the prediction based on the physics (or the governing equation) of the viscoelastic turbulent flow should be improved for applying the practical flows. Although the use of DNS can be extremely computationally expensive for most engineering applications, DNS databases are useful as calibration/comparison database when considering the application of Reynolds-averaged Navier-Stokes (RANS) simulations. The development of mathematical and computational models for viscoelastic flow is not sufficiently well advanced to allow RANS simulation to be undertaken. With this background, it is then a matter of interest to build more physically-based turbulence closures for the prediction of flows with drag-reducing additives with the aid of the DNS database and the already-proposed mechanism of DR.

Several researchers have attempted to develop more general eddy viscosity or Reynolds stress transport closures for viscoelastic fluids. Leighton et al. [11] developed a Reynolds stress model for the drag-reducing viscoelastic fluids described by the FENE-P rheological constitutive equation, as motivated by the need for robust predictive tools. Their modeling approach focused on two kinds of expected polymer effects on the turbulent flow: one is what they call the “explicit” effect that would be expressed as a new term containing the interactions between fluctuations of the polymer stress and kinematic quantities; the other is the “implicit” effect pertaining to the slow pressure-strain (redistribution) term in the Reynolds-stress transport equation. Pinho and coworkers [12–14] have developed a 𝑘-𝜀 model for the FENE-P fluid, using the Reynolds averaged form of the conformation tensor equation to determine the average polymer stress. They devoted their efforts to devise the closures for the new terms appearing in the governing equations, especially, the nonlinear turbulence distortion term (represented by Λ_𝑖𝑗 in this paper) in the averaged equation of the conformation tensor. Their closure models were calibrated well on the basis of DNS database, Iaccarino et al. [15] employed the 𝑣²-𝑓 approach to treat the intensed wall damping of wall-normal turbulence in the drag-reduced flow. Moreover, they modelled directly the Reynolds-averaged polymer stress without solving the constitutive equation for the conformation tensor, resulting in less coefficients and functions. The results predicted by these models were generally in good agreement with DNS database for low and intermediate DR. More recently, Resende et al. [16] proposed a 𝑘-𝜔 model based on the model of Bredberg et al. [17] for Newtonian fluids. Some improvement was achieved for predictions in the intermediate DR region, but still not enough to estimate 𝑘 and the high DR (in which the drag reduction rate is more than 60%).

The present objective is to construct RANS that can be applied for the viscoelastic turbulent flow by adding corrections to a low Reynolds number 𝑘-𝜀 model of Newtonian-fluid flow (here, 𝑘 and 𝜀 are the turbulent kinetic energy and its dissipation rate, resp.). Two major corrections to be applied to the low Reynolds number 𝑘-𝜀 model, which has been well designed and employed widely in commercial softwares, are the addition of a damping function based on the Lumley’s time criteria and the simultaneous calculation of the constitutive equation of the Giesekus model. The first correction is proposed from the fact that one of the most important effects of the fluid viscoelasticity on turbulence is suppression of small-scale vortical motions near the wall. Hence, the knowledge of scales of eddies that are suppressed by additives is indispensable to discuss the form of the present model. We have obtained this information from existing DNS database [18] enabling a priori development of turbulence closures. To analyze the data, we performed an investigation of frequency response of the Giesekus model and considered the relationship between the shift of the dissipation scale and the viscoelastic behavior of the fluid.

The following section provides a brief introduction of the currently adopted model, which is originally for the Newtonian flow. Thereafter, the scale shifts of the dissipation range and of the viscoelastic behavior are considered and the present constitutive equation as well as the 𝑘 and 𝜀 equations are then developed. Finally, the proposed model will be tested in comparison with DNS results.

2. Introduction of Low-Re 𝑘-𝜀 Model

The present section commences with a brief review of established models for turbulent flows of Newtonian fluid using two-equation turbulent viscosity closures. In the following, the vector 𝑥_𝑖=(𝑥₁,𝑥₂,𝑥₃)≡(𝑥,𝑦,𝑧) denotes the coordinate directions; 𝑢_𝑖=(𝑢₁,𝑢₂,𝑢₃)≡(𝑢,𝑣,𝑤) is the velocity fluctuating component, while those capital letters represent the mean components. The overbar indicates an ensemble average over all homogeneous directions and time.

The low Reynolds number 𝑘-𝜀 model for a turbulent flow relates the Reynolds stress to the mean strain field as follows: −𝑢_𝑖𝑢_𝑗=𝜈_𝑡𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝑈_𝑗𝜕𝑥_𝑖−23𝑘𝛿_𝑖𝑗, (1) where 𝛿_𝑖𝑗 is the Kronecker delta and the turbulent kinematic viscosity is expressed as 𝜈_𝑡=𝐶_𝜇𝑓_𝜇𝑘²𝜀̃. (2) Here, 𝐶_𝜇 is a constant with value equal to 0.09. A number of proposals for the damping function 𝑓_𝜇 have been made; for instance, Jones and Launder [19] functionalized 𝑓_𝜇 using a turbulent Reynolds number Re_𝑡 as follows: 𝑓_𝜇=exp−2.51+Re_𝑡./50 (3) On the other hand, some researchers introduced the wall-normal distance normalized in wall unit, 𝑦⁺=𝑦𝜌𝑢_𝜏/𝜇₀, as a parameter of 𝑓_𝜇. A typical example is the one proposed by Nagano and Tagawa [20] as follows: 𝑓_𝜇=−𝑦1−exp⁺261+4.1Re_𝑡^3/4. (4) Their damping function is based on the van Driest wall-damping formula [21]. They considered the behavior of turbulence on an infinite flat plate with a simple harmonic oscillation parallel to the plate. The amplitude of the motion of the fluid in each position was examined to describe the damping function. Figure 1 shows the 𝑓_𝜇 distribution of these individual models when applied to the turbulent channel flow at a given Reynolds number. The model proposed by Nagano and Tagawa [20] is in better agreement with the DNS data for the Newtonian flow [22], especially, in the near-wall region.

Figure 1 

Comparison of damping function 𝑓𝜇  for turbulent channel flow at Re𝜏=395.

As for a viscoelastic-fluid flow accompanied by DR, both timescale and lengthscale of turbulent-eddy motions are known to be different from those in Newtonian flows. Therefore, a modification of 𝑓_𝜇 for viscoelastic flow is necessarily required. Here, we have investigated the predictive values of 𝑓_𝜇 for viscoelastic fluids using our DNS database [18], which provides various statistical data of the drag-reducing turbulent channel flow for several different rheological properties. Figure 1 shows that 𝑓_𝜇 is dramatically suppressed in any drag-reducing condition. Such aspect is attributed to diminishing turbulent-eddy motions due to the elasticity of fluid. Cruz and Pinho [12] proposed a damping function for a pipe flow of the viscoelastic fluid. They deduced the damping function in the similar way with van Driest [21] and devised an equation including two damping functions to take into account both the fluid rheology and the wall effect into the eddy viscosity. Their function is written as 𝑓_𝜇=𝑓_𝜇𝑣⋅𝑓_𝜇𝑒=1−1+1−𝑛𝑦1+𝑛^{+−(1/26.5)|(1+𝑛)/(1−𝑛)|}⋅1−1+𝑞−1𝑦3−𝑞⁺𝐶_𝑛^{(1−𝑞)/(2−𝑞)−(1/26.5)|(3−𝑞)/(𝑞−1)|}, (5) where 𝑛 is the viscosity power index of the shear behavior, 𝑞 is the viscosity power index of the Trouton-ratio behavior, and 𝐶_𝑛 is a parameter dependent on the additive. In the present paper, we also propose a new additional damping function based on the Lumley’s time criteria [6] for the Giesekus viscoelastic-fluid model.

3. Shift of Dissipative Scale

3.1. Lumley’s Consideration of Time Criteria

As many existing studies show, the turbulent eddy motions and ordered structures in drag-reducing viscoelastic flows are known to exhibit time and spatial scales that are different from those of Newtonian flows [3, 4]. As inspired by the work of Kim and Sureshkumar [5], it is clear that a knowledge of scales (especially, dissipation scale of turbulent kinetic energy) of turbulent eddy motions should be useful and indispensable to construct a turbulence model for predicting the drag-reduced flow. The theoretical concepts [6, 23] about the eddy-scale range that can be shifted in the drag-reduced turbulent flow have not yet been utilized fully and, as far as the authors know, none has been tested against DNS quantitatively. This provides the motivation for re-examining this issue. According to Lumley’s consideration in terms of turbulent energy spectrum [6], we discuss here the energy spectra and the scale shift for the turbulent channel flow.

It is assumed that the mean pressure gradient is constant and the turbulent flow field is fully developed, where the channel width is 2𝛿. The one-dimensional energy spectrum of the streamwise velocity fluctuation is defined as ^∞₀𝐸_𝑢𝑢(𝜅)d𝜅=𝑢𝑢, (6) where 𝜅 is a wave number. In the logarithmic layer, there is no relevant length other than the wall-normal height 𝑦 according to the mixing length theory [24] and, hence, the energy-containing large-scale eddies should be scaled roughly with 𝑦. With respect to the energy spectrum, its peak is expected to occur at a wave number 𝜅_𝐿 given by 𝜅_𝐿⋅𝑦≃1, (7) which then represents the limit of the production range of the spectrum. On the other hand, the scale of dissipative small-scale eddies can be estimated from the energy dissipation spectrum, described as 𝜅²_𝑖𝐸_𝑢𝑢(𝜅) [25]. Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonian-fluid turbulence would occur at 𝜅_𝜂𝜂≃0.2, (8) where 𝜂=(𝜈³/𝜀)^1/4 is the Kolmogorov scale. Here, if the turbulence production and the viscous dissipation rate 𝜀 are locally balanced for each flow unit, we have 𝑢𝜀≃^3𝜏𝜅_𝐾,⋅𝑦 (9) using the von Kármán constant 𝜅_𝐾 and the friction velocity 𝑢_𝜏 [24]. Therefore, the relation between 𝑦 and 𝜅_𝜂 is rewritten from (8) as 𝜅⁺_𝜂𝜅_𝐾⋅𝑦^+1/4≃0.2. (10) The size of turbulent fluctuations (or eddies) ranges between 𝜅_𝐿 and 𝜅_𝜂. Figure 2 shows both wave numbers of 𝜅_𝐿 and 𝜅_𝜂 as a function of 𝑦⁺ in the ordinate. If 𝜅_𝐾=0.4, the lines of (7) and (10) intersect at 𝑦⁺=6.3. Here, we call this point as A. Above the point A, the momentum transport by turbulent vortical motions is active. As 𝑦⁺ being very close to A, the energy productive scale and the dissipative one become comparable with each other, resulting in the unsustainable turbulence at the relevant heights. It can be believed that the upper bound of the viscous sublayer corresponds to the point A, below which the viscous stress is dominant in the total stress.

Figure 2 

Scaling relation in a wall-bounded turbulence: Lumley’s concept [6].

In the drag-reducing viscoelastic flow, the relationship of 𝜅_𝐿 and 𝜅_𝜂 is expected to be changed as shown with broken lines in Figure 2, according to Lumley [6]. He reasoned that less difference between 𝜅⁺_𝐿 and 𝜅⁺_𝜂 as a result of decreased 𝜅⁺_𝜂 might lead to the expanding of the viscous sublayer. The intersection of 𝜅⁺_𝐿 and 𝜅⁺_𝜂, therefore, shifts from A to a higher position labeled as B. Assuming that B is the shifted bound of the viscous sublayer in the relevant flow, the viscous sublayer of the viscoelastic fluid should be enlarged. In other words, the new region, so-called the elastic layer [26], arises in the vicinity of the wall with the upper bound at the height of B.

In order to verify this concept of DR, we reexamine the existing DNS data of the turbulent channel flow accompanied by DR [18], as shown in Figure 3. The value of 𝜅_𝜂 is determined from the peak of the spanwise energy-dissipation spectrum function of 𝜅²_𝑧𝐸_𝑢𝑢(𝜅_𝑧). For a viscoelastic flow at the friction Weissenberg number Wi_𝜏=30, the 𝜅_𝜂 is remarkably decreased at every height, especially in the near-wall region, compared to the case of Wi_𝜏=10. As a result, the intersection of the 𝜅_𝐿 line and the 𝜅_𝜂 line moves to a higher position of about 𝑦⁺=12, as given in Figure 4. The 𝜅_𝜂 is prone to shift leftwards as the Weissenberg number is increasing consistent to the above Lumley’s drag-reduction model.

Figure 3 

Premultiplied one-dimensional spanwise energy spectra of streamwise velocity fluctuation in drag-reducing viscoelastic flow at Re𝜏=395. (a) Wi𝜏=10, 𝛽=0.8 and (b) Wi𝜏=30, 𝛽=0.5.

Figure 4 

(a) Dependence of the storage modulus and loss modulus on a frequency for the viscoelastic fluid of Giesekus model. (b) Scaling relation in the turbulent channel flow with/without drag reduction: shifting in the wave number of the energy-dissipative range by viscoelastic effects.

3.2. Dynamic Characterization of Viscoelastic Fluid

In this section, we investigate the dynamic characterization of the viscoelastic fluid, in order to discuss the effect of its viscoelasticity on turbulent fluctuations (eddies). A model fluid for this study is based on the Giesekus model [27], which is one of the typical viscoelastic models to describe the relationship between stress and strain exerted on the fluid of interest. The rheological constitutive equation with respect to the viscoelastic stress tensor 𝜏^𝑝_𝑖𝑗 is written as 𝜏^𝑝_𝑖𝑗+𝜆𝐷𝜏^𝑝_𝑖𝑗𝐷𝑡−𝜏^𝑝_𝑖𝑘𝜕𝑢_𝑗𝜕𝑥_𝑘−𝜕𝑢_𝑖𝜕𝑥_𝑘𝜏^𝑝_𝑘𝑗+𝛼𝜇_𝛼𝜏^𝑝_𝑖𝑘𝜏^𝑝_𝑘𝑗=𝜇₀𝜕𝑢_𝑖𝜕𝑥_𝑗+𝜕𝑢_𝑗𝜕𝑥_𝑖, (11) where 𝜆 is the relaxation time, 𝜇₀ is the zero shear rate viscosity of the solution, and 𝛼 is the mobility factor to indicate the level of nonlinearity.

It is generally known that a measurement under a small amplitude oscillatory flow is a powerful tool for describing the microscopic state of the test fluid. Therefore, we consider here a small fluid volume under a simple shear flow whose shear rate varies periodically with time as 𝛾̇=𝛾̇_maxexp(𝑖𝜅𝑡), (12) where 𝜅 is the frequency of the oscillatory shear stress. The dimensionless complex modulus as one of rheological material functions is obtained from 𝑐𝐺=¹²𝛾̇_max, (13) where the nondimensional conformation tensor 𝑐₁₂ is derived explicitly from the viscoelastic stress tensor; that is, 𝜏^𝑝_𝑖𝑗=𝜇₀(𝑐_𝑖𝑗−𝛿_𝑖𝑗)/𝜆. Note that the complex shear modulus 𝐺 consists of the dynamic storage modulus 𝐺^′ and the dynamic loss modulus 𝐺^′′ as follows: 𝐺=𝐺^′(𝜅)+𝑖𝐺^′′(𝜅). (14) In the case of an elastic body, there is no phase difference between strain and stress, while in a viscous body the phase difference is 𝜋/2. Therefore, 𝐺^′ and 𝐺^′′ can be used as barometers of elasticity and viscosity, respectively.

Typical 𝐺^′ and 𝐺^′′ versus frequency data for the viscoelastic fluid described by the Giesekus model are shown in the top panel of Figure 4, where the horizontal axis represents the wave number 𝜅. Both moduli clearly depend only on the Weissenberg number, and they would not depend on any other flow parameter, such as the Reynolds number and the maximum shear rate 𝛾_max, for the laminar simple shear flow. Here, the Weissenberg number Wi_𝜉 is defined as the product of 𝜆 and 𝛾_max. The most obvious consequence of the increasing Weissenberg number on the frequency sweep in the figure is to decrease the crossover frequency, at which 𝐺^′=𝐺^′′, thereby extending the high plateau in 𝐺^′(𝜅) to lower frequencies. As Wi_𝜉 varies from 11 to 30, the frequency of the crossover point shifts from 𝜅⁺₁≈0.090 to 𝜅⁺₂≈0.033. This shift factor corresponds to the ratio of the two Weissenberg numbers; that is, Δ=𝜅₂/𝜅₁=11/30. The magnitude of the plateau modulus is essentially unchanged against the Weissenberg number.

From the above facts, it can be conjectured that, as the Weissenberg number increases, the behavior of relevant fluid becomes elastic against high-frequency inputs associated with small-scale eddies in the near-wall turbulence. Correspondingly, the relationship between the energy-dissipation range and the frequency response of the Giesekus model reveals that the shift in wave number of the energy dissipation is consistent with the change of viscoelastic behavior, as shown in Figure 4. From this result, it can be said that viscoelastic fluid should behave more elastically against high-frequency eddy motions in turbulence and suppress them.

4. Model Development

While some details of the present model differ from Jones and Launder [19], their low-Reynolds-number 𝑘-𝜀 model will be applied unaltered in most respects. As for the effect of viscoelastic fluid, the above observation has driven us to consider the time-scale ratio between the characteristic time scale of turbulent eddy motions and the relaxation time 𝜆 of the relevant fluid, and we employ it to describe the diminished values of 𝜈_𝑡. In the drag-reducing viscoelastic fluid that has a relaxation time comparable to some scale of turbulent eddies, the smaller eddies, if larger than the Kolmogorov scale, should be damped due to the elastic behavior of the fluid. Using the Kolmogorov time 𝑡_𝜂 as the characteristic time scale in turbulence, the time-scale ratio can be proposed as 𝜆𝑡_𝜂=𝜆𝜈_𝑠=/𝜀Wi_𝜏𝛽Re_𝜏/𝜀^∗, (15) where 𝜈_𝑠, Re_𝜏, and 𝛽 are defined as the solvent kinematic viscosity (𝜇_𝑠/𝜌), the Reynolds number based on the friction velocity and the channel half width, and the ratio of the solvent viscosity 𝜇_𝑠 to the zero-shear-rate viscosity of the solution 𝜇₀, respectively. Then, replacing (2), let us modify 𝜈_𝑡 by introducing an additional damping function 𝑓_𝑣 based on this time-scale ratio, as follows: 𝜈_𝑡=𝐶_𝜇𝑓_𝜇𝑓_𝑣𝑘²𝜀̃,𝑓_𝑣=exp−𝜓₀Wi_𝜏𝛽Re_𝜏/𝜀^∗̃1−𝛽𝛽^0.1. (16)(17) Here, the “homogeneous” dissipation rate 𝜀̃(=𝜀−𝜀)̂ is introduced, while 𝜀̂=2𝜈_𝑠𝜕√𝑘𝜕𝑦²,when𝜕√𝑘𝜕𝑦≥0,0,when𝜕√𝑘𝜕𝑦<0, (18) is defined with an assumption that 𝑦 is normal to the wall. To accommodate the effect of 𝛽 on DR, the form of (𝛽⁻¹−1)^0.1 was chosen ad hoc and adopted to 𝑓_𝑣. The presently optimized value of the additional model coefficient 𝜓₀ is 0.36.

The conventional damping function 𝑓_𝜇, that is used for the low-Reynolds-number turbulence model to comply with the damping effect of the wall, remains the same as for a Newtonian fluid. For 𝑓_𝜇 in (15), we adopt the function proposed by Kawashima and Kawamura [28] as follows: 𝑓_𝜇−𝑦=1.0−exp^𝜂−𝑦28520^𝜂3, (19) where 𝑦_𝜂 is the wall-normal distance normalized by the Kolmogorov length scale, 𝜂=(𝜈³_𝑠/𝜀)^1/4. The model of 𝜈_𝑡 by (17) and (19) satisfies the following two assumptions: first, very small-scale eddies are damped by the elastic behavior of the fluid and thereby the eddy viscosity decreases (𝜈_𝑡→0 when 𝑡_𝜂≪𝜆); secondly, the eddy viscosity for less elastic fluid is close to be that for the Newtonian case (𝜈_𝑡→ (2) when 𝑡_𝜂≫𝜆).

The following equations are the Reynolds equation and the ensemble-averaged constitutive equation, which govern the incompressible viscoelastic-fluid flow as follows: 𝐷𝑈_𝑖1𝐷𝑡=−𝜌𝜕𝑃𝜕𝑥_𝑖+𝜈_𝑠𝜕²𝑈_𝑖𝜕𝑥²_𝑗−𝜕𝑢_𝑖𝑢_𝑗𝜕𝑥_𝑗+𝜈_𝑎𝜆𝜕𝐶_𝑖𝑗𝜕𝑥_𝑗,𝐷𝐶_𝑖𝑗+𝜕𝐷𝑡𝑐_𝑖𝑗𝑢_𝑘𝜕𝑥_𝑘−𝐶_𝑖𝑘𝜕𝑈_𝑗𝜕𝑥_𝑘−𝜕𝑈_𝑖𝜕𝑥_𝑘𝐶_𝑘𝑗−𝜕𝑢_𝑖𝜕𝑥_𝑘𝑐_𝑘𝑗+𝑐_𝑖𝑘𝜕𝑢_𝑗𝜕𝑥_{𝑘Λ𝑖𝑗}+1𝜆𝐶_𝑖𝑗−𝛿_𝑖𝑗+𝛼𝐶_𝑖𝑘𝐶_𝑘𝑗+𝑐_𝑖𝑘𝑐_𝑘𝑗−2𝐶_𝑖𝑗+𝛿_𝑖𝑗=0. (20)(21) In (21), the terms which include the 𝑐_𝑖𝑗, which represent the fluctuating part of the conformation tensor, should be modeled. The second term and the term of 𝑐_𝑖𝑘𝑐_𝑘𝑗 (included in the last term) in the left-hand side of (21) are small enough to be ignored. The term labeled as Γ_𝑖𝑗 in the left-hand side is modeled in the following form: Λ_𝑖𝑗≡𝑐_𝑖𝑘𝜕𝑢_𝑗𝜕𝑥_𝑘+𝜕𝑢_𝑖𝜕𝑥_𝑘𝑐_𝑘𝑗=𝜓₁1𝑡_𝜂𝜆𝜈_𝑎𝑢_𝑖𝑢_𝑗+𝜓₂1𝑡_𝜂𝐶_𝑖𝑗−𝛿_𝑖𝑗, (22) where 𝜓₁=0.001(1−exp(−𝑦_𝜂/100)) and 𝜓₂=−0.001(1−exp(−𝑦_𝜂/5)) are given based on the DNS data [18]. This form was originally proposed by Leighton et al. [11].

From the Navier-Stokes equation, the transport equation of turbulent kinetic energy is obtained and modeled as 𝐷𝑘𝐷𝑡=𝐷_𝑘+𝑃_𝑘−(𝜀̃+𝜀)̂+𝜋_𝑘+𝐸^𝑑_𝑘+𝐸^𝜀_𝑘,𝐷_𝑘=𝜕𝜕𝑥_𝑗𝜈_𝑠+𝜈_𝑡𝜎_𝑘𝜕𝑘𝜕𝑥_𝑗,𝑃_𝑘=−𝜈_𝑡𝜕𝑈_𝑖𝜕𝑥_𝑗+𝜕𝑈_𝑗𝜕𝑥_𝑖+23𝛿_𝑖𝑗𝑘𝜕𝑈_𝑖𝜕𝑥_𝑗,𝜋_𝑘1=−2𝜈_𝑠𝜕𝜕𝑥_𝑘𝑘𝜀𝜀̂𝜀²𝜕𝜀̂𝜕𝑥_𝑘,𝐸^𝑑_𝑘=𝜕𝜕𝑥_𝑘𝜈_𝑎𝜆𝑢_𝑖𝑐_𝑖𝑘𝐸≈0,_𝜀𝑘𝜈=−^𝑎𝜆𝑐_𝑖𝑘𝜕𝑢_𝑖𝜕𝑥_𝑘𝜈=−^𝑎Λ2𝜆^𝑖𝑖, (23) and the transport equation of the dissipation of turbulent kinetic energy is written as 𝐷𝜀̃𝐷𝑡=𝐷_𝜀+𝑃_𝜀−𝑌+𝐸+𝐸_𝜀=𝜕𝜕𝑥_𝑗𝜈_𝑠+𝜈_𝑡𝜎_𝜀𝜕𝜀𝜕𝑥_𝑗+𝜀̃𝑘𝐶_𝜀1𝑃_𝑘−𝐶_𝜀2𝜀̃−𝜀̂+𝐸+𝐸_𝜀,𝐸=𝐶_𝜀3𝜈_𝑠𝜈_𝑡𝜕²𝑈𝜕𝑥²_𝑘²,𝐸_𝜀𝜈=2^𝑠𝜈_𝑎𝜆𝜕𝑢_𝑖𝜕𝑥_𝑘𝜕𝜕𝑥_𝑘𝜕𝑐_𝑖𝑗𝜕𝑥_𝑗≈0, (24) where 𝐸^𝑑_𝑘 and 𝐸_𝜀 are assumed to be negligible compared to the other terms. Numerical studies of these terms showed them indeed to be small [7, 18]. Their omission in this model simplified the formulation of a stable numerical scheme. We can obtain the value of 𝐸^𝜀_𝑘 by applying (22). Following Kawashima and Kawamura [28], the model constants in the above equations are taken as 𝐶_𝜀1=1.44, 𝐶_𝜀2=1.92, 𝐶_𝜀3=0.6, and 𝜎_𝑘=𝜎_𝜀−𝑦=1.0−0.5exp^𝜂20². (25)

Prior to discussion on the feasibility of the present model for the drag-reducing viscoelastic fluid, Figure 5 shows the computed budget of the turbulent kinetic energy for the Newtonian fluid. The results of DNS by Abe et al. [22] are also shown for comparison. The agreement can be regarded as excellent everywhere in the turbulent channel flow at Re_𝜏=395.

Figure 5 

Budget of turbulent kinetic energy for Newtonian fluid in turbulent channel flow at Re𝜏=395 with emphasis on the near-wall region.

5. Results

Results are presented below for drag-reduced turbulent channel flows. The tested range of the flow and rheological parameters were set to be Re_𝜏=150 and 395, 𝛽=0.3–0.8, 𝛼=0.001, and Wi_𝜏=10–40. Note that we tested four different values of the Weissenberg number (10, 11, 30, and 40), at which Tsukahara et al. [18] executed a series of DNS at the same Reynolds number. The numerical solutions of the proposed model equations were obtained using the second-order central-difference scheme for discretization, while convergence calculation was performed by the successive over-relaxation (SOR) method. We used nonuniform grids, which divide the half length of channel into 128 parts. The wall boundary conditions are 𝑈=𝑘=𝜀̃=0 and the zero gradient boundary condition is applied to all variables at the channel center (except 𝑉=0). In this paper, the superscript (⁺) describes a nondimensionalized quantity by the friction velocity 𝑢_𝜏, the effective viscosity [1], and the density of solution 𝜌, while (^∗) describes a nondimensionalized quantity by the channel half-width 𝛿, 𝑢_𝜏, and 𝜌. The friction Reynolds and Weissenberg numbers to be used in the following are defined as Re_𝜏=𝜌𝑢_𝜏𝛿/𝜇₀ and Wi_𝜏=𝜌𝑢²_𝜏𝜆/𝜇₀, respectively.

As an illustration of the performance of the present model, we present some results of the predicted rate of DR in several cases, in which the DNS database is available, see Table 1. Here, the percent drag reduction is determined as DR𝐶%=1−^𝑓visc𝐶_𝑓Newt×100(%), (26) where the suffixes “visc” and “Newt” stand for the friction-coefficient values in the viscoelastic and Newtonian flows, respectively, at the same bulk Reynolds number. The model predicts DR%and its parameter dependence quite well in the region from low to high DR flows, as shown in Table 1. At both Reynolds numbers, the parameter combination of (Wi_𝜏,𝛽)=(11,0.5) results in quite low DR%, while the DNS demonstrated moderate DR. This defect occurs in very low Weissenberg numbers with a low value of 𝛽 and, moreover, the model has predicted that a near-zero Wi_𝜏 would produce a negative DR%, as shown later. This issue is not serious, because such a combination of low Wi_𝜏 and low 𝛽 should not be impractical rheological condition; usually, dilute polymer solutions should be equivalent to viscoelastic fluids with low Wi_𝜏 and high 𝛽 (≈1).

Table 1 

Percentage drag reduction DR%.

5.1. Mean Flow Statistics

Figure 6 shows the mean velocity profiles, as a function of 𝑦⁺, computed using the present model together with DNS data. Note again that, for the drag-reduced flows, 𝑦⁺ is normalized by the effective viscosity of each case, whose values are described in Table 2. Also shown in Figure 6 is the MDR asymptote, where the log-law profile for maximum drag reduction is that proposed by Virk [26]. At each Weissenberg number, the present model was in good agreement with the DNS data. Especially, the model applied for Wi_𝜏=10 and 30 exhibited a better performance as 𝑦 is close to the channel center. Small discrepancies are noticeable in the elastic layer and at the channel center, as Wi_𝜏 increases to 40, at which the achieved drag-reduction rate is close to MDR. However, in the elastic layer, the present model successfully demonstrates the profile that coincides with the MDR asymptote. The vertical displacements of the velocity profiles for high Weissenberg numbers are consistent with the DNS and the aforementioned results of DR%.

Table 2 

Effective wall viscosity normalized by 𝜇₀.

(a) Re𝜏 = 150

(b) Re𝜏 = 395

Figure 6 

Mean velocity profiles in wall coordinates. Comparison between the model prediction and the corresponding DNS data [18].

From the above assessments, it can be said that the Weissenberg-number dependence of the mean velocity is well described by the present model in a range of parameters including high DR cases, except for the combination of low Wi_𝜏 and low 𝛽. The model provided an underestimated velocity for Wi_𝜏=11 and 𝛽=0.5, as shown in Figure 6(b), while the result of a set of Wi_𝜏=10 and 𝛽=0.8 seems to be predicted well, as given in Figure 6(a).

Figures 7 and 8 compare the predictions of 𝐶_𝑓-Re_𝑚 and DR%-Re_𝑚 with data from DNS [18]. The computations were made in a range of Wi_𝜏 from low DR to the maximum DR for the three values of 𝛽 (=0.3, 0.5, and 0.8), at two Reynolds numbers. For reference, the MDR asymptote of Virk [26] is also plotted. As seen in the figures, the Weissenberg number dependences of 𝐶_𝑓 and DR%are demonstrated, enabling us to estimate the maximum Wi_𝜏 that gives rise to the maximum DR. For example, according to the model predictions at Re_𝜏=395, the MDR would be achieved at Wi_𝜏≈45, 50, and 60 for 𝛽=0.3, 0.5, and 0.8, respectively. However, it is not necessarily the case that those critical states of MDR should be corresponding to the limiting conditions of the model prediction. The present model may provide a result with both lower 𝐶_𝑓 and higher DR%than those for MDR as the Weissenberg number increases further (not shown in Figures 7 and 8). With too high Wi_𝜏, the mean velocity obtained by the model reveals that it exceeds the profile for the MDR asymptote [26] and, in addition, the Reynolds shear stress of −𝑢𝑣 still remains to be nonzero value. This result argues against the experimental fact that the contribution of −𝑢𝑣 would be almost absent in the highly drag-reduced turbulent channel flows by polymer/surfactant [4, 29, 30]. The reproducibility of the MDR state by means of RANS would be an interesting and challenging issue.

Figure 7 

Frictional coefficient.

Figure 8 

Drag reduction rate.

As can be seen in Figures 9 and 10, the predicted 𝑘⁺ and 𝐶₁₁ follow those trends in the DNS data, although the near-wall peak value reveals clearly 𝑘⁺ to be overestimated. The discrepancy becomes large as the Weissenberg number increases. As for the wall-normal height of the peak, the model shows almost good agreements with the DNS data. It would be difficult to improve this aspect in the framework of the isotropic 𝑘-𝜀 model, because highly drag-reduced viscoelastic flows should be inevitably accompanied by enhanced anisotropy in the near-wall region.

(a) Re𝜏 = 150

(b) Re𝜏 = 395

Figure 9 

Same with Figure 6, but with respect to turbulent kinetic energy, 𝑘+.

(a) Re𝜏 = 150

(b) Re𝜏 = 395

Figure 10 

Same with Figure 6, but with respect to a conformation normal stress, 𝐶+11.

5.2. Reynolds and Conformation Stresses

The total shear stress can be obtained as 𝜏⁺_total=𝛽Re_𝜏𝜕𝑈⁺𝜕𝑦^∗−𝑢⁺𝑣⁺+1−𝛽Wi_𝜏𝐶₁₂. (27) Here, the terms are (in order from left to right of the right-hand side) the viscous shear stress, the Reynolds shear stress, and the viscoelastic shear stress. Figure 11 shows the wall-normal distributions of these stresses in several conditions chosen from Table 1. In all the test cases, the obtained total shear stress is in accordance with the theoretical value as follows: 𝜏⁺_total=1−𝑦^∗, (28) proving the validity of the model calculation.

(a) (Re𝜏, Wi𝜏, 𝛽) = (150, 10, 0.8)

(b) (Re𝜏, Wi𝜏, 𝛽) = (150, 30, 0.5)

(c) (Re𝜏, Wi𝜏, 𝛽) = (150, 40, 0.5)

(d) (Re𝜏, Wi𝜏, 𝛽) = (395, 30, 0.5)

Figure 11 

Shear stress profiles of drag-reducing viscoelastic flow in channel.

The viscous shear stress is well represented by the present model throughout the entire channel, as shown in Figure 11. All the shear stresses for Wi_𝜏=10 are considerably close to the DNS data, indicating that nearly Newtonian fluid calculations would be also executed adequately, as can be seen in Figure 11(a). However, the other three cases shown in the figure have resulted in the overestimated peaks of  −𝑢⁺𝑣⁺ and their shifts towards the channel center, compared with those by the DNS. As for the viscoelastic stress, the model exhibits tendencies to overestimate it in 10<𝑦⁺<70 in the buffer and elastic layers and to underestimate it significantly in the outer layer. The model does, however, correctly show the dependence of the drag-reduction rate on the Weissenberg number, as given in Table 1.

To obtain a better agreement with DNS, it is necessary to further improve the modeling of the viscoelastic contribution in the region above the buffer layer for high Weissenberg numbers, where the drag-reduction rate is being close to the maximum drag reduction. A typical result about the profiles of 𝐶_𝑖𝑗 computed by the present model is shown in Figure 12. As noted above with respect to 𝐶₁₂, the comparison with DNS data shows the present model to underestimate the conformation stresses throughout the channel except for the near-wall region.

Figure 12 

Mean conformation tensor for Re𝜏=395, Wi𝜏=30, and 𝛽=0.5.

Figure 13 shows the budget in the 𝑘 transport equation of the viscoelastic fluid flow at Wi_𝜏=30 and 𝛽=0.5, with a moderate DR, for both the present model and the DNS data. The level of the production term calculated by the model agrees well with that of the DNS data, both in terms of its magnitude and the general shape of the curve. There is a similar agreement for the molecular diffusion, but its negative peak at 𝑦⁺≈18 is slightly overestimated. Due to both this overestimation and the false negative peak in the pressure diffusion, the turbulent diffusion is enhanced considerably around 𝑦⁺=15. The dissipation rate 𝜀 is in good agreement near the wall, but in the other region it is two or three times too high. This should be a deficiency due to the characteristic of the 𝑘-𝜀 model to be applied to the strongly anisotropic turbulence, as similar to that observed already in the overprediction of 𝑘. The viscoelastic contribution 𝐸^𝜀_𝑘 is found to be inadequately predicted as being very small compared to the other terms, This also might cause the overestimated value of 𝜀 that would compensate the unbalance (in the budget of 𝑘) due to too small 𝐸^𝜀_𝑘.

Figure 13 

Budget of turbulent kinetic energy for drag-reducing viscoelastic fluid in turbulent channel flow at Re𝜏=395, Wi𝜏=30, and 𝛽=0.5.

5.3. Damping Function

It is well known that DR is associated with a decrease in the Reynolds shear stress [4]. By adopting (16), this reduction in the Reynolds shear stress can be achieved through a decrease in 𝑘, an increase in 𝜀̃, and/or a decrease of 𝑓_𝑣. Since, as shown already, both variations of 𝑘 and 𝜀̃ might not provide a significant cause of the reduction, the DR predicted well by the model must be accompanied by a large decrease of the damping function 𝑓_𝑣. Figure 14 shows the damping function, the product of 𝑓_𝜇 and 𝑓_𝑣, both for the present model and the estimated value via DNS, indicating that the present model has yielded fairly close agreement with the DNS result. It can be confirmed that the near-wall asymptotic behavior of the product of the damping functions has been well predicted via the present model. In terms of the Reynolds shear stress (Figure 11), the deviations from the DNS results are detected in the outer layer, that is, 𝑦/𝛿>0.1-0.2, and probably attributed to the inaccuracy (overestimation) of the turbulent kinetic energy in the highly turbulent region.

(a) Re𝜏 = 150

(b) Re𝜏 = 395

Figure 14 

Validation of damping function, the product of 𝑓𝜇 and 𝑓𝑣.

6. Conclusion

A new low-Reynolds-number 𝑘-𝜀 model including an additional damping function in the eddy viscosity was applied to predict the viscoelastic flow accompanied by turbulent drag reduction. The proposed model was constructed, based on the energy-dissipative range and the dynamic characterization of the viscoelastic fluid, and was tested by comparison with the DNS data for the drag-reduced channel flow.

The analysis of frequency response of the Giesekus model revealed that, for a high Weissenberg number, the elastic behavior of fluid should have influence on the diminishing turbulent eddy motions. Consequently, with the increasing relaxation time of the fluid, the dissipation scale of turbulent kinetic energy is increased and the viscous sublayer is thickened. The mean velocity and the drag-reduction rate were well reproduced by the present model (at least for the present parameters). The present model, however, showed some significant deviations from the DNS data with respect to the Reynolds shear and viscoelastic stresses, when applied for high Weissenberg numbers.

The most serious deficiency found in the present comparisons with DNS data was the overpredicted dissipation rate, which compensated the underestimated term pertaining to the viscoelastic contribution. This error in the overpredictions of the dissipation rate as well as the turbulent kinetic energy would occur in other situations relating to strongly anisotropic turbulence. One cause of this may be the insufficient model for the viscoelastic contribution in the 𝑘 transport equation. It seems that the deficiency is associated with also the inadequacy in the isotropic 𝑘-𝜀 model for highly drag-reducing viscoelastic flows, although the drag-reduction phenomena in terms of bulk flows were predicted well by the present model. Hence, the model predictions should be improved further also in the framework of other models including the anisotropic/nonlinear 𝑘-𝜀 model and the Reynolds stress model. In addition, the reproducibility of the MDR state and the application of the present model to practical flows and situations will be investigated in the near future.

In this paper, our discussion is restricted to the case of 𝛼=0.001. The present dumping function model of (17) takes into account the inertia (through the Reynolds number), the viscoelastic relaxation effect (through the dependence on the Weissenberg number), and the polymer/surfactant concentration in solution (through the viscosity ratio). However, any dependence on the maximum polymer effect on extensional thickening has not been explicitly involved in the model formulation. For the Giesekus model used here, the dependence on the shear thinning/thickening is controlled through the mobility factor 𝛼. Only one value of that parameter has been used throughout this work because of the limited DNS database. Although not shown in this paper, some comparative verifications with DNS data at different values of 𝛼 [31, 32] were carried out. We confirmed the present model would work well and drew the same conclusions, at least, in the range of 𝛼=0.001–0.01. We should consider the key effect depending on 𝛼further to improve the model applicable to dilute-to-dense solutions.

Appendix

For fully-developed channel flows, some terms including the spatial derivatives in the horizontal directions can be neglected. In the following, the reduced equations for obtaining time-averaged values in the channel flow are given with being nondimensionalized by the channel half width 𝛿, the friction velocity 𝑢_𝜏=−(𝛿/𝜌)𝜕𝑃/𝜕𝑥, and the solution zero-shear-rate kinematic viscosity 𝜈.

The momentum equation 𝛽1+Re_𝜏𝜕²𝑈⁺𝜕𝑦^∗2−𝜕𝑢⁺𝑣⁺𝜕𝑦^∗+1−𝛽Wi_𝜏𝜕𝐶₁₂𝜕𝑦^∗=0, (A.1)

the constitutive equation 𝐶₁₁=1()×)𝐶1−2𝛼(1−𝛼−𝛼_1𝑘𝐶_𝑘1+Wi_𝜏Re_𝜏2𝐶₁₂𝜕𝑈⁺𝜕𝑦^∗+Λ^∗₁₁,𝐶₁₂=1()𝐶1−2𝛼−𝛼_1𝑘𝐶_𝑘2+Wi_𝜏Re_𝜏𝐶₂₂𝜕𝑈⁺𝜕𝑦^∗+Λ^∗₁₂,𝐶₁₃=1()𝐶1−2𝛼−𝛼_1𝑘𝐶_𝑘3+Wi_𝜏Re_𝜏𝐶₂₃𝜕𝑈⁺𝜕𝑦^∗+Λ^∗₁₃,𝐶₂₂=1())𝐶1−2𝛼(1−𝛼−𝛼_2𝑘𝐶_𝑘2+Wi_𝜏Re_𝜏Λ^∗₂₂,𝐶₂₃=1()𝐶1−2𝛼−𝛼_2𝑘𝐶_𝑘3+Wi_𝜏Re_𝜏Λ^∗₂₃,𝐶₃₃=1())𝐶1−2𝛼(1−𝛼−𝛼_3𝑘𝐶_𝑘3+Wi_𝜏Re_𝜏Λ^∗₃₃, (A.2) where Λ^∗_𝑖𝑗=𝜓₁Wi_𝜏1−𝛽𝛽/Re_𝜏𝜀^∗̃𝑢⁺_𝑖𝑢⁺_𝑗+𝜓₂𝛽/Re_𝜏𝜀^∗̃𝐶_𝑖𝑗−𝛿_𝑖𝑗, (A.3) the transport equation of 𝑘, normalized by 𝛿/𝑢³_𝜏, 𝜕𝜕𝑦^∗𝛽Re_𝜏+𝜈^∗_𝑡𝜎_𝑘𝜕𝑘⁺𝜕𝑦^∗+𝜈^∗_𝑡𝜕𝑈⁺𝜕𝑦^∗2−𝜀^∗̃+𝜀^∗̂−12𝛽Re_𝜏𝜕𝜕𝑦^∗×𝑘⁺𝜀^∗𝜀^∗̂𝜀^∗2𝜕𝜀^∗̂𝜕𝑦^∗−1−𝛽2Wi_𝜏Λ_𝑖𝑖=0, (A.4) and the transport equation of 𝜀^∗̃, normalized by 𝛿²/𝑢⁴_𝜏, 𝜕𝜕𝑦^∗𝛽Re_𝜏+𝜈^∗_𝑡𝜎_𝜀𝜕𝜀^∗𝜕𝑦^∗+𝜀^∗̃𝑘⁺𝐶_𝜀1𝜈^∗_𝑡𝜕𝑈⁺𝜕𝑦^∗2−𝐶_𝜀2𝜀̃−𝜀̂+𝐶_𝜀3𝛽Re_𝜏𝜈^∗_𝑡𝜕²𝑈⁺𝜕𝑦^∗22=0. (A.5) The Reynolds stress 𝑢⁺_𝑖𝑢⁺_𝑗 is given by (1) and (16).

Nomenclature

𝐶𝑓:	Skin friction coefficient, =2𝜏𝑤/[𝜌(∫10𝑢𝑑𝑦∗)2]
𝐶𝑖𝑗:	Mean conformation tensor
𝑐𝑖𝑗:	Fluctuation of conformation tensor
𝐶𝜇, 𝐶𝜀𝑖:	Model coefficients
DR%:	Percentage of drag reduction
𝑓𝑣, 𝑓𝜇:	Damping functions
𝐺:	Complex shear modulus
𝐺′:	Dynamic storage modulus
𝐺′′:	Dynamic loss modulus
𝑘:	Turbulent kinetic energy, =𝑢𝑖𝑢𝑖/2
𝑝:	Pressure
Re𝑚:	Bulk Reynolds number, =2𝜌(∫10𝑢𝑑𝑦∗)𝛿/𝜇0
Re𝑡:	Turbulence Reynolds number, =𝜌𝑘2/(𝜇0𝜀)̃
Re𝜏:	Friction Reynolds number, =𝜌𝑢𝜏𝛿/𝜇0
𝑡:	Time
𝑡𝜂:	Kolmogorov time scale
𝑈𝑖:	Mean velocity vector
𝑢𝑖:	Fluctuation of velocity vector
𝑢𝜏:	Friction velocity, =𝜏𝑤/𝜌
Wi𝜉:	Weissenberg number, =𝜆𝛾̇max
Wi𝜏:	Friction Weissenberg number, =𝜌𝑢2𝜏𝜆/𝜇0
𝑥:	Streamwise coordinate
𝑦:	Wall-normal coordinate
𝑧:	Spanwise coordinate
𝛼:	Anisotropic mobility factor
𝛽:	Ratio of shear viscosities, =𝜇𝑠/𝜇0
𝛾̇:	Shear rate
𝛿:	Channel half width
𝜀:	Dissipation rate of 𝑘
𝜂:	Kolmogorov length scale
𝜅:	Wave number or frequency
𝜅𝐾:	Wave number or frequency
𝜆:	Relaxation time
𝜇0:	Solution zero-shear-rate shear viscosity, =𝜇𝑠+𝜇𝑎
𝜇𝑎, (𝜈𝑎):	Shear (kinematic) viscosity of additive contribution
𝜇𝑠, (𝜈𝑠):	Shear (kinematic) viscosity of solvent contribution
𝜇𝑡, (𝜈𝑡):	Turbulent (kinematic) viscosity
𝜌:	Density
𝜎𝑘, 𝜎𝜀:	Model coefficients
𝜓𝑖:	Model functions
𝜏𝑖,𝑗:	Viscoelastic stress tensor
𝜏𝑤:	Mean wall shear stress
()∗:	Normalized by 𝑢𝜏 and 𝛿
()+:	Normalized by 𝑢𝜏 and the effective viscosity
():	Statistically averaged.

Acknowledgments

The authors are grateful to Professor H. Kawamura (Tokyo University of Science, Suwa) for fruitful discussions and comments, and to Mr. H. Kawamoto (who was a master’s course student at Tokyo University of Science) for his work. The DNS database we utilized here to assess results predicted by the present model was acquired using supercomputing resources at the Cyberscience Center at Tohoku University. This work was partially supported by Grants-in-Aid for Scientific Research (C) no. 25420131 from Japan Society for the Promotion of Science (JSPS).