A FENE-P k–ε Viscoelastic Turbulence Model Valid up to High Drag Reduction without Friction Velocity Dependence

Abstract

A viscoelastic turbulence model in a fully-developed drag reducing channel flow is improved, with turbulent eddies modelled under a $k–\epsilon$ representation, along with polymeric solutions described by the finitely extensible nonlinear elastic-Peterlin (FENE-P) constitutive model. The model performance is evaluated against a wide variety of direct numerical simulation data, described by different combinations of rheological parameters, which is able to predict all drag reduction (low, intermediate and high) regimes with good accuracy. Three main contributions are proposed: one with a simplified viscoelastic closure for the ${NLT}_{ij}$ term (which accounts for the interactions between the fluctuating components of the conformation tensor and the velocity gradient tensor), by removing additional damping functions and reducing complexity compared with previous models; second through a reformulation for the closure of the viscoelastic destruction term, $E_{{\tau }_p}$, which removes all friction velocity dependence; lastly by an improved modified damping function capable of predicting the reduction in the eddy viscosity and thus accurately capturing the turbulent kinetic energy throughout the channel. The main advantage is the capacity to predict all flow fields for low, intermediate and high friction Reynolds numbers, up to high drag reduction without friction velocity dependence.

1. Introduction

Since the pioneering experiment by Toms [1], it is known that the additions of small (parts per million) amounts of long-chain flexible polymers to a turbulent flow can drastically reduce the transport energy by decreasing the turbulent drag. The effects are most evident in turbulent shear flow, in which dissolving the polymers in solution can reduce friction losses by as much as 80% compared to the solvent alone [2]. After the discovery of the drag reduction (DR) phenomena, several comprehensive studies were carried out to understand the physical mechanisms of the interactions between the turbulent structures and polymer chains. Early comprehensive studies in this area come from Lumley [3,4], Hoyt [5] and Virk [6]. Lumley suggests that the DR phenomenon is the result of an increase in effective viscosity in an area outside the viscous sub-layer and in the buffer layer, caused by polymer chains stretching in a turbulent flow.

More recent studies have been proposed for the theory of the mechanisms of drag-reducing polymer additives [7,8]. Several Direct Numerical Simulation (DNS) studies were conducted to understand further the energy exchanges between the polymer chains and turbulence structures [9,10,11,12,13,14]. It is now known to at least low to moderate levels of DR, that the mechanism is the suppression of the near-wall streamwise vortices by polymers that stretch in the extensional flow, and then relax as they are rolled into other vortices, generating forces that tend to weaken these vortices. Quantitatively, this can be expressed as a polymer body force [9], which is positive in the streamwise direction, with an opposite sign (anti-correlation) in the wall direction.

DNS is a great resource to explore the underlying mechanics of drag-reducing viscoelastic turbulent flows. However, for the majority of engineering motivations, DNS is not practical because of the high number of variables which requires a substantial expense of memory and CPU-time. This cost is more prevalent in high DR (HDR) schemes in which the near-wall velocity streaks become more elongated, requiring an increased demand on computational resources.

An alternative approach in capturing flow features at much less computational demands is the application of Reynolds-averaged Navier–Stokes (RANS) models, whose interest has increased in recent years. One of the original implementations of elastic effects within turbulence models was achieved by Pinho [15] and Cruz et al. [16]. Their work focused on low-Reynolds number $k–\epsilon$ turbulence models, applying a Generalised Newtonian Fluid (GNF) constitutive equation involving dependency of the fluid strain hardening on the third invariant of the rate of deformation tensor. Following these studies, an anisotropic version was also developed which included an increased Reynolds stress anisotropy [17], along with a Reynolds stress turbulence model [18], both able to satisfactorily predict drag-reducing behavior. Nevertheless, the models are constrained because of the inelastic formulation of the GNF constitutive equations.

Further developments in viscoelastic RANS models became possible owing to the emergence of DNS data regarding turbulent viscoelastic fluids. The first elastic model was developed by Leighton et al. [19], which was based on the finitely extensible nonlinear elastic with Peterlin closure (FENE-P) dumbbell constitutive equation model. Their study involved the development of a polymer strain–stress coupling based on the tensor expansion, which incorporated the conformation tensor and Reynolds stress. From this work, more attention arose to the FENE-P model given the molecular roots of the equations. Later, based on a-priori analysis of DNS data, Pinho et al. [20] developed a low-Reynolds number $k–\epsilon$ model for FENE-P fluids which could predict flow features up to the low drag reduction regime (LDR < 20%). Turbulent viscoelastic closures were proposed, including the non-linear term involving the conformation tensor and the strain rate fluctuations within the conformation tensor equation (denoted ${NLT}_{ij}$ following the nomenclature of Housiadas et al. [21] and Li et al. [22]); along with the viscoelastic turbulent transport term of the turbulent kinetic energy. One of the key difficulties that arose in this initial study was the decrease in the magnitude of turbulent kinetic energy as viscoelasticity increased, opposite to that found in the DNS literature [23]. Subsequently, the model closures were improved by Resende et al. [24] and the capacity of the model predictions were extended to the intermediate drag reduction regime (20% < IDR < 40%). However, the model closures involved complex damping functions and model constants which gave spurious results for the high drag reduction regime (HDR > 40%). Resende et al. [25] applied the same viscoelastic closures to a low-Reynolds number $k–\omega$ model with only a mathematical transformation of the governing terms involving $\omega$. The closures had identical limitations as the $k–\epsilon$ model for predicting DR behaviour but demonstrated great versatility and robustness given its application to alternative two-equation models.

During this time, a $k–\epsilon –\overline{v^2}–f$ model for FENE-P fluids in fully developed channel flow was proposed by Iaccarino et al. [26], following the initial studies of Dubief et al. [27]. They introduced the idea of a turbulent polymer viscosity which accounts for the effects of viscoelasticity and turbulence on the polymer stress within the momentum equation. The reduction in the Reynolds shear stress is assumed by a-priori DNS data analysis from the decreasing $v^2$ shown within the DNS studies [9]. The model closure for the ${NLT}_{ij}$ is much simpler than the one developed by Resende et al. [24], but contains only the trace and not the individual components. The model was later improved by Masoudian et al. [28] and can predict flow features up to maximum DR (MDR). The key advancement of the closures were an ${NLT}_{ij}$ closure based on DNS analysis and comparisons to the local eddy viscosity peaks; the viscoelastic stress work in the turbulent kinetic energy equations; viscoelastic stress in the momentum equation; and a viscoelastic destruction term in the dissipation transport equation. The viscoelastic turbulent closures within the $v^2$ equation (transverse viscoelastic stress work, ${\epsilon }_{yy}^V$) should be strictly a function of ${NLT}_{yy}$, which is a key component in the formulation of an effective polymer viscosity. However, because only the trace of the ${NLT}_{ij}$ term is present within the model, the closure had to be formulated using DNS analysis of alternative parameters.

