# Proposal of Damping Function for Low-Reynolds-Numberk-ε Model Applicable in Prediction of Turbulent Viscoelastic-Fluid Flow

Authors: Takahiro Tsukahara; Yasuo Kawaguchi
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

A low-Reynolds-numberk-ε 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.


---

## Body

## 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 k-ε 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 Λij 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 v2-f 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 k-ω 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 k 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 numberk-ε model of Newtonian-fluid flow (here, k and ε are the turbulent kinetic energy and its dissipation rate, resp.). Two major corrections to be applied to the low Reynolds number k-ε 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 thek and ε equations are then developed. Finally, the proposed model will be tested in comparison with DNS results.


## 2. Introduction of Low-Rek-ε 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 vectorxi=(x1,x2,x3)≡(x,y,z) denotes the coordinate directions; ui=(u1,u2,u3)≡(u,v,w) 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 numberk-ε model for a turbulent flow relates the Reynolds stress to the mean strain field as follows:
(1)-uiuj¯=νt(∂Ui∂xj+∂Uj∂xi)-23kδij,
where δij is the Kronecker delta and the turbulent kinematic viscosity is expressed as
(2)νt=Cμfμk2ε~.
Here, Cμ is a constant with value equal to 0.09. A number of proposals for the damping function fμ have been made; for instance, Jones and Launder [19] functionalized fμ using a turbulent Reynolds number Ret as follows:
(3)fμ=exp(-2.51+Ret/50).
On the other hand, some researchers introduced the wall-normal distance normalized in wall unit, y+=yρuτ/μ0, as a parameter of fμ. A typical example is the one proposed by Nagano and Tagawa [20] as follows:
(4)fμ={1-exp(-y+26)}{1+4.1Ret3/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 fμ 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 functionfμ  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 offμ for viscoelastic flow is necessarily required. Here, we have investigated the predictive values of fμ 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 fμ 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
(5)fμ=fμv·fμe={1-[1+|1-n1+n|y+]-(1/26.5)|(1+n)/(1-n)|}·{1-[1+|q-13-q|y+Cn(1-q)/(2-q)]-(1/26.5)|(3-q)/(q-1)|},
where n is the viscosity power index of the shear behavior, q is the viscosity power index of the Trouton-ratio behavior, and Cn 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 is2δ. The one-dimensional energy spectrum of the streamwise velocity fluctuation is defined as
(6)∫0∞‍Euu(κ)dκ=uu¯,
where κ is a wave number. In the logarithmic layer, there is no relevant length other than the wall-normal height y according to the mixing length theory [24] and, hence, the energy-containing large-scale eddies should be scaled roughly with y. With respect to the energy spectrum, its peak is expected to occur at a wave number κL given by
(7)κL·y≃1,
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 κi2Euu(κ) [25]. Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonian-fluid turbulence would occur at
(8)κηη≃0.2,
where η=(ν3/ε)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
(9)ε≃uτ3κK·y,
using the von Kármán constant κK and the friction velocity uτ [24]. Therefore, the relation between y and κη is rewritten from (8) as
(10)κη+(κK·y+)1/4≃0.2.
The size of turbulent fluctuations (or eddies) ranges between κL and κη. Figure 2 shows both wave numbers of κL and κη as a function of y+ in the ordinate. If κK=0.4, the lines of (7) and (10) intersect at y+=6.3. Here, we call this point as A. Above the point A, the momentum transport by turbulent vortical motions is active. As y+ 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κL and κη is expected to be changed as shown with broken lines in Figure 2, according to Lumley [6]. He reasoned that less difference between κL+ and κη+ as a result of decreased κη+ might lead to the expanding of the viscous sublayer. The intersection of κL+ 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 κz2Euu(κz). 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 κL line and the κη line moves to a higher position of about y+=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 atReτ=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 τijp is written as
(11)τijp+λ{DτijpDt-τikp∂uj∂xk-∂ui∂xkτkjp+αμα(τikpτkjp)}=μ0(∂ui∂xj+∂uj∂xi),
where λ is the relaxation time, μ0 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(12)γ˙=γ˙maxexp(iκt),
where κ is the frequency of the oscillatory shear stress. The dimensionless complex modulus as one of rheological material functions is obtained from
(13)G=c12γ˙max,
where the nondimensional conformation tensor c12 is derived explicitly from the viscoelastic stress tensor; that is, τijp=μ0(cij-δij)/λ. Note that the complex shear modulus G consists of the dynamic storage modulus G′ and the dynamic loss modulus G′′ as follows:
(14)G=G′(κ)+iG′′(κ).
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, G′ and G′′ can be used as barometers of elasticity and viscosity, respectively.TypicalG′ and G′′ 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 G′=G′′, thereby extending the high plateau in G′(κ) to lower frequencies. As Wiξ varies from 11 to 30, the frequency of the crossover point shifts from κ1+≈0.090 to κ2+≈0.033. This shift factor corresponds to the ratio of the two Weissenberg numbers; that is, Δ=κ2/κ1=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 Figure4. From this result, it can be said that viscoelastic fluid should behave more elastically against high-frequency eddy motions in turbulence and suppress them.


## 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 is2δ. The one-dimensional energy spectrum of the streamwise velocity fluctuation is defined as
(6)∫0∞‍Euu(κ)dκ=uu¯,
where κ is a wave number. In the logarithmic layer, there is no relevant length other than the wall-normal height y according to the mixing length theory [24] and, hence, the energy-containing large-scale eddies should be scaled roughly with y. With respect to the energy spectrum, its peak is expected to occur at a wave number κL given by
(7)κL·y≃1,
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 κi2Euu(κ) [25]. Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonian-fluid turbulence would occur at
(8)κηη≃0.2,
where η=(ν3/ε)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
(9)ε≃uτ3κK·y,
using the von Kármán constant κK and the friction velocity uτ [24]. Therefore, the relation between y and κη is rewritten from (8) as
(10)κη+(κK·y+)1/4≃0.2.
The size of turbulent fluctuations (or eddies) ranges between κL and κη. Figure 2 shows both wave numbers of κL and κη as a function of y+ in the ordinate. If κK=0.4, the lines of (7) and (10) intersect at y+=6.3. Here, we call this point as A. Above the point A, the momentum transport by turbulent vortical motions is active. As y+ 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κL and κη is expected to be changed as shown with broken lines in Figure 2, according to Lumley [6]. He reasoned that less difference between κL+ and κη+ as a result of decreased κη+ might lead to the expanding of the viscous sublayer. The intersection of κL+ 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 κz2Euu(κz). 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 κL line and the κη line moves to a higher position of about y+=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 atReτ=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 τijp is written as
(11)τijp+λ{DτijpDt-τikp∂uj∂xk-∂ui∂xkτkjp+αμα(τikpτkjp)}=μ0(∂ui∂xj+∂uj∂xi),
where λ is the relaxation time, μ0 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(12)γ˙=γ˙maxexp(iκt),
where κ is the frequency of the oscillatory shear stress. The dimensionless complex modulus as one of rheological material functions is obtained from
(13)G=c12γ˙max,
where the nondimensional conformation tensor c12 is derived explicitly from the viscoelastic stress tensor; that is, τijp=μ0(cij-δij)/λ. Note that the complex shear modulus G consists of the dynamic storage modulus G′ and the dynamic loss modulus G′′ as follows:
(14)G=G′(κ)+iG′′(κ).
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, G′ and G′′ can be used as barometers of elasticity and viscosity, respectively.TypicalG′ and G′′ 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 G′=G′′, thereby extending the high plateau in G′(κ) to lower frequencies. As Wiξ varies from 11 to 30, the frequency of the crossover point shifts from κ1+≈0.090 to κ2+≈0.033. This shift factor corresponds to the ratio of the two Weissenberg numbers; that is, Δ=κ2/κ1=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 Figure4. 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 k-ε 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 νt. 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 tη as the characteristic time scale in turbulence, the time-scale ratio can be proposed as
(15)λtη=λνs/ε=WiτβReτ/ε*,
where νs, Reτ, and β are defined as the solvent kinematic viscosity (μs/ρ), the Reynolds number based on the friction velocity and the channel half width, and the ratio of the solvent viscosity μs to the zero-shear-rate viscosity of the solution μ0, respectively. Then, replacing (2), let us modify νt by introducing an additional damping function fv based on this time-scale ratio, as follows:
(16)νt=Cμfμfvk2ε~,(17)fv=exp[-ψ0WiτβReτ/ε~*(1-ββ)0.1].
Here, the “homogeneous” dissipation rate ε~(=ε-ε^) is introduced, while
(18)ε^={2νs(∂k∂y)2,when(∂k∂y)≥0,0,when(∂k∂y)<0,
is defined with an assumption that y is normal to the wall. To accommodate the effect of β on DR, the form of (β-1-1)0.1 was chosen ad hoc and adopted to fv. The presently optimized value of the additional model coefficient ψ0 is 0.36.The conventional damping functionfμ, 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 fμ in (15), we adopt the function proposed by Kawashima and Kawamura [28] as follows:
(19)fμ=1.0-exp[-(yη285)-(yη20)3],
where yη is the wall-normal distance normalized by the Kolmogorov length scale, η=(νs3/ε)1/4. The model of νt 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 (νt→0 when tη≪λ); secondly, the eddy viscosity for less elastic fluid is close to be that for the Newtonian case (νt→ (2) when tη≫λ).The following equations are the Reynolds equation and the ensemble-averaged constitutive equation, which govern the incompressible viscoelastic-fluid flow as follows:(20)DUiDt=-1ρ∂P∂xi+νs∂2Ui∂xj2-∂uiuj¯∂xj+νaλ∂Cij∂xj,(21)DCijDt+∂cijuk¯∂xk-Cik∂Uj∂xk-∂Ui∂xkCkj-(∂ui∂xkckj¯+cik∂uj∂xk¯)︸Λij+1λ[Cij-δij+α(CikCkj+cikckj¯-2Cij+δij)]=0.
In (21), the terms which include the cij, which represent the fluctuating part of the conformation tensor, should be modeled. The second term and the term of cikckj¯ (included in the last term) in the left-hand side of (21) are small enough to be ignored. The term labeled as Γij in the left-hand side is modeled in the following form:
(22)Λij≡cik∂uj∂xk¯+∂ui∂xkckj¯=ψ11tηλνauiuj¯+ψ21tη(Cij-δij),
where ψ1=0.001(1-exp(-yη/100)) and ψ2=-0.001(1-exp(-yη/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)DkDt=Dk+Pk-(ε~+ε^)+πk+Ekd+Ekε,Dk=∂∂xj[(νs+νtσk)∂k∂xj],Pk=-[νt(∂Ui∂xj+∂Uj∂xi)+23δijk]∂Ui∂xj,πk=-12νs∂∂xk[kε(ε^ε)2∂ε^∂xk],Ekd=∂∂xk(νaλuicik¯)≈0,Ekε=-νaλ(cik∂ui∂xk¯)=-νa2λΛii,
and the transport equation of the dissipation of turbulent kinetic energy is written as
(24)Dε~Dt=Dε+(Pε-Y)+E+Eε=∂∂xj[(νs+νtσε)∂ε∂xj]+ε~k(Cε1Pk-Cε2ε~-ε^)+E+Eε,E=Cε3νsνt(∂2U∂xk2)2,Eε=2νsνaλ∂ui∂xk∂∂xk(∂cij∂xj)¯≈0,
where Ekd and Eε 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 Ekε by applying (22). Following Kawashima and Kawamura [28], the model constants in the above equations are taken as Cε1=1.44, Cε2=1.92, Cε3=0.6, and
(25)σk=σε=1.0-0.5exp[-(yη20)2].Prior to discussion on the feasibility of the present model for the drag-reducing viscoelastic fluid, Figure5 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 atReτ=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 beReτ=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 U=k=ε~=0 and the zero gradient boundary condition is applied to all variables at the channel center (except V=0). In this paper, the superscript (+) describes a nondimensionalized quantity by the friction velocity uτ, the effective viscosity [1], and the density of solution ρ, while (*) describes a nondimensionalized quantity by the channel half-width δ, uτ, and ρ. The friction Reynolds and Weissenberg numbers to be used in the following are defined as Reτ=ρuτδ/μ0 and Wiτ=ρuτ2λ/μ0, 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 Table1. Here, the percent drag reduction is determined as
(26)DR%=(1-CfviscCfNewt)×100(%),
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%.
Reynolds numberReτ
150
150
150
150
150
395
395
Weissenberg numberWiτ
10
11
30
30
40
11
30
Shear-viscosity ratioβ
0.8
0.5
0.5
0.3
0.5
0.5
0.5
Present model
15%
1%
57%
63%
68%
4%
60%
DNS [18]
13%
23%
58%
64%
67%
27%
58%


### 5.1. Mean Flow Statistics

Figure6 shows the mean velocity profiles, as a function of y+, computed using the present model together with DNS data. Note again that, for the drag-reduced flows, y+ 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 y 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μ0.
Reynolds numberReτ
150
150
150
150
150
395
395
Weissenberg numberWiτ
10
11
30
30
40
11
30
Shear-viscosity ratioβ
0.8
0.5
0.5
0.3
0.5
0.5
0.5
Present model
0.961
0.880
0.697
0.527
0.654
0.880
0.697
DNS [18]
0.962
0.870
0.690
0.518
0.654
0.870
0.690Mean velocity profiles in wall coordinates. Comparison between the model prediction and the corresponding DNS data [18].
(a)
Re
τ = 150
(b)
Re
τ = 395From 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 lowWiτ 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).Figures7 and 8 compare the predictions of Cf-Rem and DR%-Rem 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 Cf 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 Cf 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 -uv¯ still remains to be nonzero value. This result argues against the experimental fact that the contribution of -uv¯ 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 Figures9 and 10, the predicted k+ and C11 follow those trends in the DNS data, although the near-wall peak value reveals clearly k+ 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 k-ε model, because highly drag-reduced viscoelastic flows should be inevitably accompanied by enhanced anisotropy in the near-wall region.Same with Figure6, but with respect to turbulent kinetic energy, k+.
(a)
Re
τ = 150
(b)
Re
τ = 395Same with Figure6, but with respect to a conformation normal stress, C11+.
(a)
Re
τ = 150
(b)
Re
τ = 395


### 5.2. Reynolds and Conformation Stresses

The total shear stress can be obtained as(27)τtotal+=βReτ∂U+∂y*-u+v+¯+1-βWiτC12.
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:
(28)τtotal+=1-y*,
proving the validity of the model calculation.Shear stress profiles of drag-reducing viscoelastic flow in channel.
(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)The viscous shear stress is well represented by the present model throughout the entire channel, as shown in Figure11. 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  -u+v+¯ 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<y+<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 ofCij computed by the present model is shown in Figure 12. As noted above with respect to C12, 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 forReτ=395, Wiτ=30, and β=0.5.Figure13 shows the budget in the k 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 y+≈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 y+=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 k-ε model to be applied to the strongly anisotropic turbulence, as similar to that observed already in the overprediction of k. The viscoelastic contribution Ekε 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 k) due to too small Ekε.Figure 13
Budget of turbulent kinetic energy for drag-reducing viscoelastic fluid in turbulent channel flow atReτ=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 k, an increase in ε~, and/or a decrease of fv. Since, as shown already, both variations of k 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 fv. Figure 14 shows the damping function, the product of fμ and fv, 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, y/δ>0.1-0.2, and probably attributed to the inaccuracy (overestimation) of the turbulent kinetic energy in the highly turbulent region.Validation of damping function, the product offμ and fv.
(a)
Re
τ = 150
(b)
Re
τ = 395


## 5.1. Mean Flow Statistics

Figure6 shows the mean velocity profiles, as a function of y+, computed using the present model together with DNS data. Note again that, for the drag-reduced flows, y+ 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 y 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μ0.
Reynolds numberReτ
150
150
150
150
150
395
395
Weissenberg numberWiτ
10
11
30
30
40
11
30
Shear-viscosity ratioβ
0.8
0.5
0.5
0.3
0.5
0.5
0.5
Present model
0.961
0.880
0.697
0.527
0.654
0.880
0.697
DNS [18]
0.962
0.870
0.690
0.518
0.654
0.870
0.690Mean velocity profiles in wall coordinates. Comparison between the model prediction and the corresponding DNS data [18].
(a)
Re
τ = 150
(b)
Re
τ = 395From 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 lowWiτ 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).Figures7 and 8 compare the predictions of Cf-Rem and DR%-Rem 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 Cf 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 Cf 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 -uv¯ still remains to be nonzero value. This result argues against the experimental fact that the contribution of -uv¯ 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 Figures9 and 10, the predicted k+ and C11 follow those trends in the DNS data, although the near-wall peak value reveals clearly k+ 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 k-ε model, because highly drag-reduced viscoelastic flows should be inevitably accompanied by enhanced anisotropy in the near-wall region.Same with Figure6, but with respect to turbulent kinetic energy, k+.
(a)
Re
τ = 150
(b)
Re
τ = 395Same with Figure6, but with respect to a conformation normal stress, C11+.
(a)
Re
τ = 150
(b)
Re
τ = 395


