###### Abstract
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 [PERSON] 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, [PERSON] et al. [20] developed a low-Reynolds number \(k-\varepsilon\) 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 [PERSON] et al. [21] and [PERSON] 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 [PERSON] 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%). [PERSON] 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-\varepsilon\) model for predicting DR behaviour but demonstrated great versatility and robustness given its application to alternative two-equation models. During this time, a \(k-\varepsilon-\overline{v^{2}}-f\) model for FENE-P fluids in fully developed channel flow was proposed by [PERSON] et al. [26], following the initial studies of [PERSON] 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 [PERSON] et al. [24], but contains only the trace and not the individual components. The model was later improved by [PERSON] 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, \(e_{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 [PERSON] et al. [29], extending on the idea of a correlation between the Reynolds stresses and the \(NLT_{ij}\) components, similar to [PERSON] et al. [19].