Research Article  Open Access
Takahiro Tsukahara, Yasuo Kawaguchi, "Proposal of Damping Function for LowReynoldsNumber  Model Applicable in Prediction of Turbulent ViscoelasticFluid Flow", Journal of Applied Mathematics, vol. 2013, Article ID 197628, 15 pages, 2013. https://doi.org/10.1155/2013/197628
Proposal of Damping Function for LowReynoldsNumber  Model Applicable in Prediction of Turbulent ViscoelasticFluid Flow
Abstract
A lowReynoldsnumber k model applicable for viscoelastic fluid was proposed to predict the frictionaldrag reduction and the turbulence modification in a wallbounded 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 dragreduced turbulent channel flows demonstrated good agreement with DNS (direct numerical simulation) data.
1. Introduction
It is known that minute amount of additives, such as polymer and surfactant, can reduce the frictional drag in 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 19971998 [1, 2]: see reviews [3, 4]. They used several kinds of viscoelastic models for a rheological constitutive equation, such as the FENEP (finitely extensible nonlinear elasticPeterlin) model, the OldroydB 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 turbulenteddy (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, nearwall 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 outerlayer 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 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 Reynoldsaveraged NavierStokes (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 physicallybased turbulence closures for the prediction of flows with dragreducing additives with the aid of the DNS database and the alreadyproposed 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 dragreducing viscoelastic fluids described by the FENEP 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 Reynoldsstress transport equation. Pinho and coworkers [12–14] have developed a  model for the FENEP fluid, using the Reynolds averaged form of the conformation tensor equation to determine the average polymer stress. They devoted their efforts to devise the closures for the new terms appearing in the governing equations, especially, the nonlinear turbulence distortion term (represented by in this paper) in the averaged equation of the conformation tensor. Their closure models were calibrated well on the basis of DNS database, Iaccarino et al. [15] employed the  approach to treat the intensed wall damping of wallnormal turbulence in the dragreduced flow. Moreover, they modelled directly the Reynoldsaveraged 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 smallscale vortical motions near the wall. Hence, the knowledge of scales of eddies that are suppressed by additives is indispensable to discuss the form of the present model. We have obtained this information from existing DNS database [18] enabling a priori development of turbulence closures. To analyze the data, we performed an investigation of frequency response of the Giesekus model and considered the relationship between the shift of the dissipation scale and the viscoelastic behavior of the fluid.
The following section provides a brief introduction of the currently adopted model, which is originally for the Newtonian flow. Thereafter, the scale shifts of the dissipation range and of the viscoelastic behavior are considered and the present constitutive equation as well as the and equations are then developed. Finally, the proposed model will be tested in comparison with DNS results.
2. Introduction of LowRe  Model
The present section commences with a brief review of established models for turbulent flows of Newtonian fluid using twoequation turbulent viscosity closures. In the following, the vector denotes the coordinate directions; is the velocity fluctuating component, while those capital letters represent the mean components. The overbar indicates an ensemble average over all homogeneous directions and time.
The low Reynolds number  model for a turbulent flow relates the Reynolds stress to the mean strain field as follows: 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 as follows: On the other hand, some researchers introduced the wallnormal distance normalized in wall unit, , 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 nearwall region.
As for a viscoelasticfluid flow accompanied by DR, both timescale and lengthscale of turbulenteddy 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 dragreducing 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 turbulenteddy 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 Troutonratio 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 viscoelasticfluid 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 dragreducing 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 eddyscale range that can be shifted in the dragreduced 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 reexamining 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 . The onedimensional 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 wallnormal height according to the mixing length theory [24] and, hence, the energycontaining largescale 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 [25]. Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonianfluid turbulence would occur at where 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 , the lines of (7) and (10) intersect at . 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 dragreducing 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, socalled 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 energydissipation spectrum function of . For a viscoelastic flow at the friction Weissenberg number , the is remarkably decreased at every height, especially in the nearwall region, compared to the case of . As a result, the intersection of the line and the line moves to a higher position of about , as given in Figure 4. The is prone to shift leftwards as the Weissenberg number is increasing consistent to the above Lumley’s dragreduction model.
3.2. Dynamic Characterization of Viscoelastic Fluid
In this section, we investigate the dynamic characterization of the viscoelastic fluid, in order to discuss the effect of its viscoelasticity on turbulent fluctuations (eddies). A model fluid for this study is based on the Giesekus model [27], which is one of the typical viscoelastic models to describe the relationship between stress and strain exerted on the fluid of interest. The rheological constitutive equation with respect to the viscoelastic stress tensor is written as where is the relaxation time, is the zero shear rate viscosity of the solution, and is the mobility factor to indicate the level of nonlinearity.
It is generally known that a measurement under a small amplitude oscillatory flow is a powerful tool for describing the microscopic state of the test fluid. Therefore, we consider here a small fluid volume under a simple shear flow whose shear rate varies periodically with time as 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 is derived explicitly from the viscoelastic stress tensor; that is, . Note that the complex shear modulus consists of the dynamic storage modulus and the dynamic loss modulus as follows: 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 . 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 , for the laminar simple shear flow. Here, the Weissenberg number is defined as the product of and . 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 varies from 11 to 30, the frequency of the crossover point shifts from to . This shift factor corresponds to the ratio of the two Weissenberg numbers; that is, . 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 highfrequency inputs associated with smallscale eddies in the nearwall turbulence. Correspondingly, the relationship between the energydissipation 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.
4. Model Development
While some details of the present model differ from Jones and Launder [19], their lowReynoldsnumber  model will be applied unaltered in most respects. As for the effect of viscoelastic fluid, the above observation has driven us to consider the timescale 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 dragreducing 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 timescale ratio can be proposed as where , , 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 zeroshearrate viscosity of the solution , respectively. Then, replacing (2), let us modify by introducing an additional damping function based on this timescale 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 was chosen ad hoc and adopted to . The presently optimized value of the additional model coefficient is 0.36.
The conventional damping function , that is used for the lowReynoldsnumber 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 wallnormal distance normalized by the Kolmogorov length scale, . The model of by (17) and (19) satisfies the following two assumptions: first, very smallscale eddies are damped by the elastic behavior of the fluid and thereby the eddy viscosity decreases ( 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 ensembleaveraged constitutive equation, which govern the incompressible viscoelasticfluid 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 lefthand side of (21) are small enough to be ignored. The term labeled as in the lefthand side is modeled in the following form: where and are given based on the DNS data [18]. This form was originally proposed by Leighton et al. [11].
From the NavierStokes 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 , , , and
Prior to discussion on the feasibility of the present model for the dragreducing 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 .
5. Results
Results are presented below for dragreduced turbulent channel flows. The tested range of the flow and rheological parameters were set to be and 395, –0.8, , and –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 centraldifference 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 and the zero gradient boundary condition is applied to all variables at the channel center (except ). 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 halfwidth , , and . The friction Reynolds and Weissenberg numbers to be used in the following are defined as and , 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 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 results in quite low , 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 nearzero would produce a negative , as shown later. This issue is not serious, because such a combination of low and low should not be impractical rheological condition; usually, dilute polymer solutions should be equivalent to viscoelastic fluids with low and high (1).

5.1. Mean Flow Statistics
Figure 6 shows the mean velocity profiles, as a function of , computed using the present model together with DNS data. Note again that, for the 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 loglaw 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 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 increases to 40, at which the achieved dragreduction 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 .

(a) = 150
(b) = 395
From the above assessments, it can be said that the Weissenbergnumber 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 and low . The model provided an underestimated velocity for and , as shown in Figure 6(b), while the result of a set of and seems to be predicted well, as given in Figure 6(a).
Figures 7 and 8 compare the predictions of  and  with data from DNS [18]. The computations were made in a range of 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 are demonstrated, enabling us to estimate the maximum that gives rise to the maximum DR. For example, according to the model predictions at , the MDR would be achieved at , 50, and 60 for , 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 than those for MDR as the Weissenberg number increases further (not shown in Figures 7 and 8). With too high , the mean velocity obtained by the model reveals that it exceeds the profile for the MDR asymptote [26] and, in addition, the Reynolds shear stress of still remains to be nonzero value. This result argues against the experimental fact that the contribution of would be almost absent in the highly dragreduced 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 follow those trends in the DNS data, although the nearwall peak value reveals clearly to be overestimated. The discrepancy becomes large as the Weissenberg number increases. As for the wallnormal 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 nearwall region.
(a) = 150
(b) = 395
(a) = 150
(b) = 395
5.2. 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.
(a) ( , , ) = (150, 10, 0.8)
(b) ( , , ) = (150, 30, 0.5)
(c) ( , , ) = (150, 40, 0.5)
(d) ( , , ) = (395, 30, 0.5)
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 are considerably close to the DNS data, indicating that nearly Newtonian fluid calculations would be also executed adequately, as can be seen in Figure 11(a). However, the other three cases shown in the figure have resulted in the overestimated peaks of and their shifts towards the channel center, compared with those by the DNS. As for the viscoelastic stress, the model exhibits tendencies to overestimate it in 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 dragreduction 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 dragreduction rate is being close to the maximum drag reduction. A typical result about the profiles of computed by the present model is shown in Figure 12. As noted above with respect to , the comparison with DNS data shows the present model to underestimate the conformation stresses throughout the channel except for the nearwall region.
Figure 13 shows the budget in the transport equation of the viscoelastic fluid flow at and , 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 is slightly overestimated. Due to both this overestimation and the false negative peak in the pressure diffusion, the turbulent diffusion is enhanced considerably around . 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 .
5.3. Damping Function
It is well known that DR is associated with a decrease in the Reynolds shear stress [4]. By adopting (16), this reduction in the Reynolds shear stress can be achieved through a decrease in , an increase in , and/or a decrease of . Since, as shown already, both variations of and might not provide a significant cause of the reduction, the DR predicted well by the model must be accompanied by a large decrease of the damping function . Figure 14 shows the damping function, the product of and , both for the present model and the estimated value via DNS, indicating that the present model has yielded fairly close agreement with the DNS result. It can be confirmed that the nearwall 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, , and probably attributed to the inaccuracy (overestimation) of the turbulent kinetic energy in the highly turbulent region.
(a) = 150
(b) = 395
6. Conclusion
A new lowReynoldsnumber  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 energydissipative range and the dynamic characterization of the viscoelastic fluid, and was tested by comparison with the DNS data for the dragreduced 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 dragreduction rate were well reproduced by the present model (at least for the present parameters). The present model, however, showed some significant deviations from the DNS data with respect to the Reynolds shear and viscoelastic stresses, when applied for high Weissenberg numbers.
The most serious deficiency found in the present comparisons with DNS data was the overpredicted dissipation rate, which compensated the underestimated term pertaining to the viscoelastic contribution. This error in the overpredictions of the dissipation rate as well as the turbulent kinetic energy would occur in other situations relating to strongly anisotropic turbulence. One cause of this may be the insufficient model for the viscoelastic contribution in the transport equation. It seems that the deficiency is associated with also the inadequacy in the isotropic  model for highly dragreducing viscoelastic flows, although the dragreduction phenomena in terms of bulk flows were predicted well by the present model. Hence, the model predictions should be improved further also in the framework of other models including the anisotropic/nonlinear  model and the Reynolds stress model. In addition, the reproducibility of the MDR state and the application of the present model to practical flows and situations will be investigated in the near future.
In this paper, our discussion is restricted to the case of . 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.01. We should consider the key effect depending on further to improve the model applicable to dilutetodense solutions.
Appendix
For fullydeveloped channel flows, some terms including the spatial derivatives in the horizontal directions can be neglected. In the following, the reduced equations for obtaining timeaveraged values in the channel flow are given with being nondimensionalized by the channel half width , the friction velocity , and the solution zeroshearrate kinematic viscosity .
The momentum equation
the constitutive equation where the transport equation of , normalized by , and the transport equation of , normalized by , The Reynolds stress is given by (1) and (16).
Nomenclature
Skin friction coefficient,  
Mean conformation tensor  
Fluctuation of conformation tensor  
,  Model coefficients 
DR%  Percentage of drag reduction 
,  Damping functions 
Complex shear modulus  
Dynamic storage modulus  
Dynamic loss modulus  
Turbulent kinetic energy,  
Pressure  
Bulk Reynolds number,  
Turbulence Reynolds number,  
Friction Reynolds number,  
Time  
Kolmogorov time scale  
Mean velocity vector  
Fluctuation of velocity vector  
Friction velocity,  
Weissenberg number,  
Friction Weissenberg number,  
Streamwise coordinate  
Wallnormal coordinate  
Spanwise coordinate  
Anisotropic mobility factor  
Ratio of shear viscosities,  
Shear rate  
Channel half width  
Dissipation rate of  
Kolmogorov length scale  
Wave number or frequency  
Wave number or frequency  
Relaxation time  
Solution zeroshearrate shear viscosity,  
,  Shear (kinematic) viscosity of additive contribution 
,  Shear (kinematic) viscosity of solvent contribution 
,  Turbulent (kinematic) viscosity 
Density  
,  Model coefficients 
Model functions  
Viscoelastic stress tensor  
Mean wall shear stress  
Normalized by and  
Normalized by and the effective viscosity  
Statistically averaged. 
Acknowledgments
The authors are grateful to Professor H. Kawamura (Tokyo University of Science, Suwa) for fruitful discussions and comments, and to Mr. H. Kawamoto (who was a master’s course student at Tokyo University of Science) for his work. The DNS database we utilized here to assess results predicted by the present model was acquired using supercomputing resources at the Cyberscience Center at Tohoku University. This work was partially supported by GrantsinAid for Scientific Research (C) no. 25420131 from Japan Society for the Promotion of Science (JSPS).
References
 R. Sureshkumar, A. N. Beris, and R. A. Handler, “Direct numerical simulation of the turbulent channel flow of a polymer solution,” Physics of Fluids, vol. 9, no. 3, pp. 743–755, 1997. View at: Google Scholar
 C. D. Dimitropoulos, R. Sureshkumar, and A. N. Beris, “Direct numerical simulation of viscoelastic turbulent channel flow exhibiting drag reduction: effect of the variation of rheological parameters,” Journal of NonNewtonian Fluid Mechanics, vol. 79, no. 23, pp. 433–468, 1998. View at: Publisher Site  Google Scholar
 C. M. White and M. G. Mungal, “Mechanics and prediction of turbulent drag reduction with polymer additives,” Annual Review of Fluid Mechanics, vol. 40, pp. 235–256, 2008. View at: Google Scholar
 F. C. Li, B. Yu, J. J. Wei, and Y. Kawaguchi, Turbulent Drag Reduction by Surfactant Additives, Wiley, 2012.
 K. Kim and R. Sureshkumar, “Spatiotemporal evolution of hairpin eddies, Reynolds stress, and polymer torque in polymer dragreduced turbulent channel flows,” Physical Review E, vol. 87, Article ID 063002, 11 pages, 2013. View at: Google Scholar
 J. Lumley, “Drag reduction by additives,” Annual Review of Fluid Mechanics, vol. 1, pp. 367–384, 1969. View at: Google Scholar
 L. Thais, T. B. Gatski, and G. Mompean, “Some dynamical features of the turbulent flow of a viscoelastic fluid for reduced drag,” Journal of Turbulence, vol. 13, no. 12, 26 pages, 2012. View at: Publisher Site  Google Scholar  MathSciNet
 L. Thais, G. Mompean, and T. B. Gatski, “Spectral analysis of turbulent viscoelastic and Newtonian channel flows,” Journal of NonNewtonian Fluid Mechanics, vol. 200, pp. 165–176, 2003. View at: Google Scholar
 T. Tsukahara, T. Kawase, and Y. Kawaguchi, “DNS of viscoelastic turbulent channel flow with rectangular orifice at low Reynolds number,” International Journal of Heat and Fluid Flow, vol. 32, no. 3, pp. 529–538, 2011. View at: Publisher Site  Google Scholar
 T. Tsukahara, M. Motozawa, D. Tsurumi, and Y. Kawaguchi, “PIV and DNS analyses of viscoelastic turbulent flows behind a rectangular orifice,” International Journal of Heat and Fluid Flow, vol. 41, pp. 66–79, 2013. View at: Publisher Site  Google Scholar
 R. Leighton, D. T. Walker, T. Stephens, and G. Garwood, “Reynolds stress modeling for drag reducing viscoelastic flows,” in Proceedings of the 4th ASME/JSME Joint Fluids Engineering Conference, pp. 735–744, Honolulu, Hawaii, USA, July 2003. View at: Google Scholar
 D. O. A. Cruz and F. T. Pinho, “Turbulent pipe flow predictions with a low Reynolds number kε model for drag reducing fluids,” Journal of NonNewtonian Fluid Mechanics, vol. 114, no. 23, pp. 109–148, 2003. View at: Publisher Site  Google Scholar
 F. T. Pinho, C. F. Li, B. A. Younis, and R. Sureshkumar, “A low Reynolds number turbulence closure for viscoelastic fluids,” Journal of NonNewtonian Fluid Mechanics, vol. 154, no. 23, pp. 89–108, 2008. View at: Publisher Site  Google Scholar
 P. R. Resende, K. Kim, B. A. Younis, R. Sureshkumar, and F. T. Pinho, “A FENEP kε turbulence model for low and intermediate regimes of polymerinduced drag reduction,” Journal of NonNewtonian Fluid Mechanics, vol. 166, no. 1213, pp. 639–660, 2011. View at: Publisher Site  Google Scholar
 G. Iaccarino, E. S. G. Shaqfeh, and Y. Dubief, “Reynoldsaveraged modeling of polymer drag reduction in turbulent flows,” Journal of NonNewtonian Fluid Mechanics, vol. 165, no. 78, pp. 376–384, 2010. View at: Publisher Site  Google Scholar
 P. R. Resende, F. T. Pinho, B. A. Younis, K. Kim, and R. Sureshkumar, “Development of a LowReynoldsnumber kω Model for FENEP Fluids,” Flow Turbulence and Combustion, vol. 90, pp. 89–108, 2013. View at: Google Scholar
 J. Bredberg, S.H. Peng, and L. Davidson, “An improved kω turbulence model applied to recirculating flows,” International Journal of Heat and Fluid Flow, vol. 23, no. 6, pp. 731–743, 2002. View at: Publisher Site  Google Scholar
 T. Tsukahara, T. Ishigami, B. Yu, and Y. Kawaguchi, “DNS study on viscoelastic effect in dragreduced turbulent channel flow,” Journal of Turbulence, vol. 12, article N13, pp. 1–25, 2011. View at: Publisher Site  Google Scholar
 W. P. Jones and B. E. Launder, “The prediction of laminarization with a twoequation model of turbulence,” International Journal of Heat and Mass Transfer, vol. 15, no. 2, pp. 301–314, 1972. View at: Google Scholar
 Y. Nagano and M. Tagawa, “An improved k$\epsilon $ model for boundary layer flows,” Journal of Fluids Engineering, vol. 112, no. 1, pp. 33–39, 1990. View at: Google Scholar
 E. R. van Driest, “On turbulent flow near a wall,” Journal of the Aeronautical Sciences, vol. 23, no. 11, pp. 1007–1011, 1956. View at: Google Scholar
 H. Abe, H. Kawamura, and Y. Matsuo, “Direct numerical simulation of a fully developed turbulent channel flow with respect to the Reynolds number dependence,” Journal of Fluids Engineering, vol. 123, pp. 382–393, 2001. View at: Google Scholar
 K. R. Sreenivasan and C. M. White, “The onset of drag reduction by dilute polymer additives, and the maximum drag reduction asymptote,” Journal of Fluid Mechanics, vol. 409, pp. 149–164, 2000. View at: Google Scholar
 H. Tennekes and J. L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, Mass, USA, 1972.
 S. B. Pope, Turbulent Flows, Cambridge University Press, Cambridge, Mass, USA, 2000. View at: MathSciNet
 P. S. Virk, “An elastic sublayer model for drag reduction by dilute solutions of linear macromolecules,” Journal of Fluid Mechanics, vol. 45, pp. 417–440, 1971. View at: Google Scholar
 H. Giesekus, “A simple constitutive equation for polymer fluids based on the concept of deformationdependent tensorial mobility,” Journal of NonNewtonian Fluid Mechanics, vol. 11, no. 12, pp. 69–109, 1982. View at: Google Scholar
 N. Kawashima and H. Kawamura, “Reconstruction of kε model with relevance to the wall asymptotic behavior,” Transactions of the Japan Society of Mechanical Engineers B, vol. 62, no. 597, pp. 1700–1708, 1996 (Japanese). View at: Google Scholar
 M. D. Warholic, D. K. Heist, M. Katcher, and T. J. Hanratty, “A study with particleimage velocimetry of the influence of dragreducing polymers on the structure of turbulence,” Experiments in Fluids, vol. 31, no. 5, pp. 474–483, 2001. View at: Publisher Site  Google Scholar
 M. Motozawa, T. Watanabe, W. Gu, and Y. Kawaguchi, “Experimental investigation on zonal structure in dragreducing channel flow with surfactant additives,” Advances in Mechanical Engineering, vol. 2011, Article ID 120438, 12 pages, 2011. View at: Publisher Site  Google Scholar
 B. Yu and Y. Kawaguchi, “Direct numerical simulation of viscoelastic dragreducing flow: a faithful finite difference method,” Journal of NonNewtonian Fluid Mechanics, vol. 116, no. 23, pp. 431–466, 2004. View at: Publisher Site  Google Scholar
 B. Yu, F. Li, and Y. Kawaguchi, “Numerical and experimental investigation of turbulent characteristics in a dragreducing flow with surfactant additives,” International Journal of Heat and Fluid Flow, vol. 25, no. 6, pp. 961–974, 2004. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2013 Takahiro Tsukahara and Yasuo Kawaguchi. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.