## 5.2. Reynolds and Conformation Stresses

The total shear stress can be obtained as(27)τtotal+=βReτ∂U+∂y*-u+v+¯+1-βWiτC12.
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:
(28)τtotal+=1-y*,
proving the validity of the model calculation.Shear stress profiles of drag-reducing viscoelastic flow in channel.
(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)The viscous shear stress is well represented by the present model throughout the entire channel, as shown in Figure11. 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  -u+v+¯ 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<y+<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 ofCij computed by the present model is shown in Figure 12. As noted above with respect to C12, 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 forReτ=395, Wiτ=30, and β=0.5.Figure13 shows the budget in the k 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 y+≈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 y+=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 k-ε model to be applied to the strongly anisotropic turbulence, as similar to that observed already in the overprediction of k. The viscoelastic contribution Ekε 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 k) due to too small Ekε.Figure 13
Budget of turbulent kinetic energy for drag-reducing viscoelastic fluid in turbulent channel flow atReτ=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 k, an increase in ε~, and/or a decrease of fv. Since, as shown already, both variations of k 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 fv. Figure 14 shows the damping function, the product of fμ and fv, 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, y/δ>0.1-0.2, and probably attributed to the inaccuracy (overestimation) of the turbulent kinetic energy in the highly turbulent region.Validation of damping function, the product offμ and fv.
(a)
Re
τ = 150
(b)
Re
τ = 395


## 6. Conclusion

A new low-Reynolds-numberk-ε 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 thek transport equation. It seems that the deficiency is associated with also the inadequacy in the isotropic k-ε 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 k-ε 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.



---

*Source: 197628-2013-11-27.xml*