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

: A viscoelastic turbulence model in a fully-developed drag reducing channel ﬂow is improved, with turbulent eddies modelled under a k − ε representation, along with polymeric solutions described by the ﬁnitely 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 simpliﬁed viscoelastic closure for the NLT ij term (which accounts for the interactions between the ﬂuctuating 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 τ p , which removes all friction velocity dependence; lastly by an improved modiﬁed 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 ﬂow ﬁelds for low, intermediate and high friction Reynolds numbers, up to high drag reduction without friction velocity dependence.


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 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, ε V yy ) 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 − ε − 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, ε, 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 − ε − 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 − ε 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 − ε 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 Re τ 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 − ε model for FENE-P fluids is proposed, validated for all drag reduction regimes (low, intermediate, high) and up to the largest friction Reynolds number (Re τ 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 τ p , which removes all friction velocity present in the previous k − ε models. The model is assessed against DNS data covering a wide range of flow conditions in terms of the friction Weissenberg number, Wi τ 0 , maximum polymer extension, L 2 , viscosity ratio, β, and friction Reynolds number, Re τ 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.

Governing Equations
The governing equations for incompressible turbulent flow of dilute polymer solutions are the continuity and momentum equations respectively: where the hat represents instantaneous quantities of velocityû i , pressurep, stress tensorτ ij , and fluid density ρ. The stress tensor is the sum of the Newtonian solvent which obeys Newton's law of viscosity, τ s ij = 2µ sŝij , with µ s representing the solvent viscosity coefficient, and polymeric contributions,τ p ik , The kinematic viscosity is used alternatively throughout this study and is defined as ν = µ/ρ. The instantaneous rate of strain tensor,ŝ ij , is defined aŝ The instantaneous polymer contributions are based on the FENE-P rheological dumbbell model [33], with closure given byτ with known as the Peterlin function, andĉ kk is the trace of the instantaneous conformation tensor. The other parameters that are associated with the FENE-P model are: λ, the relaxation time of the polymeric fluid; L 2 , the maximum extensibility of the dumbbell model; and µ p , the polymer viscosity coefficient. The behaviour of the instantaneous conformation tensor follows a hyperbolic differential equation of the form, The Oldroyd's upper convective derivative of the instantaneous conformation tensor is here denoted with ∇ 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,û 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: referred to as the Reynolds-averaged Navier-Stokes (RANS). The Reynolds stress tensor is u i u k and requires a closure model. The Reynolds-averaged polymer stress is τ p ik and written fully as where the additional term on the right requires a closure. The Peterlin function becomes After Reynolds averaging, the instantaneous conformation tensor equation becomes 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 CT 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, 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 CT 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.

Model for the Reynolds Stress Tensor
The Reynolds stress tensor is computed by adopting the Boussinesq turbulent stress strain relationship, where k is the turbulent kinetic energy, S ij is the mean rate of strain tensor and µ T = ρν T is the eddy viscosity. ν T is modelled by the typical isotropic k − ε turbulence model for low Reynolds numbers, which includes a damping function f µ to account for near-wall effects: is the viscous dissipation of k by the Newtonian solvent, and a µ = 26.5. The dimensionless wall scaling is y + = u τ 0 y/ν 0 , where u τ 0 is the friction velocity, y is the distance to the nearest wall, and ν 0 is the sum of solvent and polymer viscosity coefficients (ν 0 = ν s + ν 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).

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, with (20) 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 ε N =ε 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 /σ k = 1 + 3.5 exp(−(Re T /150) 2 ) with Re T = k 2 /(ν sε ) and model constant σ 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: 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 ε 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 − ε − v 2 − f model and is also not included here as well.

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, with As mentioned in the previous sub-section, all terms are modelled in the Newtonian context (excluding E τ p ). The damping functions of Equation (22) The last term in Equation (22) is the viscoelastic contribution to the overallε N balance. This term acts as a Newtonian destruction to the dissipation and is given by, It has non-negligible effects on flow predictions for all DR regimes and thus requires a suitable model.

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), 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 Re τ 0 = hu τ /ν 0 is based on the friction velocity (u τ ), 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 (ν 0 = ν s + ν p ); The Weissenberg number Wi τ 0 = λu 2 τ /ν 0 ; and the ratio between the solvent viscosity and the solution viscosity at zero shear rate is β = ν s /ν 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 µ (Equation (18)), named f ν , which accounts for the reduction of the Reynolds shear stress due to viscoelastic effects; the viscoelastic stress work, ε V of Equation (19); and the viscoelastic destruction, E τ p , of Equation (22).

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 f (ĉ mm )c kj 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].
where f N = ν T /ν 0 is the local eddy viscosity,γ = 2S pq S pq is the shear rate invariant,L = √ 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. (26) is modelled in three parts: parts I and I I are modeled in the same fashion as the model of Resende et al. [32], part I I I is greatly improved and is the main contribution to the NLT ij closure.