Subsequently, after this study, a second-order Reynolds stress model for FENE-P fluids was proposed by Masoudian et al. [29], extending on the idea of a correlation between the Reynolds stresses and the ${NLT}_{ij}$ components, similar to Leighton et al. [19]. The model can predict all DR regimes but is generally unattractive due to the higher number of Newtonian terms resulting from higher-order modeling. Masoudian et al. [30] then further improved the $k–\epsilon –\overline{v^2}–f$ model capabilities via the ${NLT}_{ij}$ term by introducing a simple extension to include heat transfer, along with removing wall dependence via the friction velocity. There are concerning features when one examines the Bousinesq-type ${NLT}_{ij}$ term, which has a zero ${NLT}_{yy}$ component, along with an opposite sign for ${NLT}_{xy}$, both terms being crucial for the polymer shear stress in the momentum balance (see Appendix 1 in Pinho et al. [20]). Further, the increase of k in the buffer layer is small, meaning the decoupling of the $v^2$ component may not be enough to decrease the eddy viscosity. This is compensated by an opposite trend in the dissipation rate, $\epsilon$, for increasing DR, which subsequently balances the momentum equations and causes the necessary increase in the velocity profiles.

An alternative approach in predicting DR flow features other than the use of higher-order models such as the $k–\epsilon –v^2–f$ and Reynolds stress models mentioned previously, is that of a modified damping function or polymer eddy viscosity, accounting for the effect the polymer has on reducing the Reynolds shear stress. Tsukahara and Kawaguchi [31] proposed a modified damping function for a low-Reynolds $k–\epsilon$ model for fluids described by the Giesekus constitutive equation, following the same ideas as Pinho [15] and Cruz. The closure was developed based on the energy-dissipative range and the dynamic characterization of the viscoelastic fluid. The model successfully captures the increase in the magnitude of the turbulent kinetic energy, along with the shift through the buffer layer. Although the magnitudes of k are largely over-predicted in many cases, which is counterbalanced by a lack of closure for the viscoelastic destruction term. In some instances, the model predicts a DR of $1\%$ with a DNS result of $23\%$. Resende et al. [32] proposed a modified damping function for a low-Reynolds number $k–\epsilon$ model for FENE-P fluids which can capture the increase of turbulent kinetic energy as flow viscoelasticity is increased, improving on model predictions made previously by Pinho et al. [20] and Resende et al. [24]. The study also improved largely on the ${NLT}_{ij}$ closure accuracy and simplicity formulated in the previous work. The model is able to predict flow features for a large range of rheological parameters but is limited to a friction Reynolds number of $R{e}_{{\tau }_0}=395$, along with friction velocity still present in the model. For model applicability in flows with reattachment, the friction velocity dependence poses a problem as the values become null at these points and lead to spurious results or floating point errors within computational solvers.

In the present study, an improved $k–\epsilon$ model for FENE-P fluids is proposed, validated for all drag reduction regimes (low, intermediate, high) and up to the largest friction Reynolds number ($R{e}_{{\tau }_0}=1000$) available in the DNS data. The important contribution to the current model is improved and simplified ${NLT}_{ij}$ term that removes complexity from the most recent model developed by Resende et al. [32]; along with a modified damping function which accurately predicts the viscoelastic contributions near and away from the wall, effectively reducing the eddy viscosity and thickening the buffer-layer as DR increases. Further, a reformulation of the viscoelastic destruction term, $E_{{\tau }_p}$, which removes all friction velocity present in the previous $k–\epsilon$ models. The model is assessed against DNS data covering a wide range of flow conditions in terms of the friction Weissenberg number, $W{i}_{{\tau }_0}$, maximum polymer extension, $L^2$, viscosity ratio, $\beta$, and friction Reynolds number, $R{e}_{{\tau }_0}$; along with comparisons against other turbulent FENE-P models within the literature.

The paper is organized as follows: Section 2 introduces the instantaneous and time-averaged governing equations and identifies the viscoelastic terms that will require modeling; Section 3 explains in detail the development of the viscoelastic turbulent closures; Section 4 summarises the model; Section 5 gives the numerical procedure applied; Section 6 presents the results of the flow fields in fully developed channel flow, covering all range of DR and flow conditions; and finally in Section 7, the main conclusions are presented.

2. Governing Equations

The governing equations for incompressible turbulent flow of dilute polymer solutions are the continuity and momentum equations respectively:

$$\frac{\partial {\widehat{u}}_k}{\partial x_k}=0,$$

(1)

$$\rho \frac{\mathrm{D}{\widehat{u}}_i}{\mathrm{D}t}\equiv \rho \left(\frac{\partial {\widehat{u}}_i}{\partial t}+{\widehat{u}}_k\frac{\partial {\widehat{u}}_i}{\partial x_k}\right)=-\frac{\partial \widehat{p}}{\partial x_i}+\frac{\partial {\widehat{\tau }}_{ik}}{\partial x_k},$$

(2)

where the hat represents instantaneous quantities of velocity ${\widehat{u}}_i$, pressure $\widehat{p}$, stress tensor ${\widehat{\tau }}_{ij}$, and fluid density $\rho$. The stress tensor is the sum of the Newtonian solvent which obeys Newton’s law of viscosity, ${\widehat{\tau }}_{ij}^s=2{\mu }_s{\widehat{s}}_{ij}$, with ${\mu }_s$ representing the solvent viscosity coefficient, and polymeric contributions, ${\widehat{\tau }}_{ik}^p$,

$${\widehat{\tau }}_{ij}={\widehat{\tau }}_{ij}^s+{\widehat{\tau }}_{ij}^p.$$

(3)

The kinematic viscosity is used alternatively throughout this study and is defined as $\nu =\mu /\rho$. The instantaneous rate of strain tensor, ${\widehat{s}}_{ij}$, is defined as

$${\widehat{s}}_{ij}=\frac{1}{2}\left(\frac{\partial {\widehat{u}}_i}{\partial x_j}+\frac{\partial {\widehat{u}}_j}{\partial x_i}\right).$$

(4)

The instantaneous polymer contributions are based on the FENE-P rheological dumbbell model [33], with closure given by

$${\widehat{\tau }}_{ij}^p=\frac{{\mu }_p}{\lambda }\left(f\left({\widehat{c}}_{kk}\right){\widehat{c}}_{ij}-{\delta }_{ij}\right),$$

(5)

with

$$f\left({\widehat{c}}_{kk}\right)=\frac{L^2-3}{L^2-{\widehat{c}}_{kk}},$$

(6)

known as the Peterlin function, and ${\widehat{c}}_{kk}$ is the trace of the instantaneous conformation tensor. The other parameters that are associated with the FENE-P model are: $\lambda$, the relaxation time of the polymeric fluid; $L^2$, the maximum extensibility of the dumbbell model; and ${\mu }_p$, the polymer viscosity coefficient.

The behaviour of the instantaneous conformation tensor follows a hyperbolic differential equation of the form,

$$\frac{\partial {\widehat{c}}_{ij}}{\partial t}+{\widehat{u}}_k\frac{\partial {\widehat{c}}_{ij}}{\partial x_k}-\left({\widehat{c}}_{kj}\frac{\partial {\widehat{u}}_i}{\partial x_k}+{\widehat{c}}_{ik}\frac{\partial {\widehat{u}}_j}{\partial x_k}\right)={\stackrel{\nabla }{\widehat{c}}}_{ij}=-\frac{{\widehat{\tau }}_{ij}^p}{{\mu }_p}.$$

(7)

The Oldroyd’s upper convective derivative of the instantaneous conformation tensor is here denoted with ${\stackrel{\nabla }{\widehat{c}}}_{ij}$. The local and advective derivatives are the first and second terms respectively. The bracketed term accounts for the effect of polymer stretching by the instantaneous flow.

