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

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 fluidmodel. Assessment of the proposedmodel’s performance in several rheological conditions for drag-reduced turbulent channel flows demonstrated good agreement with DNS (direct numerical simulation) data.


Introduction
It is known that minute amount of additives, such as polymer and surfactant, can reduce the frictional drag in wallbounded 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 energysaving 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 pressurestrain (redistribution) term in the Reynolds-stress transport equation.Pinho and coworkers [12][13][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 V 2 - 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 Newtonianfluid 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.

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   = ( 1 ,  2 ,  3 ) ≡ (, , ) denotes the coordinate directions;   = ( 1 ,  2 ,  3 ) ≡ (, V, ) 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: Jones and Launder [19] Nagano and Tagawa [20] DNS by Abe et al. [22] Viscoelastic fluid DNS by Tsukahara et al. [18] Wi  = 11,  = 0.5 Wi  = 30,  = 0.5 where   is the Kronecker delta and the turbulent kinematic viscosity is expressed as 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: ) .
On the other hand, some researchers introduced the wallnormal distance normalized in wall unit,  + =   / 0 , as a parameter of   .A typical example is the one proposed by Nagano and Tagawa [20] as follows: Their damping function is based on the van Driest walldamping 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.
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 dragreducing 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 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.

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 dragreduced 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 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 which then represents the limit of the production range of the spectrum.On the other hand, the scale of dissipative smallscale eddies can be estimated from the energy dissipation spectrum, described as  2    () [25].Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonianfluid turbulence would occur at 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 using the von Kármán constant   and the friction velocity   [24].Therefore, the relation between  and   is rewritten from (8) as 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.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  2    (  ).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.

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 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 where  is the frequency of the oscillatory shear stress.The dimensionless complex modulus as one of rheological material functions is obtained from where the nondimensional conformation tensor  12 is derived explicitly from the viscoelastic stress tensor; that is,    =  0 (  −   )/.Note that the complex shear modulus  consists of the dynamic storage modulus   and the dynamic loss modulus   as follows: 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  + 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 Figure 4. From this result, it can be said that viscoelastic fluid should behave more elastically against highfrequency eddy motions in turbulence and suppress them.

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 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  0 , respectively.Then, replacing (2), let us modify ]  by introducing an additional damping function  V based on this time-scale ratio, as follows: Here, the "homogeneous" dissipation rate ε (=  − ε) is introduced, while is defined with an assumption that  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  V .The presently optimized value of the additional model coefficient  0 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: where   is the wall-normal distance normalized by the Kolmogorov length scale,  = (] 3  /) 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: 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: where  1 = 0.001(1 − exp(−  /100)) and  2 = −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 and the transport equation of the dissipation of turbulent kinetic energy is written as 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 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.

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 secondorder central-difference scheme for discretization, while convergence calculation was performed by the successive overrelaxation (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  =   / 0 and Wi  =  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 Table 1.Here, the percent drag reduction is determined as where the suffixes "visc" and "Newt" stand for the frictioncoefficient 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).

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 dragreduced 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%.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 −V still remains to be nonzero value.This result argues against the experimental fact that the contribution of −V 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.
As can be seen in Figures 9 and 10, the predicted  + and  11 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 dragreduced viscoelastic flows should be inevitably accompanied by enhanced anisotropy in the near-wall region.

Reynolds and Conformation
Stresses.The total shear stress can be obtained as Here, the terms are (in order from left to right of the righthand side) the viscous shear stress, the Reynolds shear stress, and the viscoelastic shear stress.Figure 11 shows the wallnormal 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: proving the validity of the model calculation.
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 − + 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 <  + < 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  12 , the comparison with DNS data shows the present model to underestimate the conformation stresses throughout the channel except for the near-wall region.
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    .

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  V .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  V .Figure 14 shows the damping function, the product of   and  V , 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.

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 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 zeroshear-rate kinematic viscosity ].
The momentum equation

Figure 1 :
Figure 1: Comparison of damping function   for turbulent channel flow at Re  = 395.

Figure 4 :
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.

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

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

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

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

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