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

1. Introduction

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

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

As for high Reynolds number simulations, the available information on the viscoelastic turbulent flow is even more scarce than for the Newtonian turbulence, despite of the development of computer resources in the recent decades. A noticeable exception is the DNS performed by Thais et al. [7, 8], where high DR flows at the friction Reynolds number 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 [1214] have developed a - model for the FENE-P fluid, using the Reynolds averaged form of the conformation tensor equation to determine the average polymer stress. They devoted their efforts to devise the closures for the new terms appearing in the governing equations, especially, the nonlinear turbulence distortion term (represented by in this paper) in the averaged equation of the conformation tensor. Their closure models were calibrated well on the basis of DNS database, Iaccarino et al. [15] employed the - approach to treat the intensed wall damping of wall-normal turbulence in the drag-reduced flow. Moreover, they modelled directly the Reynolds-averaged polymer stress without solving the constitutive equation for the conformation tensor, resulting in less coefficients and functions. The results predicted by these models were generally in good agreement with DNS database for low and intermediate DR. More recently, Resende et al. [16] proposed a - model based on the model of Bredberg et al. [17] for Newtonian fluids. Some improvement was achieved for predictions in the intermediate DR region, but still not enough to estimate and the high DR (in which the drag reduction rate is more than 60%).

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

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

2. Introduction of Low-Re - Model

The present section commences with a brief review of established models for turbulent flows of Newtonian fluid using two-equation turbulent viscosity closures. In the following, the vector denotes the coordinate directions; is the velocity fluctuating component, while those capital letters represent the mean components. The overbar indicates an ensemble average over all homogeneous directions and time.

The low Reynolds number - model for a turbulent flow relates the Reynolds stress to the mean strain field as follows: 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 wall-normal distance normalized in wall unit, , as a parameter of . A typical example is the one proposed by Nagano and Tagawa [20] as follows: Their damping function is based on the van Driest wall-damping formula [21]. They considered the behavior of turbulence on an infinite flat plate with a simple harmonic oscillation parallel to the plate. The amplitude of the motion of the fluid in each position was examined to describe the damping function. Figure 1 shows the distribution of these individual models when applied to the turbulent channel flow at a given Reynolds number. The model proposed by Nagano and Tagawa [20] is in better agreement with the DNS data for the Newtonian flow [22], especially, in the near-wall region.

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

3. Shift of Dissipative Scale

3.1. Lumley’s Consideration of Time Criteria

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

It is assumed that the mean pressure gradient is constant and the turbulent flow field is fully developed, where the channel width is . 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 small-scale eddies can be estimated from the energy dissipation spectrum, described as [25]. Lumley [6] assumed that the peak of the dissipation spectrum for the Newtonian-fluid turbulence would occur at 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 drag-reducing viscoelastic flow, the relationship of and is expected to be changed as shown with broken lines in Figure 2, according to Lumley [6]. He reasoned that less difference between and as a result of decreased might lead to the expanding of the viscous sublayer. The intersection of and , therefore, shifts from A to a higher position labeled as B. Assuming that B is the shifted bound of the viscous sublayer in the relevant flow, the viscous sublayer of the viscoelastic fluid should be enlarged. In other words, the new region, so-called the elastic layer [26], arises in the vicinity of the wall with the upper bound at the height of B.

In order to verify this concept of DR, we reexamine the existing DNS data of the turbulent channel flow accompanied by DR [18], as shown in Figure 3. The value of is determined from the peak of the spanwise energy-dissipation spectrum function of . For a viscoelastic flow at the friction Weissenberg number , the is remarkably decreased at every height, especially in the near-wall 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 drag-reduction 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 high-frequency inputs associated with small-scale eddies in the near-wall turbulence. Correspondingly, the relationship between the energy-dissipation range and the frequency response of the Giesekus model reveals that the shift in wave number of the energy dissipation is consistent with the change of viscoelastic behavior, as shown in Figure 4. From this result, it can be said that viscoelastic fluid should behave more elastically against high-frequency eddy motions in turbulence and suppress them.

4. Model Development

While some details of the present model differ from Jones and Launder [19], their low-Reynolds-number - model will be applied unaltered in most respects. As for the effect of viscoelastic fluid, the above observation has driven us to consider the time-scale ratio between the characteristic time scale of turbulent eddy motions and the relaxation time of the relevant fluid, and we employ it to describe the diminished values of . In the drag-reducing viscoelastic fluid that has a relaxation time comparable to some scale of turbulent eddies, the smaller eddies, if larger than the Kolmogorov scale, should be damped due to the elastic behavior of the fluid. Using the Kolmogorov time as the characteristic time scale in turbulence, the time-scale ratio can be proposed as 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 zero-shear-rate viscosity of the solution , respectively. Then, replacing (2), let us modify by introducing an additional damping function based on this time-scale ratio, as follows: 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 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, . 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 ( 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 and 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 , , , 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 .

5. Results

Results are presented below for drag-reduced turbulent channel flows. The tested range of the flow and rheological parameters were set to be 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 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 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 half-width , , 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 friction-coefficient 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 near-zero 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 drag-reduced flows, is normalized by the effective viscosity of each case, whose values are described in Table 2. Also shown in Figure 6 is the MDR asymptote, where the log-law profile for maximum drag reduction is that proposed by Virk [26]. At each Weissenberg number, the present model was in good agreement with the DNS data. Especially, the model applied for 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 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 .

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 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 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 follow those trends in the DNS data, although the near-wall peak value reveals clearly to be overestimated. The discrepancy becomes large as the Weissenberg number increases. As for the wall-normal height of the peak, the model shows almost good agreements with the DNS data. It would be difficult to improve this aspect in the framework of the isotropic - model, because highly drag-reduced viscoelastic flows should be inevitably accompanied by enhanced anisotropy in the near-wall region.

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 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: 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 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 drag-reduction rate on the Weissenberg number, as given in Table 1.

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

Figure 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 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, -, and probably attributed to the inaccuracy (overestimation) of the turbulent kinetic energy in the highly turbulent region.

6. Conclusion

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

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

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

In this paper, our discussion is restricted to the case of . 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 dilute-to-dense solutions.


For fully-developed channel flows, some terms including the spatial derivatives in the horizontal directions can be neglected. In the following, the reduced equations for obtaining time-averaged values in the channel flow are given with being nondimensionalized by the channel half width , the friction velocity , and the solution zero-shear-rate kinematic viscosity .

The momentum equation

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).


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,
Bulk Reynolds number,
Turbulence Reynolds number,
Friction Reynolds number,
Kolmogorov time scale
Mean velocity vector
Fluctuation of velocity vector
Friction velocity,
Weissenberg number,
Friction Weissenberg number,
Streamwise coordinate
Wall-normal 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 zero-shear-rate shear viscosity,
, Shear (kinematic) viscosity of additive contribution
, Shear (kinematic) viscosity of solvent contribution
, Turbulent (kinematic) viscosity
, Model coefficients
Model functions
Viscoelastic stress tensor
Mean wall shear stress
Normalized by and
Normalized by and the effective viscosity
Statistically averaged.


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