The Reynolds averaging process [34] is applied to the governing equations via a Reynolds decomposition of the flow fields such that, ${\widehat{u}}_i=U_i+u_i$; where the use of overbars or upper-case represents the averaged quantity; and primes or lower-case represent the instantaneous quantities. The continuity and momentum equations now take the form:

$$\frac{\partial U_k}{\partial x_k}=0,$$

(8)

$$\rho \frac{\partial U_i}{\partial t}+\rho U_k\frac{\partial U_i}{\partial x_k}=-\frac{\partial \overline{p}}{\partial x_i}+{\mu }_s\frac{{\partial }^2{U}_i}{\partial x_k\partial x_k}-\frac{\partial }{\partial x_k}\left(\rho \overline{u_i{u}_k}\right)+\frac{\partial {\overline{\tau }}_{ik}^p}{\partial x_k},$$

(9)

referred to as the Reynolds-averaged Navier–Stokes (RANS). The Reynolds stress tensor is $\overline{u_i{u}_k}$ and requires a closure model. The Reynolds-averaged polymer stress is ${\overline{\tau }}_{ik}^p$ and written fully as

$${\overline{\tau }}_{ij}^p=\frac{{\mu }_p}{\lambda }\left[f\left(C_{kk}\right)C_{ij}-{\delta }_{ij}\right]+\frac{{\mu }_p}{\lambda }\overline{f(C_{kk}+c_{kk})c_{ij}},$$

(10)

where the additional term on the right requires a closure. The Peterlin function becomes

$$f\left(C_{kk}\right)=\frac{L^2-3}{L^2-C_{kk}}.$$

(11)

After Reynolds averaging, the instantaneous conformation tensor equation becomes

$$\frac{\mathrm{D}{C}_{ij}}{\mathrm{D}t}-M_{ij}+C{T}_{ij}-{NLT}_{ij}=\frac{{\overline{\tau }}_{ij}^p}{{\mu }_p},$$

(12)

$$M_{ij}=C_{jk}\frac{\partial U_i}{\partial x_k}+C_{ik}\frac{\partial U_j}{\partial x_k},$$

(13)

$$C{T}_{ij}=\overline{u_k\frac{\partial c_{ij}}{\partial x_k}},$$

(14)

$${NLT}_{ij}=\overline{c_{jk}\frac{\partial u_i}{\partial x_k}}+\overline{c_{ik}\frac{\partial u_j}{\partial x_k}},$$

(15)

which is referred to as the Reynolds-averaged conformation evolution (RACE). $M_{ij}$ is the mean flow distortion term; it is non-zero, but requires no closure. The remaining two terms are named following the nomenclature of Li et al. [22] and Housiadas et al. [21]. They are labelled with $C{T}_{ij}$; representing the contribution to the transport of the conformation tensor due to the fluctuating advective terms; and ${NLT}_{ij}$, which accounts for the interactions between the fluctuating components of the conformation tensor and the velocity gradient tensor.

Following the analysis of Pinho et al. [20], the nonlinear fluctuating correlation of the average polymeric stress, $\overline{f(C_{kk}+c_{kk})c_{ij}}$ in Equation (10) was shown to be negligible for LDR and HDR when compared with the linear part. This was later neglected in the models of Resende et al. [24] and Masoudian et al. [28] and is also neglected here. The $C{T}_{ij}$ term can also be omitted for all DR regimes following a budget analysis of the RACE carried out by Housiadas et al. [21] and Li et al. [22]. The ${NLT}_{ij}$ term cannot be neglected since it is a significant contributor to the RACE and therefore requires a suitable closure.

2.1. Model for the Reynolds Stress Tensor

The Reynolds stress tensor is computed by adopting the Boussinesq turbulent stress strain relationship,

$$-\rho \overline{u_i{u}_j}=2\rho {\nu }_T{S}_{ij}-\frac{2}{3}\rho k{\delta }_{ij},$$

(16)

where k is the turbulent kinetic energy, $S_{ij}$ is the mean rate of strain tensor and ${\mu }_T=\rho {\nu }_T$ is the eddy viscosity. ${\nu }_T$ is modelled by the typical isotropic $k–\epsilon$ turbulence model for low Reynolds numbers, which includes a damping function $f_{\mu }$ to account for near-wall effects:

$${\nu }_T=C_{\mu }{f}_{\mu }\frac{k^2}{{\tilde{\epsilon }}^N},$$

(17)

where ${\tilde{\epsilon }}^N={\nu }_s\overline{\frac{\partial u_i}{\partial x_j}\frac{\partial u_i}{\partial x_j}}$ is the viscous dissipation of k by the Newtonian solvent,

$$f_{\mu }={\left[1-exp\left(-\frac{y^{+}}{a_{\mu }}\right)\right]}^2,$$

(18)

and $a_{\mu }=26.5$. The dimensionless wall scaling is $y^{+}=u_{{\tau }_0}y/{\nu }_0$, where $u_{{\tau }_0}$ is the friction velocity, y is the distance to the nearest wall, and ${\nu }_0$ is the sum of solvent and polymer viscosity coefficients (${\nu }_0={\nu }_s+{\nu }_p$). The damping function requires additional modelling to capture the anisotropy of the drag reducing flow as a result of viscoelastic flow effects, to be discussed further in this study (Section 3.2).

2.2. Transport Equation for the Turbulent Kinetic Energy

The governing transport equation for the turbulent kinetic energy of turbulent flow with FENE-P fluids is given by,