The closure of Equation
Part I is approached by introducing the Taylor's longitudinal micro-scale, λ 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 √ L and f (C mm ) 0.8 to f (C mm ). Part I I 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 1/4 N . 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 I I I 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 ∼ u x u x ∼ 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, L M nn /γ is included here with k. Note that for fully developed channel flow, this term reduces to 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 NLT Resende 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].

Model for the Modified Damping Function, f ν
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 − ε 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 DR ∼ C kk /L, giving rise to the correct damping of near-wall eddies as DR increases. In the k − ε − 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 ν T = C µ v 2 k/ε. 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 Previous µ → 1 as y → 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(VE), where g(VE) is a function of the viscoelastic terms, VE.
The final model presented for the modified damping function, f ν , is 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 τ /ν w , where ν 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 Re y ≡ ky/ν 0 or the turbulent Reynolds number, Re t . Wallin and Johansson [44] formulated an alternative scaling, y * , in terms of Re y so that y * ≈ y + for y + ≤ 100 in channel flows. The form proposed is where C y1 = 2.4 and C y2 = 0.003. The Re y -term is motivated by the fact that the near-wall asymptotic behaviour for Re 1/2 y is ∼ y 2 . The Re 2 y -term is artificially introduced to obtain a near linear relation in the buffer region.
The performance of the f ν 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 ν 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.

Development of Closures for ε V and E τ p
The closure model for ε 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 The closure model derived for E τ p assumes that it depends on the same quantities as the classical Newtonian destruction term of the transport equation of ε, 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 ε V contains a negative part close to the wall via the NLT mm contribution, it is not feasible to include ε V in a suitable model for E τ P , based on the DNS analysis of ε being strictly decreasing near the wall for increasing DR.
The closure derived by Resende et al. [32] is complex with Wi τ 0 dependence to force the correct trend in ε. 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 with model constant C N4 = 0.083. The effect of Equation (33) on ε 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.

Summary of the Present Model
The governing equations with complete closure models that were developed in the previous sections are presented here.
Momentum equation: where the eddy viscosity is given by with modified damping function with constants a µ = 26.5, C A = 0.071 and C B = 0.44. y * is given by Equation (31). Conformation tensor equation: with where f N = ν T /ν 0 is the local eddy viscosity,γ = 2S pq S pq is the shear rate invariant,L = √ 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:
where is the rate of production of k. Dissipation transport equation: with model constant C N4 = 0.083. The remaining constants are from the Newtonian model and are C ε 1 = 1.45, C ε 2 = 1.90, C µ = 0.09, σ k = 1.1 and σ ε = 1.3.

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 − ε 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, κ∂ 2 k C ij , where κ denotes a constant, isotropic, artificial numerical diffusivity. In earlier studies [10], the dimensionless artificial numerical diffusivity is taken to be κ/hu τ ∼ O(10 −2 ). Here, κ/hu τ ∼ O(10 −3 ) 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ε set to zero (or very small, ∼ 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, The conformation tensor is already in dimensionless form.

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. 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 , Wi τ 0 and Re τ 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 Re τ 0 , L 2 and Wi τ 0 variations with much greater simplicity, especially for increased Reynolds number (Re τ 0 = 590) compared with the model of Resende et al. [32].   (33)) which directly impacts ε 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.

Analysis of Conformation Tensor
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 τ 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, ε and ν 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 ν 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 τ p closure formulated (See Equation (33)) which decreases ε 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 τ 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 ν , k and ε 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).  Figure 6b shows the mean stream-wise velocity profiles for all ranges of DR at Re τ 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.

Analysis of Velocity Profiles
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.

Conclusions
A viscoelastic turbulence model in fully-developed drag-reducing channel flow is improved, with turbulent eddies modeled under a k − ε 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 − ε 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 Re τ 0 = 125, 180, 300, 395, 590, 1000; Wiessenberg number Wi τ 0 = 25, 36, 50, 100, 116, 200; and maximum molecular extensibility of the dumbbell chain L 2 = 900, 1800, 3600, 10, 000, 14, 400. The DNS data case (19) in Table 1 (Re τ 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 τ p , within the dissipation rate transport equation. Modeled here with dependence on k and viscoelastic quantities, showing the ability to predict ε for low and high drag reduction.
An improved modified damping function, f ν , 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 − ε 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 − ε 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 − ω 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.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The FENE-P equations simplify considerably if we consider 1D, laminar, parallel flow: u x = u z = 0 and u x = u x (y) ≡ u. The system becomes: where f (C kk ) = L 2 . Introducing the Weissenberg number as Wi = λ du dy and solving the system of equations above, one finds the following cubic equation in f (C kk ) ≡ f : which omits one real solution, f R , that satisfies the laminar equations applicable at the wall: where f R = 1 3 B 2 1/3 + (A11) In our numerical simulations, the explicit definition of the wall value using Equations (A6)-(A9) as a Dirichlet boundary condition considerably improves the stability of the solution procedure and is preferred over a zero-flux Neumann boundary condition for the conformation tensor.