$$\begin{array}{cc}\hfill \rho \frac{\partial k}{\partial t}+\rho U_i\frac{\partial k}{\partial x_i}=& \rho \frac{\partial }{\partial x_i}\left[\left({\nu }_s+\frac{f_t{\nu }_T}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+P_k-\rho ({\tilde{\epsilon }}^N+D)\hfill \\ & +Q^V-\rho {\epsilon }^V,\hfill \end{array}$$

(19)

with

$$D=2{\nu }_s{\left(\frac{d\sqrt{k}}{d{x}_i}\right)}^2.$$

(20)

$P_k=-\rho \overline{u_i{u}_j}\frac{\partial U_i}{\partial x_j}$ is the rate of production of k.

The Newtonian closures of Equation (19) are those present in the Nagano et al. [35,36] models. To increase numerical stability, a modified Newtonian rate of dissipation of k is applied instead of the true dissipation, which are related by ${\epsilon }^N={\tilde{\epsilon }}^N+D$. For better model performance and to correct for the turbulent diffusion near walls, a turbulent variable Prandtl number is added of the form, $f_t/{\sigma }_k=1+3.5exp(-{(R{e}_T/150)}^2)$ with $R{e}_T=k^2/\left({\nu }_s\tilde{\epsilon }\right)$ and model constant ${\sigma }_k=1.1$. More details of the form of Equation (19) can be found in Pinho et al. [20] and Resende et al. [24].

The last two terms on the right side of the Equation (19) are:

$$Q^V=\frac{\partial \overline{{\tau }_{ik,p}^{\prime }{u}_i}}{\partial x_k}\mathrm{and}{\epsilon }^V=\frac{1}{\rho }\overline{{\tau }_{ik,p}^{\prime }\frac{\partial u_i}{\partial x_k}},$$

(21)

which are the viscoelastic turbulent transport and the viscoelastic stress work, respectively. They represent the fluctuating viscoelastic turbulent part of the k transport equation and require suitable closure models.

A budget analysis for each term in the k transport equation was performed by Pinho et al. [20] for different regimes of DR. They demonstrated that the magnitude of $Q^V$ has more impact on the overall budget in the IDR, and also developed a closure. In the HDR, the amplitude of $Q^V$ is the same as ${\epsilon }^V$ but has a different location in the buffer layer, in which the effects of $Q^V$ are overcome by turbulent diffusion, thus, revealing negligible effects to overall flow predictions. Masoudian et al. [28] had chosen to neglect the $Q^V$ contributions in the $k–\epsilon –v^2–f$ model and is also not included here as well.

2.3. Transport Equation for the Rate of Dissipation of Turbulent Kinetic Energy

The corresponding governing transport equation for the modified Newtonian rate of dissipation of k is given by,

$$\begin{array}{cc}\hfill \rho \frac{\partial {\tilde{\epsilon }}^N}{\partial t}+\rho U_i\frac{\partial {\tilde{\epsilon }}^N}{\partial x_i}=& \rho \frac{\partial }{\partial x_i}\left[\left({\nu }_s+\frac{f_t{\nu }_T}{{\sigma }_{\epsilon }}\right)\frac{\partial {\tilde{\epsilon }}^N}{\partial x_i}\right]+f_1{C}_{{\epsilon }_1}\frac{{\tilde{\epsilon }}^N}{k}{P}_k\hfill \\ & -f_2{C}_{{\epsilon }_2}\rho \frac{{\left({\tilde{\epsilon }}^N\right)}^2}{k}+\rho E+E_{{\tau }_p},\hfill \end{array}$$

(22)

with

$$E={\nu }_s{\nu }_T(1-f_{\mu }){\left(\frac{{\partial }^2{U}_j}{\partial x_i\partial x_k}\right)}^2.$$

(23)

As mentioned in the previous sub-section, all terms are modelled in the Newtonian context (excluding $E_{{\tau }_p}$). The damping functions of Equation (22) are $f_1=1$ and $f_2=1-0.3exp(-{\left(R{e}_T\right)}^2)$; with model coefficients ${\sigma }_{\epsilon }=1.3$, $C_{{\epsilon }_1}=1.45$ and $C_{{\epsilon }_2}=1.90$.

The last term in Equation (22) is the viscoelastic contribution to the overall ${\tilde{\epsilon }}^N$ balance. This term acts as a Newtonian destruction to the dissipation and is given by,

$$E_{{\tau }_p}=2{\mu }_s\frac{{\mu }_p}{\lambda (L^2-3)}\overline{\frac{\partial u_i}{\partial x_m}\frac{\partial }{\partial x_k}\left\{\frac{\partial }{\partial x_m}\left[f\left(C_{nn}\right)f({\widehat{C}}_{pp})c_{qq}{C}_{ik}\right]\right\}}.$$

(24)

It has non-negligible effects on flow predictions for all DR regimes and thus requires a suitable model.

3. Development of Viscoelastic Closures

In this section, the turbulent viscoelastic cross-correlations that were isolated in the previous section are presented with model closures. The closures are developed on the basis of the DNS data case (19) (

Table 1. Independent Direct Numerical Simulation data for turbulent channel flow of finitely extensible nonlinear elastic-Peterlin (FENE-P) fluids at $\beta =0.9$, with drag reduction (DR) model predictions.

Case	Reference	Rheological Parameters			Drag Reduction (%)
		${\mathit{Re}}_{{\mathit{\tau }}_0}$	${\mathit{Wi}}_{{\mathit{\tau }}_0}$	${\mathit{L}}^2$	DNS	Current Model	Model [32]
(1)	Li et al. [23]	125	25	900	19	20	-
(2)	Li et al. [23]	125	25	3600	22	23	-
(3)	Li et al. [23]	125	25	14,400	24	25	-
(4)	Li et al. [23]	125	50	900	31	30	35
(5)	Li et al. [23]	125	100	900	37	36	39
(6)	Li et al. [23]	125	100	1800	45	43	-
(7)	Li et al. [23]	125	100	3600	56	51	51
(8)	Masoudian et al. [28]	180	25	900	19	19	-
(9)	Li et al. [23]	180	50	900	31	30	34
(10)	Masoudian et al. [28]	180	100	900	38	38	39
(11)	Masoudian et al. [28]	180	100	3600	54	53	51
(12)	Thais et al. [37]	180	116	10,000	64	60	-
(13)	Iaccarino et al. [26]	300	36	3600	33	32	34
(14)	Iaccarino et al. [26]	300	36	10,000	35	35	32
(15)	Iaccarino et al. [26]	300	120	10,000	59	59	58
(16)	Masoudian et al. [30]	395	25	900	19	22	19
(17)	Masoudian et al. [30]	395	50	900	30	30	-
(18)	Masoudian et al. [30]	395	50	3600	38	38	-
(19)	Masoudian et al. [30]	395	100	900	37	37	38
(20)	Masoudian et al. [30]	395	100	3600	48	47	52
(21)	Masoudian et al. [39]	395	100	10,000	55	55	-
(22)	Masoudian et al. [30]	395	100	14,400	61	60	62
(23)	Thais et al. [37]	395	116	10,000	62	60	-
(24)	Li et al. [23]	395	200	14,400	75	69	67
(25)	Masoudian et al. [30]	590	50	3600	39	40	64
(26)	Thais et al. [37]	590	116	10,000	61	59	74
(27)	Thais et al. [38]	1000	50	900	30	33	60

), and then subsequently compared with other DNS data sets for accurate model predictions. The DNS data in Table 1 pertain to all DR regimes with a large variation in rheological parameters and flow viscosity for fully-developed channel flow established by: Li et al. [23]; Thais et al. [37,38]; Masoudian et al. [28,30,39] and Iaccarino et al. [26].

The non-dimensional numbers that define the different DNS data sets are defined as follows: the friction Reynolds number $R{e}_{{\tau }_0}=h{u}_{\tau }/{\nu }_0$ is based on the friction velocity ($u_{\tau }$), the channel half-height (h), the zero shear-rate kinematic viscosity of the solution, which is the sum of the kinematic viscosity of the solvent and polymer (${\nu }_0={\nu }_s+{\nu }_p$); The Weissenberg number $W{i}_{{\tau }_0}=\lambda u_{\tau }^2/{\nu }_0$; and the ratio between the solvent viscosity and the solution viscosity at zero shear rate is $\beta ={\nu }_s/{\nu }_0$.

In the following sub-sections, closures are developed for: the ${NLT}_{ij}$ term of Equation (12) with focus on the dominant ${NLT}_{xx}$ component; a modification of the damping function $f_{\mu }$ (Equation (18)), named $f_{\nu }$, which accounts for the reduction of the Reynolds shear stress due to viscoelastic effects; the viscoelastic stress work, ${\epsilon }^V$ of Equation (19); and the viscoelastic destruction, $E_{{\tau }_p}$, of Equation (22).

3.1. Closure for ${NLT}_{ij}$

The ${NLT}_{ij}$ exact transport equation is greatly simplified based on the DNS analysis of Pinho et al. [20]: Following the transport equation of $\overline{f\left({\widehat{c}}_{mm}\right)c_{kj}\frac{\partial u_i}{\partial x_k}}$ + $\overline{f\left({\widehat{c}}_{mm}\right)c_{ik}\frac{\partial u_j}{\partial x_k}}$, it is assumed that

$$\overline{f\left({\widehat{c}}_{mm}\right)c_{kj}\frac{\partial u_i}{\partial x_k}}+\overline{f\left({\widehat{c}}_{mm}\right)c_{ik}\frac{\partial u_j}{\partial x_k}}\approx f\left(C_{mm}\right)\left(\overline{c_{kj}\frac{\partial u_i}{\partial x_k}}+\overline{c_{ik}\frac{\partial u_j}{\partial x_k}}\right)=f\left(C_{mm}\right){NLT}_{ij}.$$

(25)

The full details of this approximation and the exact transport equation of ${NLT}_{ij}$ can be found in Pinho et al. [20] and Resende et al. [24].

The complete closure of ${NLT}_{ij}$ is presented below and was developed to improve model predictions based on better physical modeling compared with the most recent model developed by Resende et al. [32].

$$\begin{array}{cc}\hfill {NLT}_{ij}& =\overline{c_{kj}\frac{\partial u_i}{\partial x_k}}+\overline{c_{ik}\frac{\partial u_j}{\partial x_k}}\hfill \\ & \approx \underset{I}{\underbrace{f_N{C}_{N1}\frac{\lambda \sqrt{\tilde{L}}{\epsilon }^N}{{\nu }_{0}f\left(C_{mm}\right)}{\delta }_{ij}}}-\underset{II}{\underbrace{f_N^{1/4}{C}_{N2}{M}_{ij}}}+\underset{III}{\underbrace{C_{N3}\frac{k}{{\nu }_0}\sqrt{\frac{\tilde{L}{M}_{nn}}{\dot{\gamma }}}\frac{\frac{\partial U_i}{\partial x_k}\frac{\partial U_j}{\partial x_k}}{{\dot{\gamma }}^2}}},\hfill \end{array}$$

(26)

where $f_N={\nu }_T/{\nu }_0$ is the local eddy viscosity, $\dot{\gamma }=\sqrt{2S_{pq}{S}_{pq}}$ is the shear rate invariant, $\tilde{L}=\sqrt{L^2/900}$ is the normalised maximum extension with the lowest DR, with model constants $C_{N1}=0.11$, $C_{N2}=0.3$ and $C_{N3}=0.3$.

The closure of Equation (26) is modelled in three parts: parts I and $II$ are modeled in the same fashion as the model of Resende et al. [32], part $III$ is greatly improved and is the main contribution to the ${NLT}_{ij}$ closure.

Part I is approached by introducing the Taylor’s longitudinal micro-scale, ${\lambda }_f$, to the relationship between the double correlation of fluctuating strain rates and the turbulent kinetic energy in homogeneous isotropic turbulence. More details can be found in Resende et al. [32], with adjustments $L^{0.42}$ to $\sqrt{L}$ and $f{\left(C_{mm}\right)}^{0.8}$ to $f\left(C_{mm}\right)$.

Part $II$ is primarily responsible for capturing the shear component, ${NLT}_{xy}$. The correlation here is with the exact term, $M_{ij}$ (see Equation (13)), and by the local eddy viscosity, $f_N^{1/4}$. The $L^{0.15}$ variation is removed from the model developed by Resende et al. [32]. The negative part of the ${NLT}_{xx}$ component is also captured here via the $M_{xx}$ term, which according to Dubief et al. [27], is the region where polymers inject energy into turbulence.

Part $III$ is developed to predict the ${NLT}_{xx}$ component which is the dominant term in the trace of ${NLT}_{ij}$, responsible for the stretch of the polymer chains due to turbulent fluctuations. Following the same assumption as Masoudian et al. [29], one can see that ${NLT}_{xx}\sim \overline{u_x^{\prime }{u}_x^{\prime }}\sim k$. In physical terms, the turbulent stretching terms represent the ability of the turbulent fluctuations to act on the polymer chains. This stretching is effective if the polymer shear and maximum extensibility are large enough. So, $\sqrt{\tilde{L}{M}_{nn}/\dot{\gamma }}$ is included here with k. Note that for fully developed channel flow, this term reduces to $\sqrt{\tilde{L}{C}_{xy}}$ which increases proportionally to drag reduction. This new term includes the same physical assumption as Masoudian et al. [28,30], and is simplified from the very complex ad-hoc approach of Resende et al. [32], viz

$$\begin{array}{cc}\hfill {NLT}_{III}^{\mathrm{Resende}}=& f_N^{0.9}exp\left(-\frac{-f_N}{1.05\sqrt{\beta }(10+0.3L+\tilde{L}-{(\tilde{L}-1)}^2)}\right)\hfill \\ & \times {\left(\frac{C_{mm}}{{(\beta /0.9)}^{0.7\beta }}\left(2-{\left[1-exp\left(-\frac{2U_{b}h/{\nu }_s}{3500}\right)\right]}^4\right)\right)}^{0.7}\frac{\frac{d{U}_i}{d{x}_k}\frac{d{U}_j}{d{x}_k}}{\sqrt{\frac{d{U}_p}{d{x}_q}\frac{d{U}_p}{d{x}_q}}}.\hfill \end{array}$$

(27)

The performance of the ${NLT}_{ij}$ closure can be analysed in Figure 1 by comparing the predictions with DNS data case (19) in Table 1, and with the model of Resende et al. [32]. Figure 1a–c plots each normal component of ${NLT}_{ij}$, with the predictions as accurate as the previous model [32]. The new ${NLT}_{xx}$ component is capable of predicting the maximum value and peak location of the destruction effect away from the wall along with the negative part near the wall, but requires a much simpler closure. The closure performance becomes more noticeable at higher Reynolds numbers, in which the polymer extension is largely overestimated previously. The ${NLT}_{yy}$ component is the leading order term in the $C_{yy}$ component away from the wall, which is the dominant contributor to an effective polymeric viscosity. This strongly influences the turbulent dynamics according to Thais et al. [40] and Benzi et al. [41] with their DNS and toy model analysis respectively. This term is represented by the first term in Equation (26), along with ${NLT}_{zz}$. The ${NLT}_{zz}$ component was shown by Pinho et al. [20] to have low impact, and thus ${NLT}_{zz}={NLT}_{yy}$ is an appropriate approximation. The shear component, ${NLT}_{xy}$, can be viewed in Figure 1e, where the predictions omit similar results compared with the previous model [32], but do not require additional $L^2$ variation via $L^{0.15}$.

Overall, all main features of ${NLT}_{ij}$ are well captured such as the peak locations and magnitudes, but with a much simpler closure for the dominant contributor of polymer stretch, ${NLT}_{xx}$. Further, the ${NLT}_{xy}$ and ${NLT}_{yy}$ terms responsible for the polymer shear stress contribution in the momentum balance are featured, which were previously represented ad-hoc with friction velocity dependence [26,28] or misrepresented [30].

3.2. Model for the Modified Damping Function, $f_{\nu }$

There have been many attempts to predict the eddy viscosity reduction as flow viscoelasticity increases for drag-reducing flows. In the case of low-Reynolds $k–\epsilon$ models for FENE-P fluids, this was examined firstly by Pinho et al. [20] for the LDR regime; then later by Resende et al. [24] for the IDR regime. In both cases, there was a consistent reduction in the magnitude of k as DR increased, contrary to the DNS findings [23]. Similar attempts to model a modified damping function were made by Pinho [15]; Cruz et al. [16]; Resende et al. [17] and Tsukahara and Kawaguchi [31] to develop viscoelastic turbulence models using different constitutive equations.

Recently, Resende et al. [32] proposed a modified damping function which was able to predict the correct behavior of the eddy viscosity close to the wall, leading to the appropriate increase for the magnitude of k, and the shift away from the wall into the buffer layer as DR increased. This proposal was founded from the a-priori DNS data analysis by Resende et al. [42], demonstrating the necessary increase to the production of k close to the wall. The model derived by Resende et al. [32] is based on the DNS analysis of Li et al. [23], with an approximation of the form $\mathrm{DR}\sim C_{kk}/L$, giving rise to the correct damping of near-wall eddies as DR increases. In the $k–\epsilon –v^2–f$ models proposed by Iaccarino et al. [26] and Masoudian et al. [28], the near-wall eddy viscosity damping effect is achieved by $v^2$, as ${\nu }_T=C_{\mu }{v}^{2}k/\epsilon$. However, the reduction in $v^2$ is not enough to increase k as given by the DNS data.

The approach by Resende et al. [32] works well in increasing k in the buffer layer, but fails to capture the viscoelastic effects away from the wall, due to the fact that $f_{\mu }^{\mathrm{Previous}}\to 1$ as $y\to h$, which is contrary to the DNS data of Li et al. [22] and the analogous behavior of $v^2$ away from the wall. Therefore an additional model is required to capture the effect of nonequilibrium away from the wall, similarly to the Newtonian model of Park et al. [43]. Benzi et al. [41] demonstrated that the overall effect of polymer stretching is to introduce an effective viscosity proportional to $C_{yy}$, which is dominated by the ${NLT}_{yy}$ component (modeled here with the first term in Equation (26)). An additional term is multiplied to the eddy viscosity to account for the global reduction of eddy structures for increasing DR. This approach is similar to the model of Resende et al. [24] and the study using DNS data of Resende et al. [42] which multiplies the damping function by a factor of $1-g\left(\mathrm{VE}\right)$, where $g\left(\mathrm{VE}\right)$ is a function of the viscoelastic terms, VE.

The final model presented for the modified damping function, $f_{\nu }$, is

$$f_{\nu }=\left(1-A\right){\left[1-exp\left(-\frac{y^{*}}{a_{\mu }\left(1+B/a_{\mu }\right)}\right)\right]}^2,$$

(28)

$$A=C_A{\left(f_N\frac{{\lambda }^2{\tilde{L}}^{3/2}}{f{\left(C_{kk}\right)}^2}\frac{\epsilon }{{\nu }_0}\right)}^{0.3},$$

(29)

$$B=C_B{(C_{kk}-3)}^{1.25}/L,$$

(30)

with model constants $C_A=0.071$ and $C_B=0.44$. An additional contribution in the present model comes from an alternative representation of the dimensionless wall scaling $y^{+}=u_{\tau }/{\nu }_w$, where ${\nu }_w$ is the wall viscosity. The presence of the wall friction velocity poses a problem for flows with re-circulation or reattachment were the friction velocity becomes null at these points, causing floating point errors within computational solvers. Possibilities other than $y^{+}$ that solve this issue are $R{e}_y\equiv ky/{\nu }_0$ or the turbulent Reynolds number, $R{e}_t$. Wallin and Johansson [44] formulated an alternative scaling, $y^{*}$, in terms of $R{e}_y$ so that $y^{*}\approx y^{+}$ for $y^{+}\le 100$ in channel flows. The form proposed is

$$y^{*}=C_{y1}R{e}_y^{1/2}+C_{y2}R{e}_y^2,$$

(31)

where $C_{y1}=2.4$ and $C_{y2}=0.003$. The $R{e}_y$-term is motivated by the fact that the near-wall asymptotic behaviour for $R{e}_y^{1/2}$ is $\sim y^2$. The $R{e}_y^2$-term is artificially introduced to obtain a near linear relation in the buffer region.

The performance of the $f_{\nu }$ closure can be analyzed in Figure 2 by comparing the predictions with DNS data cases (16, 19, 20) for LDR, IDR and HDR respectively in Table 1 and with the model of Resende et al. [32]. The predictions offer significant improvement away from the wall compared to the previous model. The effects can be viewed for the turbulent kinetic energy and the eddy viscosity in the results section, offering improved results for various levels of DR and Reynolds numbers. The $f_{\nu }$ closure more accurately represents the anisotropic effect akin to the $v^2–f$ models of Masoudian et al. [28,30], with the thickening of the buffer layer from the stretched polymers, along with a global reduction with the new closure.

3.3. Development of Closures for ${\epsilon }^V$ and $E_{{\tau }_p}$

The closure model for ${\epsilon }^V$ is approached following the DNS budget analysis of the governing terms proposed by Pinho et al. [20]. In their work, they verified that the double correlation can be neglected with respect to the triple correlation at LDR. This was later confirmed for IDR and HDR by Resende et al. [24] and Masoudian et al. [28], respectively. Pinho et al. [20] extended this analysis and demonstrated that the triple correlation can be decoupled and modeled as a function of ${NLT}_{mm}/2$. Following this, Masoudian et al. [28] confirmed the model capabilities within 5% accuracy for all DR regimes via an extensive pdf study, and is the model used here given by

$${\epsilon }^V\approx \frac{{\nu }_{\mathrm{p}}}{2\lambda }f\left(C_{mm}\right){NLT}_{mm}.$$

(32)

The closure model derived for $E_{{\tau }_p}$ assumes that it depends on the same quantities as the classical Newtonian destruction term of the transport equation of $\epsilon$, but involving a viscoelastic quantity, typically with the viscoelastic stress work used by Resende et al. [24,32] and Masoudian et al. [28,30]. However, as ${\epsilon }^V$ contains a negative part close to the wall via the ${NLT}_{mm}$ contribution, it is not feasible to include ${\epsilon }^V$ in a suitable model for $E_{{\tau }_P}$, based on the DNS analysis of $\epsilon$ being strictly decreasing near the wall for increasing DR.

The closure derived by Resende et al. [32] is complex with $W{i}_{{\tau }_0}$ dependence to force the correct trend in $\epsilon$. Here, a much simpler approach is obtained with dependence through k and some viscoelastic quantities which increases proportional with DR. The closure is given by

$$E_{{\tau }_p}\approx -C_{N4}\frac{{\tilde{\epsilon }}^N}{k}\left[{\nu }_p\sqrt{C_{\mu }{f}_{\mu }}{\tilde{L}}^{3/4}{\left(\frac{k}{{\nu }_0}\right)}^2\right],$$

(33)

with model constant $C_{N4}=0.083$. The effect of Equation (33) on $\epsilon$ predictions can be viewed in the results section for LDR and HDR.

Overall, it is clear that all the developed viscoelastic closures presented in this study perform well compared with DNS data. Most importantly, this was achieved without the use of friction velocity dependence. The simplicity of the governing closures allows easy implementation into 3D codes and can be extended to flows with reattachment when DNS data becomes available.

4. Summary of the Present Model

The governing equations with complete closure models that were developed in the previous sections are presented here.

Momentum equation:

$$\rho \frac{\mathrm{D}{U}_i}{\mathrm{D}t}=-\frac{\partial \overline{p}}{\partial x_i}+\rho \frac{\partial }{\partial x_k}\left[({\nu }_s+{\nu }_T)\frac{\partial U_i}{\partial x_k}\right]+\rho \frac{\partial }{\partial x_k}\left(\frac{{\nu }_p}{\lambda }\left[f\left(C_{nn}\right)C_{ik}-{\delta }_{ik}\right]\right),$$

(34)

where the eddy viscosity is given by

$${\nu }_T=C_{\mu }{f}_{\nu }\frac{k^2}{{\epsilon }^N},$$

(35)

with modified damping function

$$f_{\nu }=\left(1-A\right){\left[1-exp\left(-\frac{y^{*}}{a_{\mu }\left(1+B/a_{\mu }\right)}\right)\right]}^2,$$

(36)

$$A=C_A{\left(f_N\frac{{\lambda }^2{\tilde{L}}^{3/2}}{f{\left(C_{kk}\right)}^2}\frac{\epsilon }{{\nu }_0}\right)}^{0.3},$$

(37)

$$B=C_B{(C_{kk}-3)}^{1.25}/L,$$

(38)

with constants $a_{\mu }=26.5$, $C_A=0.071$ and $C_B=0.44$. $y^{*}$ is given by Equation (31).

Conformation tensor equation:

$$\frac{\mathrm{D}{C}_{ij}}{\mathrm{D}t}-M_{ij}-{NLT}_{ij}=-\frac{1}{\lambda }[f\left(C_{kk}\right)C_{ij}-{\delta }_{ij}],$$

(39)

with

$${NLT}_{ij}\approx f_N{C}_{N1}\frac{\lambda \sqrt{\tilde{L}}{\epsilon }^N}{{\nu }_{0}f\left(C_{mm}\right)}{\delta }_{ij}-f_N^{1/4}{C}_{N2}{M}_{ij}+C_{N3}\frac{k}{{\nu }_0}\sqrt{\frac{\tilde{L}{M}_{nn}}{\dot{\gamma }}}\frac{\frac{\partial U_i}{\partial x_k}\frac{\partial U_j}{\partial x_k}}{{\dot{\gamma }}^2},$$

(40)

where $f_N={\nu }_T/{\nu }_0$ is the local eddy viscosity, $\dot{\gamma }=\sqrt{2S_{pq}{S}_{pq}}$ is the shear rate invariant, $\tilde{L}=\sqrt{L^2/900}$ is the normalised maximum extension with the lowest DR, with model constants $C_{N1}=0.11$, $C_{N2}=0.3$ and $C_{N3}=0.3$.

Transport equation of k:

$$\rho \frac{\mathrm{D}k}{\mathrm{D}t}=\rho \frac{\partial }{\partial x_i}\left[\left({\nu }_s+\frac{f_t{\nu }_T}{{\sigma }_k}\right)\frac{\partial k}{\partial x_i}\right]+P_k-\rho ({\tilde{\epsilon }}^N+D)-\frac{{\nu }_p}{\lambda }f\left(C_{mm}\right)\frac{{NLT}_{mm}}{2},$$

(41)

where $P_k=-\rho \overline{u_i{u}_j}\frac{\partial U_i}{\partial x_j}$ is the rate of production of k.

Dissipation transport equation:

$$\begin{array}{cc}\hfill \rho \frac{\mathrm{D}{\tilde{\epsilon }}^N}{\mathrm{D}t}=& \rho \frac{\partial }{\partial x_i}\left[\left({\nu }_s+\frac{f_t{\nu }_T}{{\sigma }_{\epsilon }}\right)\frac{\partial {\tilde{\epsilon }}^N}{\partial x_i}\right]-f_2{C}_{{\epsilon }_2}\rho \frac{{\left({\tilde{\epsilon }}^N\right)}^2}{k}+\rho E\hfill \\ & +\left(C_{{\epsilon }_1}{P}_k-C_{N4}{\nu }_p\sqrt{C_{\mu }{f}_{\mu }}{\tilde{L}}^{3/4}{\left(\frac{k}{{\nu }_0}\right)}^2\right)\frac{{\tilde{\epsilon }}^N}{k},\hfill \end{array}$$

(42)

with model constant $C_{N4}=0.083$.

The remaining constants are from the Newtonian model and are $C_{{\epsilon }_1}=1.45$, $C_{{\epsilon }_2}=1.90$, $C_{\mu }=0.09$, ${\sigma }_k=1.1$ and ${\sigma }_{\epsilon }=1.3$.

5. Numerical Procedure

This section presents the numerical methods applied in order to examine the viscoelastic turbulence model against the available DNS data identified within the literature. A new finite volume C++ computational solver was developed in the OpenFOAM software by modifying the $k–\epsilon$ sub-class files and introducing the FENE-P viscoelastic quantities such as: the polymer stress to the momentum equation; conformation tensor transport equation; and modified damping function to include elastic effects.

A fully-developed channel flow using half of the channel height, h, is applied given the symmetry of the governing geometry. We assigned 100 cells in the transverse (wall) direction with approximately 10 cells located inside the viscous sublayer. This is to provide mesh independent results, with errors within 0.5% for the mean velocity and the friction factor compared with a very fine mesh, similarly with [30]. The initial state of the simulation is the Newtonian solution until a steady-state solution was reached for each run case, except for HDR where a similar IDR developed case is applied to reduce computational time. Relaxation factors for the additional conformation tensor field are set to 0.2, along with residual control set to ${10}^{-5}$. To improve numerical stability, an artificial diffusion term is added to the RACE of the form, $\kappa {\partial }_k^2{C}_{ij}$, where $\kappa$ denotes a constant, isotropic, artificial numerical diffusivity. In earlier studies [10], the dimensionless artificial numerical diffusivity is taken to be $\kappa /h{u}_{\tau }\sim O\left({10}^{-2}\right)$. Here, $\kappa /h{u}_{\tau }\sim O\left({10}^{-3}\right)$ and has negligible effect on mean values.

A pressure gradient is forced in the stream-wise direction to be unity, with periodic boundary conditions for all other flow fields, mimicking the DNS procedure of Li et al. [22]. No-slip boundary conditions were imposed on the solid wall for the velocity field U, along with k and $\tilde{\epsilon }$ set to zero (or very small, $\sim {10}^{-15}$). A Dirichlet boundary condition for $C_{ij}$ is reported in Appendix A (similar to [26], but for all components), which is imposed within OpenFOAM under the swak4Foam library using the groovyBC functionality developed by Gschaider [45].

When normalizing the governing equations and inherently the various physical quantities, the velocity scale is taken to be the friction velocity (leading to the use of superscript +) and the length scale is the viscous length, $x_i=x_i^{+}{\nu }_0/U_{\tau }$. The conformation tensor is already in dimensionless form.

6. Results and Discussion

Following the numerical procedure proposed in the previous section, the model performance is assessed against a range of different flow and rheological parameters presented in the DNS data within Table 1.

6.1. Analysis of Conformation Tensor

Figure 3 compares the individual components of the conformation tensor with the present model against the model of Resende et al. [32] and selected DNS data covering $L^2$, $W{i}_{{\tau }_0}$ and $R{e}_{{\tau }_0}$ variations (cases 16, 19, 20 and 26 in Table 1). as can be viewed in Figure 3a, the $C_{xx}$ predictions are consistent with the DNS data. The new closure for ${NLT}_{xx}$ (see term III in Equation (26)) is responsible for the improved predictions and can capture the $R{e}_{{\tau }_0}$, $L^2$ and $W{i}_{{\tau }_0}$ variations with much greater simplicity, especially for increased Reynolds number ($R{e}_{{\tau }_0}=590$) compared with the model of Resende et al. [32].

Figure 3b plots the $C_{yy}$ component, showing good agreement with the DNS data and improving upon the most recent model, especially away from the wall. The important feature is the location of the value at the centre-line and the peak location which both show good improvement, especially for higher Reynolds numbers ($R{e}_{{\tau }_0}=590$). The improvements are a result of the new $E_{{\tau }_p}$ closure (see Equation (33)) which directly impacts ${\epsilon }^N$ in the ${NLT}_{yy}$ closure (see term I in Equation (26)). Figure 3c plots the $C_{zz}$ component and shows an under-prediction due to the isotropic assumption used in the model of ${NLT}_{ij}$, however, its impact is not significant.

The model predictions of the $C_{yy}$ term are important in capturing the features of the $C_{xy}$ component. As can be observed in Figure 3d, the model is able to capture the near-wall region, which, according to the findings of Li et al. [22], is the region of high chain dumbbell extension (limited to $y^{+}<50$) where the effect of $C_{xy}$ acts towards the polymer shear stress.

It is evident that the overall predictions of the individual conformation tensor components are improved compared to the model of Resende et al. [32]. This is a result of the new ${NLT}_{ij}$ and $E_{{\tau }_p}$ closures developed in the present work, which allows more scope of predictability and increased numerical stability with simpler closures.

6.2. Analysis of k, $\epsilon$ and ${\nu }_T$

The predicted k profiles are shown in Figure 4a for cases 16 and 19 in Table 1, and Figure 4b for low and high Reynolds number cases (7 and 27). There is reasonable improvement of the profile away from the wall as a result of the new $f_{\nu }$ closure for increasing drag reduction and for various Reynolds numbers.

In Figure 5, the prediction of the dissipation rate are compared with the DNS data of both LDR (case 16) and very HDR (case 22), along with predictions for the $v^2–f$ model of Masoudian et al. [28]. The predictions for LDR are captured well with the DNS for both near and far from the wall. For HDR, there is a significant improvement near the wall compared with the $v^2–f$ model. This is a result of the $E_{{\tau }_p}$ closure formulated (See Equation (33)) which decreases $\epsilon$ as flow viscoelasticity increases. The model of Resende et al. [32] shows similar results to the current model and is not plotted so that the figure is clearer. However, the complexity of the present $E_{{\tau }_p}$ closure model is reduced substantially and removes all friction velocity dependence, but can still predict all the main flow features with good performance.

The local eddy viscosity is plotted in Figure 6a for all ranges of DR. The combined performance of $f_{\nu }$, k and $\epsilon$ gives rise to the predictions shown. We observe a reduction in the eddy structures within the buffer-layer and log-layer for increasing DR, as the DNS suggests. The damping function predicts well this behavior with the near-wall polymer extension via $C_{kk}$ and the global reduction via $(1-A)$.

6.3. Analysis of Velocity Profiles

Figure 6b shows the mean stream-wise velocity profiles for all ranges of DR at $R{e}_{{\tau }_0}=395$. All of the profiles reduce to the linear distribution $u^{+}=y^{+}$ in the viscous sub-layer. Further from the wall, the velocity profiles are well-captured for all ranges of DR.

The model can also predict well a range of Reynolds numbers with varying rheological parameters as can be viewed in Figure 7a. This is extended in Figure 7b for high Reynolds numbers, where there is a significant improvement compared with the model of Resende et al. [32]. This is a result of the new closure model for ${NLT}_{xx}$ which scales well with Reynolds number and with reduced complexity.

The advantage of the current model is the ability to capture all velocity profiles well within the model limits, with more simplicity with regards to model closures and without friction velocity dependence.

7. Conclusions

A viscoelastic turbulence model in fully-developed drag-reducing channel flow is improved, with turbulent eddies modeled under a $k–\epsilon$ representation, along with polymeric solutions described by the finitely extensible nonlinear elastic-Peterlin (FENE-P) constitutive model. A new finite volume C++ computational solver was developed in the OpenFOAM software by modifying the $k–\epsilon$ sub-class files and introducing the FENE-P viscoelastic quantities such as: the polymer stress to the momentum equation; conformation tensor transport equation; and modified damping function to include elastic effects. The model performance is evaluated against a variety of rheological parameters within the DNS data literature, including: friction Reynolds number $R{e}_{{\tau }_0}=125,180,300,395,590,1000$; Wiessenberg number $W{i}_{{\tau }_0}=25,36,50,100,116,200$; and maximum molecular extensibility of the dumbbell chain $L^2=900,1800,3600,$ $\mathrm{10,000},\mathrm{14,400}$. The DNS data case (19) in Table 1 ($R{e}_{{\tau }_0}=395$, DR = 37%) is used for the calibration of the closures developed for the turbulent cross-correlations identified in Section 3. The model is capable of predicting all flow features for low and high Reynolds numbers at all regimes of DR and improves significantly on the model of Resende et al. [32], with its ability to capture higher Reynolds numbers with simpler and physical-based closures.

The main feature is the formulation of the ${NLT}_{ij}$ term which accounts for the interactions between the fluctuating components of the conformation tensor and the velocity gradient tensor. The advantage of the closure is the reduction in the complexity and use of damping functions in the dominant contribution, ${NLT}_{xx}$, modeled here to increase with turbulent kinetic energy as the flow viscoelasticity increases, demonstrating significant improvement with a range of rheological parameters and flow conditions.

Further improvements are developed for the viscoelastic destruction term, $E_{{\tau }_p}$, within the dissipation rate transport equation. Modeled here with dependence on k and viscoelastic quantities, showing the ability to predict $\epsilon$ for low and high drag reduction.

An improved modified damping function, $f_{\nu }$, is also presented, which is able to predict the global reduction of the eddy viscosity and shift away from the wall for increasing viscoelasticity, whilst also improving the profiles of turbulent kinetic energy.

Overall, predictions compare very well with a wide range of DNS data and significantly improves on capturing all flow features with simplicity and performance compared with the most recent $k–\epsilon$ model developed by Resende et al. [32]. The simplicity of the present model allows easy implementation into 3D codes and increases numerical stability. All friction velocity dependence is removed in the present model which is the first of its kind for damping function $k–\epsilon$ models, whose main advantage is the realization of simulations in geometries with reattachment. Future work to extend to this study includes the development of an improved $k-\omega$ model based on the present model [25]. This would require the same concept of the modified damping function developed in this paper to be applied, with capabilities to predict flow behavior in industrially represented geometries such as pipes and constrictions.