1. Introduction
Viscoelastic fluids are ubiquitous in many industrial sectors, including fast-moving consumer goods, food and healthcare, among others, and it is of significant importance that we are able to model these flows correctly in a wide range of geometries. The Oldroyd-B model (Oldroyd Reference Oldroyd1950) is given as
\begin{equation} \boldsymbol{\tau\!}_p + \lambda \overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p = 2\eta_p\boldsymbol{\mathsf{D}}, \end{equation}
where 
$\boldsymbol {\tau\!}_p$ is the polymeric stress, 
$\lambda$ is the viscoelastic relaxation time, 
$\eta _p$ is the polymeric viscosity, 
$\boldsymbol{\mathsf{D}}$ is the rate-of-strain tensor given by 
$\boldsymbol{\mathsf{D}} \equiv 1/2(\boldsymbol {\nabla } \boldsymbol {u}+ \boldsymbol {\nabla } \boldsymbol {u}^{\mathrm {T}})$, and 
$\overset {\kern 0em\triangledown }{\boldsymbol {\tau\!}}_p$ denotes the upper-convected time derivative of the polymeric stress tensor, which is given as 
${\overset {\kern 0em\triangledown }{\boldsymbol {\tau\!}}_p \equiv \mathop {}\!\mathrm {D} \boldsymbol {\tau\!}_p / \mathop {}\!\mathrm {D} t - \boldsymbol {\tau\!}_p \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u} - \boldsymbol {\nabla } \boldsymbol {u}^{\text {T}} \boldsymbol {\cdot } \boldsymbol {\tau\!}_p}$. For many viscoelastic models, including the Oldroyd-B model, the total stress 
$\boldsymbol {\sigma }$ appearing in the momentum equation is related to the polymeric stress by 
$\boldsymbol {\sigma } = \boldsymbol {\tau } - p\boldsymbol{\mathsf{I}} = \boldsymbol {\tau\!}_p + 2\eta _s \boldsymbol{\mathsf{D}} - p\boldsymbol{\mathsf{I}}$, where 
$\boldsymbol {\tau }$ is the extra-stress tensor, 
$p$ is the pressure, and 
$\boldsymbol{\mathsf{I}}$ is the identity tensor. The solvent contribution to the stress is expressed as a Newtonian fluid with viscosity 
$\eta _s$.
Whilst the simplicity of the Oldroyd-B model makes it particularly useful for solving problems analytically (Rajagopal & Bhatnagar Reference Rajagopal and Bhatnagar1995; Qi & Xu Reference Qi and Xu2007; Zhao, Wang & Wei Reference Zhao, Wang and Wei2013; Norouzi et al. Reference Norouzi, Davoodi, Bég and Shamshuddin2018; Ghosh, Mukherjee & Chakraborty Reference Ghosh, Mukherjee and Chakraborty2021; Boyko & Stone Reference Boyko and Stone2022) and testing and validating computational codes (Mompean & Deville Reference Mompean and Deville1997; Duarte, Miranda & Oliveira Reference Duarte, Miranda and Oliveira2008; Habla et al. Reference Habla, Tan, Haßlberger and Hinrichsen2014), it has a number of well-known shortcomings. Likely the most well-known shortcoming is that the elasticity has no limit of extensibility. During steady and homogeneous extensional flow, this causes an unphysical singularity in the extensional viscosity as the strain rate is increased (Bird, Armstrong & Hassager Reference Bird, Armstrong and Hassager1987). In transient and homogeneous extensional flow, however, the extensional viscosity grows exponentially in time, and the singularity is not present. A vast number of viscoelastic models have since been developed to overcome the problems associated with the Oldroyd-B model. Many of these models are derived from micro-structural theories in order to better capture the underlying physics observed during deformation. Two such models are the finitely extensible nonlinear elastic with Peterlin closure (FENE-P; Bird, Dotson & Johnson Reference Bird, Dotson and Johnson1980) and the simplified Phan-Thien–Tanner (sPTT; Phan-Thien & Tanner Reference Phan-Thien and Tanner1977) models.
The original FENE model (Warner Reference Warner1972) is derived using kinetic theory for bead–spring dumbbells in which each polymer molecule is assumed to take the form of two beads connected together by a finitely extensible spring. Therefore, the FENE-P model is most often employed for the modelling of dilute polymer solutions where there is no significant interaction between polymer molecules. The FENE-P model uses a self-consistent pre-averaging approximation, known as the Peterlin approximation, to close the original FENE model (Bird et al. Reference Bird, Dotson and Johnson1980; Keunings Reference Keunings1997). The springs are finitely extensible since the elastic stress increases nonlinearly during deformation as the stretching of the spring approaches its prescribed limit. The sPTT model is derived from a Lodge–Yamamoto type of network theory, where the springs are interconnected via junction points. It is therefore most applicable for concentrated polymer solutions and melts where there are strong interactions between polymer molecules. Under large deformation, the junctions in the sPTT model can be created and destroyed simultaneously, limiting the build-up of elastic stresses and providing finite extensibility.
Whilst the Oldroyd-B model has a constant shear viscosity in steady and homogeneous shear flow, both the FENE-P and sPTT models are shear-thinning. The first normal stress difference 
$N_1 \equiv \sigma _{11} - \sigma _{22}$ grows quadratically with shear rate in the Oldroyd-B model for steady and homogeneous shear flow, but in the FENE-P and sPTT models, it grows quadratically with shear rate only at low shear rates, before exhibiting shear-thinning. In steady and homogeneous extensional flow, the FENE-P and sPTT models exhibit strain-hardening for low strain rates; however, a plateau is reached in the extensional viscosity for higher strain rates due to the finite extensibility. The value of the extensional viscosity at the plateau is proportional to 
$L^2$ (
$1/\epsilon$) for the FENE-P (sPTT) model, where 
$L^2$ and 
$\epsilon$ represent the respective extensibility parameters in the FENE-P and sPTT models.
The FENE-P constitutive model is given in stress tensor form as
\begin{equation} \boldsymbol{\tau\!}_p + \lambda \left(\frac{\overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p}{F(\tau_p)}\right) = 2a\eta_p\boldsymbol{\mathsf{D}}\left(\frac{1}{F(\tau_p)}\right) - a\eta_p \boldsymbol{\mathsf{I}}\,\frac{\mathop{}\!\mathrm{D}}{\mathop{}\!\mathrm{D} t} \left(\frac{1}{F(\tau_p)}\right), \end{equation} where 
$\tau _p \equiv \mathrm {tr}(\boldsymbol {\tau\!}_p)$, or equivalently,
\begin{equation} \left.\begin{gathered} \frac{F(\tau_p)}{a}\,\boldsymbol{\tau\!}_p + \frac{\lambda_1}{a}\, \overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p = 2\eta_p\boldsymbol{\mathsf{D}} - F(\tau_p) \left[\frac{\lambda}{a}\,\boldsymbol{\tau\!}_p + \eta_p\boldsymbol{\mathsf{I}}\right] \frac{\mathop{}\!\mathrm{D}}{\mathop{}\!\mathrm{D} t} \left(\frac{1}{F(\tau_p)}\right),\\ \text{where}\ F(\tau_p) \equiv a + \frac{\lambda}{L^2 \eta_p}\, \mathrm{tr}(\boldsymbol{\tau\!}_p) \quad \text{and} \quad a \equiv \frac{L^2}{L^2 - 3} \end{gathered}\right\} \end{equation}For steady and homogeneous flows, the substantial derivative term in (1.2) and (1.3) is equal to zero, and the FENE-P model can be rewritten as
\begin{equation} \frac{F(\tau_p)}{a}\,\boldsymbol{\tau\!}_p + \frac{\lambda}{a}\, \overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p = 2\eta_p\boldsymbol{\mathsf{D}}, \quad\text{where}\ \frac{F(\tau_p)}{a} = 1 + \frac{\lambda}{aL^2 \eta_p}\,\mathrm{tr}(\boldsymbol{\tau\!}_p). \end{equation}The sPTT model is given as follows:
\begin{equation} F(\tau_p)\,\boldsymbol{\tau\!}_p + \lambda\, \overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p = 2\eta_p\boldsymbol{\mathsf{D}}, \quad \mathrm{where}\ F(\tau_p) \equiv 1 + \frac{\epsilon \lambda}{\eta_p}\,\mathrm{tr}(\boldsymbol{\tau\!}_p). \end{equation}
The scalar function 
$F(\tau _p)$ is then defined on a per-model basis as
\begin{equation} F(\tau_p) \equiv \begin{cases} a + \dfrac{\lambda}{L^2\eta_p}\,\mathrm{tr}(\boldsymbol{\tau\!}_p), & \text{FENE-P}, \\ 1 + \dfrac{\epsilon\lambda}{\eta_p}\,\mathrm{tr}(\boldsymbol{\tau\!}_p), & \text{sPTT}. \end{cases} \end{equation}
The original PTT model employs the Gordon–Schowalter derivative of the polymeric stress 
$\overset {\kern 0em \Box }{\boldsymbol {\tau\!}}_p \equiv \overset {\kern 0em \triangledown }{\boldsymbol {\tau\!}}_p + \zeta (\boldsymbol {\tau\!}_p \boldsymbol {\cdot } \boldsymbol{\mathsf{D}} + \boldsymbol{\mathsf{D}} \boldsymbol {\cdot } \boldsymbol {\tau\!}_p)$, which allows for non-affine transformations between the junction points and the solvent fluid through the slip parameter 
$\zeta$. The sPTT model refers to the case for the PTT model where 
$\zeta = 0$ and so 
$\overset {\kern 0em \Box }{\boldsymbol {\tau\!}}_p = \overset {\kern 0em \triangledown }{\boldsymbol {\tau\!}}_p$. It should also be noted that the sPTT model (1.5) uses a linear term for the destruction of the junctions, as does the original PTT model; however, there have since been modifications to this where the linear term is replaced by exponential (Phan-Thien Reference Phan-Thien1978), or even generalised (Ferrás et al. Reference Ferrás, Morgado, Rebelo, McKinley and Afonso2019) terms, which are believed to help the model perform better under strong deformations. In this study, we will use the sPTT model only with the linear function (1.5), and we will always refer to this as the sPTT model. For clarity, we often use the subscripts FP and sPTT to denote the FENE-P and sPTT models, respectively.
Upon comparison of (1.4) and (1.5), it is observed that with the parameter substitutions 
$\epsilon = 1/L^2$ and 
$\lambda _{{sPTT}} = \lambda _{{FP}} / a$, the FENE-P and sPTT models become mathematically identical for steady and homogeneous flows. The equivalence of these two models was first noted in the study of Cruz, Pinho & Oliveira (Reference Cruz, Pinho and Oliveira2005), who derived analytical solutions for fully developed pipe and channel flows with the FENE-P and sPTT models. Latreche et al. (Reference Latreche, Sari, Kezzar and Eid2021) also established analytical solutions for steady, fully developed, flows of the FENE-P and sPTT models in flat and circular ducts using the aforementioned substitution of parameter values. Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022) then investigated the FENE-P and sPTT models for a range of steady and homogeneous flows, as well as unsteady and inhomogeneous flows. Due to the presence of the Lagrangian derivative term in the stress tensor form of the FENE-P model, significant differences were observed between the FENE-P and sPTT responses for the transient flows. Notably, the FENE-P model produced pronounced shear stress overshoots in start-up shear flow, and during start-up extensional flow, the extensional viscosity grew much more sharply in time for the FENE-P model response than for the sPTT model response. One of the geometries studied by Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022) was the cross-slot. For viscoelastic flows in the cross-slot, the elastic stresses cause a symmetry-breaking instability to occur at a critical 
$Wi$, which has previously been well studied and characterised (Poole, Alves & Oliveira Reference Poole, Alves and Oliveira2007; Rocha et al. Reference Rocha, Poole, Alves and Oliveira2009; Xi & Graham Reference Xi and Graham2009; Afonso, Alves & Pinho Reference Afonso, Alves and Pinho2010; Haward et al. Reference Haward, Ober, Oliveira, Alves and McKinley2012; Cruz et al. Reference Cruz, Poole, Afonso, Pinho, Oliveira and Alves2014; Davoodi, Domingues & Poole Reference Davoodi, Domingues and Poole2019; Davoodi et al. Reference Davoodi, Houston, Downie, Oliveira and Poole2021). Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022) observed that the critical value of 
$Wi$ for the onset of the asymmetry is lower for the FENE-P model than for the sPTT model when relatively low (high) values of 
$L^2\ (\epsilon )$ are used for the FENE-P (sPTT) model, again highlighting that the complex nature of the flow causes a discrepancy between the model responses even though the flow is Eulerian steady. Many industrial processes and flows involve complex geometries that might induce Lagrangian unsteadiness, even for an Eulerian steady flow. Recently, Varchanis et al. (Reference Varchanis, Tsamopoulos, Shen and Haward2022) highlighted that even the Oldroyd-B model exhibits complex rheological behaviour in Lagrangian transient flows, which has significant consequences for, for example, the understanding of pressure drop measurements across a contraction. It is therefore of significant importance to compare and understand how these nonlinear models behave in transient flows.
An ideal way of probing the transient nonlinear response of viscoelastic materials and models, and in particular classifying complex fluids (Hyun et al. Reference Hyun, Kim, Ahn and Lee2002), is with large-amplitude oscillatory shear (LAOS), which has become a widely used technique for characterising nonlinear viscoelasticity experimentally (Leblanc Reference Leblanc2008; Hyun et al. Reference Hyun, Wilhelm, Klein, Cho, Nam, Ahn, Lee, Ewoldt and McKinley2011; Sun et al. Reference Sun, Yang, Wang, Liu, Wang and Tong2011; Szopinski & Luinstra Reference Szopinski and Luinstra2016), theoretically (Gurnon & Wagner Reference Gurnon and Wagner2012; Khair Reference Khair2016; Bae & Cho Reference Bae and Cho2017; Kammer & Castañeda Reference Kammer and Castañeda2020) and numerically (Ewoldt & McKinley Reference Ewoldt and McKinley2010; D'Avino et al. Reference D'Avino, Greco, Hulsen and Maffettone2013; Cordasco & Bagchi Reference Cordasco and Bagchi2016). In small-amplitude oscillatory shear (SAOS), the shear stress response of a material or constitutive model is approximately linear and given by 
$\tau _{p,12} = \gamma _0[G' \sin (\omega t) + G'' \cos (\omega t)]$, where 
$\gamma _0$ and 
$\omega$ are the amplitude and angular frequency of the oscillation, respectively. Here, 
$G'$ and 
$G''$ represent the storage and loss moduli, respectively. Due to the linearity of the shear stress response, SAOS is one of the most popular techniques for extracting information regarding linear viscoelasticity. For example, 
$\lambda$ is very often estimated as the inverse of the frequency at which 
$G'$ and 
$G''$ cross over in a frequency sweep. However, as 
$\gamma _0$ increases, flow-induced micro-structural changes take place during the oscillation (Gilbert & Giacomin Reference Gilbert and Giacomin2016), and the periodic response of the material (or constitutive model) deviates from linearity. This behaviour can then be interpreted in terms of higher-order harmonics in the shear-stress waveform. Therefore, in LAOS, the stress response cannot be reconstructed accurately using a single mode of 
$G'$ and 
$G''$. Multiple frameworks have been developed for quantitative analysis of the nonlinear stress response obtained from LAOS, namely Fourier transform rheology (Wilhelm, Maring & Spiess Reference Wilhelm, Maring and Spiess1998), stress decomposition (Cho et al. Reference Cho, Hyun, Ahn and Lee2005) with Chebyshev analysis (Ewoldt, Hosoi & McKinley Reference Ewoldt, Hosoi and McKinley2008), and a sequence of physical processes (Rogers et al. Reference Rogers, Erwin, Vlassopoulos and Cloitre2011). The LAOS is considered to be especially useful for the purpose of fitting constitutive models to experimental data (Bae & Cho Reference Bae and Cho2015). Calin, Wilhelm & Balan (Reference Calin, Wilhelm and Balan2010) used LAOS to fit the spectrum of the tensorial mobility parameter of the Giesekus model (Giesekus Reference Giesekus1982) to experimental data using an iterative numerical solution. Gurnon & Wagner (Reference Gurnon and Wagner2012) then derived an asymptotic solution for the Giesekus model in oscillatory shear, which they use to fit easily the tensorial mobility parameter to experimental data obtained in the medium-amplitude oscillatory shear (MAOS) regime, where the asymptotic solution is valid. Asymptotic solutions in oscillatory shear have also been derived for the pom-pom model (Hoyle et al. Reference Hoyle, Auhl, Harlen, Barroso, Wilhelm and McLeish2014), the co-rotational Maxwell model (Giacomin et al. Reference Giacomin, Gilbert, Merger and Wilhelm2015), and the White–Metzner model (Merger et al. Reference Merger, Abbasi, Merger, Giacomin, Saengow and Wilhelm2016), among others. Hyun et al. (Reference Hyun, Baik, Ahn, Lee, Sugimoto and Koyama2007) compared the responses of the exponential PTT model, the Giesekus model, and the pom-pom model in MAOS, as well as the experimental MAOS response of linear and branched polymers. For perspective of LAOS tests, the reader is referred to the comprehensive reviews by Hyun et al. (Reference Hyun, Wilhelm, Klein, Cho, Nam, Ahn, Lee, Ewoldt and McKinley2011) and Kamkar et al. (Reference Kamkar2022).
For a purely oscillatory shear flow where the strain rate is uniform in space, the strain 
$\gamma (t)$ and strain rate 
$\dot {\gamma }(t)$ are given by 
$\gamma (t) = \gamma _0 \sin (\omega t)$ and 
$\dot {\gamma }(t) = \gamma _0 \omega \cos (\omega t)$, respectively. We define here the non-dimensional polymeric stress 
$\boldsymbol {\tau\!}_p^* \equiv \boldsymbol {\tau\!}_p / (\gamma _0 \omega [\eta _p +\eta _s])$, the non-dimensional velocity gradient 
$(\boldsymbol {\nabla } \boldsymbol {u})^* \equiv \boldsymbol {\nabla } \boldsymbol {u} / (\gamma _0\omega )$, and the non-dimensional time 
${t^* \equiv t\omega }$. We also define the Weissenberg number as 
$Wi \equiv \lambda \gamma _0 \omega$, and the Deborah number as 
$De \equiv \lambda \omega$. As usual, we define the dimensionless parameter 
$\beta$ as the ratio of the solvent viscosity to the total viscosity (polymeric viscosity plus solvent viscosity), such that 
$\beta \equiv \eta _s / (\eta _s + \eta _p)$. Using these definitions to non-dimensionalise the FENE-P (1.3) and sPTT (1.5) models, and dropping the asterisks upon non-dimensionalisation, we have, respectively,
\begin{align} \frac{F(\tau_p)}{a}\,\boldsymbol{\tau\!}_p &+ \frac{De}{a}\,\frac{\partial}{\partial t}\boldsymbol{\tau\!}_p - \frac{Wi}{a}\,(\boldsymbol{\tau\!}_p \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\text{T}} \boldsymbol{\cdot} \boldsymbol{\tau\!}_p - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tau\!}_p) = 2(1-\beta)\boldsymbol{\mathsf{D}} \nonumber\\ &\quad - F(\tau_p)\left[\frac{Wi}{a}\,\boldsymbol{\tau\!}_p + (1-\beta)\boldsymbol{\mathsf{I}} \right]\left(\frac{De}{Wi}\,\frac{\partial}{\partial t} \left(\frac{1}{F(\tau_p)}\right) + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} \left(\frac{1}{F(\tau_p)}\right) \right) \end{align}and
\begin{equation} F(\tau_p)\,\boldsymbol{\tau\!}_p + De\,\frac{\partial}{\partial t}\boldsymbol{\tau\!}_p - Wi\,(\boldsymbol{\tau\!}_p \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\text{T}} \boldsymbol{\cdot} \boldsymbol{\tau\!}_p- \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\tau\!}_p) = 2(1-\beta)\boldsymbol{\mathsf{D}}, \end{equation}where
\begin{equation} F(\tau_p) \equiv \begin{cases} a + \dfrac{Wi}{L^2(1-\beta)}\,\mathrm{tr}(\boldsymbol{\tau\!}_p), & \text{FENE-P}, \\ 1 + \dfrac{\epsilon\,Wi}{1-\beta}\,\mathrm{tr}(\boldsymbol{\tau\!}_p), & \text{sPTT}. \end{cases} \end{equation} For all of the models discussed in this study, including the Oldroyd-B model, the extra-stress tensor is given as 
$\boldsymbol {\tau } = \boldsymbol {\tau\!}_p + 2\beta \boldsymbol{\mathsf{D}}$. In dimensionless form, it is upon substitution of 
$\epsilon = 1/L^2$ and 
$Wi_{{sPTT}} = Wi_{{FP}} / a$ that the FENE-P and sPTT models become mathematically identical for steady (
$De = 0$) and homogeneous (
$\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {\tau\!}_p = 0$) flows.
With regard to the definitions of 
$De$ and 
$Wi$ for LAOS, Kamani et al. (Reference Kamani, Donley, Rao, Grillet, Roberts, Shetty and Rogers2023) highlighted recently that using a time-independent value of 
$De$ might seem unphysical in some cases since the true ratio of the flow time scale (the inverse of the oscillation frequency) and the material time scale may not necessarily be constant during the oscillation in certain conditions. This requires that 
$De$, according to its physical interpretation, be a time-dependent value rather than constant value. Whilst the FENE-P and sPTT models have constant relaxation times, and constant values of 
$De$ and 
$Wi$ appear naturally from the non-dimensionalisation of the equations for LAOS, the White–Metzner model, on the contrary, contains a strain-rate-dependent relaxation time. In this case, transient values of 
$De$ and 
$Wi$ would appear naturally from the equations. In the FENE-P and sPTT models, one might also think of an ‘effective’ relaxation time based on 
$\lambda$ and 
$F(\tau _p)$. Therefore, we note that, depending on the model or material in question, one may start to question the correct choice of definition for 
$De$ and 
$Wi$ in LAOS, and whether they should be indeed constant or not during an oscillation. However, this is outside of the scope of the current study, and we use only the time-independent values for 
$De$ and 
$Wi$ defined previously.
In the limit 
$De \ (De / a) \rightarrow 0$ and 
$Wi \ (Wi / a) \rightarrow 0$, the sPTT (FENE-P) models, as well as the Oldroyd-B model, reduce to that of a Newtonian fluid. Note that 
$\lim _{Wi \rightarrow 0} (1/F(\tau _p)_{{FP}}) = 1/a$ and 
$\mathop {}\!\mathrm {D} (1/a) /\mathop {}\!\mathrm {D} t = 0$. In the limit 
$De \rightarrow 0$, the response of each model reduces to its respective steady-state response. In the case that 
$F(\tau _p)_{{sPTT}} \rightarrow 1$ for the sPTT model, or 
$F(\tau _p)_{{FP}}/a \rightarrow 1$ for the FENE-P model, the Oldroyd-B model is obtained. Note that for the FENE-P model, 
$F(\tau _p)_{{FP}}/a \rightarrow 1$ is equivalent to 
$F(\tau _p)_{{FP}} \rightarrow a$ and therefore 
$(\partial /\partial t)(1/F(\tau _p)_{{FP}}) \rightarrow 0$ and 
$\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } (1/F(\tau _p)_{{FP}}) \rightarrow 0$, so the last term on the right-hand side of (1.7) vanishes in this limit.
Viscoelastic constitutive models can also be written for the conformation tensor 
$\boldsymbol{\mathsf{A}}$. For dumbbell models such as the FENE-P model, 
$\boldsymbol{\mathsf{A}}$ can be given as 
$\boldsymbol{\mathsf{A}} \equiv \langle \boldsymbol {Q}\boldsymbol {Q} \rangle / Q_{eq}^2$, where 
$\boldsymbol {Q}$ is the end-to-end vector of an individual dumbbell (the angle brackets represent the ensemble average), and 
$Q_{eq}^2$ is the square of the magnitude at equilibrium, given as 
$Q_{eq}^2 \equiv \langle \boldsymbol {Q} \boldsymbol {\cdot } \boldsymbol {Q} \rangle _{eq} / 3$ (Alves, Oliveira & Pinho Reference Alves, Oliveira and Pinho2021). In general, for other viscoelastic models, 
$\boldsymbol {Q}$ might represent the end-to-end vector of polymer chains or subchains (Hoyle & Fielding Reference Hoyle and Fielding2016), rather than the dumbbell vector specifically. The dimensionless FENE-P model is given in conformation tensor form as
\begin{equation} \left.\begin{gathered} De\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - Wi\,(\boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{A}}) =- (F(A)\, \boldsymbol{\mathsf{A}} - a \boldsymbol{\mathsf{I}}),\\ \text{where}\ F(A) \equiv \frac{L^2}{L^2 - \mathrm{tr}(\boldsymbol{\mathsf{A}})}, \end{gathered}\right\} \end{equation}which can also be rewritten as
\begin{equation} \left.\begin{gathered} \frac{De}{a}\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - \frac{Wi}{a}\,(\boldsymbol{\mathsf{A}} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{\mathsf{A}}) =- \left( \frac{F(A)}{a}\,\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}} \right), \\ \text{where}\ \frac{F(A)}{a} = \frac{L^2 - 3}{L^2 - \mathrm{tr}(\boldsymbol{\mathsf{A}})}. \end{gathered}\right\} \end{equation}The sPTT model is given in conformation tensor form as
\begin{equation} \left.\begin{gathered} De\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - Wi\,(\boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} -\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{A}}) =- F(A)\,(\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}}), \\ \text{where}\ F(A) \equiv [1+\epsilon(\mathrm{tr}(\boldsymbol{\mathsf{A}})-3)]. \end{gathered}\right\} \end{equation}
Note that 
$A \equiv \mathrm {tr}(\boldsymbol{\mathsf{A}})$. Then 
$\boldsymbol {\tau\!}_p$ is recovered from the solutions of (1.10)–(1.12) as
\begin{equation} \boldsymbol{\tau\!}_p = \begin{cases} \dfrac{a(1-\beta)}{Wi} \left(\dfrac{F(A)}{a}\,\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}} \right), & \text{FENE-P}, \\ \dfrac{1-\beta}{Wi}\,(\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}}), & \mathrm{sPTT}. \end{cases} \end{equation} Therefore, 
$F(A)$ is also defined on a per-model basis as
\begin{equation} F(A) \equiv \begin{cases} \dfrac{L^2}{L^2 - \mathrm{tr}(\boldsymbol{\mathsf{A}})}, & \text{FENE-P}, \\ 1+\epsilon(\mathrm{tr}(\boldsymbol{\mathsf{A}}-3)), & \text{sPTT}. \end{cases} \end{equation}
As highlighted by Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022), the evolution equation for 
$\boldsymbol{\mathsf{A}}$ in network theory models follows a general form given by
\begin{equation} De\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - Wi\, (\boldsymbol{\mathsf{A}} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{A}}) =- (D(A)\,\boldsymbol{\mathsf{A}}-C(A)\,\boldsymbol{\mathsf{I}}), \end{equation}with
\begin{equation} \boldsymbol{\tau\!}_p = \frac{1-\beta}{Wi}\,(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}), \end{equation}
where 
$D(A)$ and 
$C(A)$ represent, respectively, the rates of destruction and creation of micro-structures. For the sPTT model, 
$D(A) = C(A) = F(A)_{{sPTT}}$. It is therefore observed, given (1.11), that the FENE-P model might be considered as a type of network model in which, under large deformations, the rate of destruction of micro-structures is faster than the rate of creation of micro-structures. Network models with faster destruction rates than creation rates are expected, and have been observed, to exhibit large amounts of elastic recoil (Davoodi et al. Reference Davoodi, Zografos, Oliveira and Poole2022).
Generally, the LAOS response of a viscoelastic material or model can be classified as one of four archetypes: I, strain thinning; II, strain hardening (or strain thickening); III, weak strain overshoot; and IV, strong strain overshoot (Hyun et al. Reference Hyun, Kim, Ahn and Lee2002). Physically, each classification is believed to correspond to a particular type of underlying micro-structural interaction. Sim, Ahn & Lee (Reference Sim, Ahn and Lee2003) investigated numerically the LAOS response of a general network model, and found that the classification of the LAOS response varied depending on the choice of the parameters defining the rates of creation and destruction of junctions. Townsend & Wilson (Reference Townsend and Wilson2018) simulated the LAOS response of a Newtonian solvent with suspended dumbbells, where the dumbbells are implemented in Stokesian dynamics, thus forming a viscoelastic medium. They compare the simulation results for FENE dumbbells with the LAOS response of the FENE-P constitutive model, which they obtain numerically. For 
$De = 0.56$, they observe that the FENE-P constitutive model shows purely strain thinning behaviour, whereas the FENE dumbbell simulations show some weak strain overshoot for the elastic or storage modulus 
$G'$. With increased oscillation frequency, the FENE-P response changed to a type III response where 
$G''$ exhibited a strain overshoot, whereas the FENE dumbbell simulations showed a type I response. Recently, some authors have also used a micro–macro approach for modelling standard FENE dumbbells and FENE-type networks in LAOS using a technique known as the Brownian configuration field method (Gómez-López et al. Reference Gómez-López, Ferrer, Rincón, Aguayo, Chávez and Vargas2019; Vargas et al. Reference Vargas, Gómez-López, Escandón, Mil-Martínez and Phillips2023). In the FENE-type network model response, self-intersecting secondary loops were observed in the viscous Lissajous curves when the rate of destruction of micro-structures was faster than the rate of creation of micro-structures. As already mentioned, in the context of the PTT model framework (1.15), this causes the model to appear more similar in form to the FENE-P model and likely leads to more elastic recoil in the transient model response. Self-intersecting secondary loops, which will be discussed in more detail later, are known to be related specifically to large amounts of elastic recoil (Ewoldt & McKinley Reference Ewoldt and McKinley2010). Ng, McKinley & Ewoldt (Reference Ng, McKinley and Ewoldt2011) performed LAOS experiments with a gluten dough, which they then modelled with a transient network model. The rate of destruction of the junctions was modelled by a term that is essentially a blend between 
$F(A)_{{sPTT}}$ at low stretching and 
$F(A)_{{FP}}$ at high stretching. They also include 
$F(A)_{{FP}}$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship, so the constitutive model represents a FENE-type network model. The model was able to predict at least qualitatively the experimental Lissajous curves; however, the authors note that the stress overshoots were grossly over-predicted. They attribute this to the functional form of the spring function (essentially 
$F(A)_{{FP}}$), and they introduce a modified function, which diverges to infinity before the FENE limit is approached to temper empirically the magnitude of the stress overshoots. Keunings (Reference Keunings1997) shows that these transient stress overshoots in the FENE-P model arise from the pre-averaging Peterlin approximation used to close the original FENE model. The response of the sPTT model in LAOS was obtained and studied recently by Ofei (Reference Ofei2020), who showed that with increasing 
$De$ and 
$Wi$, clearly the sPTT response deviates away from the linear upper-convected Maxwell (UCM)/Oldroyd-B response. However, in this study, no quantitative analysis of the generated waveforms was conducted.
Despite the facts that there has been significant recent interest in the similarities and differences between the FENE-P and sPTT model responses in steady and unsteady (or complex) flows, and that the LAOS responses of these models have been studied independently, there has yet to be an explicit comparison made between the responses of the models in LAOS when the parameters are chosen such that models provide the same steady and homogeneous response. The aim of this study is to compare the responses of the FENE-P and sPTT constitutive models specifically in LAOS, and to understand and highlight any differences observed in the responses.
2. Numerical methodology
2.1. Zero-dimensional modelling
The majority of the results in this study are obtained by solving constitutive equations assuming that 
$\boldsymbol{\mathsf{A}}$, 
$\boldsymbol {\tau\!}_p$, and 
$\dot {\gamma }$ are uniform in space. We denote this approach to solving the equations as the zero-dimensional (0-D) method. This methodology will now be detailed.
For an ideal oscillatory shear flow, the dimensionless constitutive model can be solved using
\begin{equation} \boldsymbol{\nabla} \boldsymbol{u}(t) = \begin{bmatrix} 0 & 0 & 0\\ \dot{\gamma}(t) & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 0 & 0\\ \mathrm{cos}(t) & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}\!. \end{equation}
Since the FENE-P model in stress tensor form (1.3) cannot be expressed easily as a set of ordinary differential equations (ODEs) for an oscillatory shear flow, we solve the models in conformation tensor form. From this point on, we work only with dimensionless variables, and we re-confirm that the asterisks denoting the dimensionless variables have been dropped for brevity. The following system of ODEs is obtained for the time evolution of 
$\boldsymbol{\mathsf{A}}$ according to the FENE-P model,
$$\begin{gather} \frac{\mathop{}\!\mathrm{d} {A}_{11}}{\mathop{}\!\mathrm{d} t} = 2\left(\frac{Wi}{De}\right){A}_{12} \cos(t) - \frac{1}{De} \left(\left[\frac{L^2}{L^2 - ( {A}_{11}+ {A}_{22}+ {A}_{33})}\right]{A}_{11} - a\right)\!, \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{12}}{\mathop{}\!\mathrm{d} t} =\left(\frac{Wi}{De}\right) {A}_{22} \cos(t) - \frac{1}{De} \left(\left[\frac{L^2}{L^2 - ( {A}_{11}+ {A}_{22}+ {A}_{33})}\right]A_{12}\right)\!, \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{22}}{\mathop{}\!\mathrm{d} t} =- \frac{1}{De} \left(\left[\frac{L^2}{L^2 - ( {A}_{11}+ {A}_{22}+ {A}_{33})}\right]{A}_{22} - a \right)\!, \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{33}}{\mathop{}\!\mathrm{d} t} =- \frac{1}{De} \left(\left[\frac{L^2}{L^2 - ( {A}_{11}+ {A}_{22}+ {A}_{33})}\right]{A}_{33} - a \right)\!, \end{gather}$$and according to the sPTT model,
$$\begin{gather} \frac{\mathop{}\!\mathrm{d} {A}_{11}}{\mathop{}\!\mathrm{d} t} = 2 \left(\frac{Wi}{De}\right) {A}_{12} \cos(t) - \frac{1}{De}\, (1+\epsilon({A}_{11}+ {A}_{22}+ {A}_{33} -3))({A}_{11} - 1), \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{12}}{\mathop{}\!\mathrm{d} t} = \left(\frac{Wi}{De}\right) {A}_{22} \cos(t) - \frac{1}{De}\, (1+\epsilon({A}_{11}+ {A}_{22}+ {A}_{33} -3)){A}_{12}, \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{22}}{\mathop{}\!\mathrm{d} t} =- \frac{1}{De}\,(1+\epsilon({A}_{11}+ {A}_{22}+ {A}_{33} -3))({A}_{22}-1), \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{33}}{\mathop{}\!\mathrm{d} t} =- \frac{1}{De}\,(1+\epsilon({A}_{11}+ {A}_{22}+ {A}_{33} -3))({A}_{33}-1), \end{gather}$$ where 
$\boldsymbol {\tau\!}_p$ is recovered from 
$\boldsymbol{\mathsf{A}}$ with (1.13). With an initial condition 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{I}}$ (i.e. 
$\boldsymbol {\tau\!}_p = \boldsymbol{0}$), 
$\mathop {}\!\mathrm {d} A_{22}/\mathop {}\!\!\mathrm {d} t = \mathop {}\!\mathrm {d} A_{33}/\mathop {}\!\!\mathrm {d} t = 0$ at all times for the sPTT model, so 
$A_{22}$ and 
$A_{33}$ remain fixed at unity. However, for the FENE-P model, under large oscillatory deformations, 
$A_{22}$ and 
$A_{33}$ will become time-dependent and lower than unity, although it is still the case that 
$A_{22}=A_{33}$. We note here that for the FENE-P model, 
$A_{22} = A_{33} = a/F(A)_{{FP}}$ in steady-state conditions, so 
$\tau _{p,22} = \tau _{p,33} = 0$ according to (1.13). Therefore, 
$A_{22}$ and 
$A_{33}$ deviate from unity in the FENE-P response in steady shear, even though the corresponding stresses are still zero. In LAOS, however, the unsteadiness of the flow implies that 
$\tau _{p,22}$ and 
$\tau _{p,33}$ are also non-zero.
For all 0-D simulations, we omit the solvent contribution to the stress by setting 
$\beta = 0$ so that we study only the response of the viscoelastic constitutive model itself. Therefore, from here on, the Oldroyd-B model is denoted as the UCM model. For the FENE-P and sPTT models, we performed simulations for five values 
$L^2=1/\epsilon =3.1, 5, 10, 100, 1000$. Equations (2.2) and (2.3) were solved in MATLAB using the ode15s solver, which uses built-in adaptive time stepping. The simulations were run until a steady periodic state was reached. For the results in § 3, data are plotted only for the final oscillation when the system is steady periodic (i.e. the limit cycle).
2.2. One-dimensional modelling
We also use a one-dimensional (1-D) modelling approach by solving both the momentum equation and the constitutive model in a 1-D gap of fluid. This is more representative of an actual shear rheometry experiment in which the velocity gradient can become non-uniform in the gap due to phenomena such as shear banding. To solve the equations in the 1-D approach, we use the method of lines (MOL) technique, in which spatial derivatives of flow variables are discretised (in this case using finite difference approximations).
The top and bottom walls of the gap are parallel to the 
$x$-direction. The first- and second-order spatial derivatives of a scalar variable 
$\phi$ are discretised using a fourth-order finite difference scheme, respectively, as
\begin{equation} \left(\frac{\partial \phi}{\partial y} \right)_i= \begin{cases} \dfrac{-3\phi_{i-1}-10\phi_i+18\phi_{i+1}-6\phi_{i+2}+\phi_{i+3}}{12\varDelta}, & i = 2, \\ \dfrac{\phi_{i-2}-8\phi_{i-1}+8\phi_{i+1}-\phi_{i+2}}{12\varDelta}, & 2 < i < (N_y-1), \\ \dfrac{-\phi_{i-3}+6\phi_{i-2}-18\phi_{i-1}+10\phi_i+3\phi_{i+1}}{12\varDelta}, & i = N_y - 1, \end{cases} \end{equation}and
\begin{equation} \left(\frac{\partial^2 \phi}{\partial y^2} \right)_i = \begin{cases} \dfrac{11\phi_{i-1}-20\phi_i+6\phi_{i+1}+4\phi_{i+2}-\phi_{i+3}}{12\varDelta^2}, & i = 2, \\ \dfrac{-\phi_{i-2}+16\phi_{i-1}-30\phi_i+16\phi_{i+1}-\phi_{i+2}}{12\varDelta^2}, & 2 < i < (N_y-1) \\ \dfrac{-\phi_{i-3}+4\phi_{i-2}+6\phi_{i-1}-20\phi_i+11\phi_{i+1}}{12\varDelta^2}, & i = N_y - 1, \end{cases} \end{equation}
where the index 
$i$ denotes the node number in a uniformly discretised domain with 
$N_y$ elements. Here, 
$\varDelta$ is the distance between neighbouring cells, given as 
$\varDelta \equiv y_i - y_{i-1}$. The (dimensional) velocity at the top boundary, 
$u_{N_y}(t)$, is varied according to 
${u_{N_y} = \gamma _0\omega H \cos (\omega t)}$, where 
$H$ is the gap height, and 
$\gamma _0\omega$ represents the strain-rate amplitude. For non-dimensionalisation, 
$H$ is used for the length scale, 
$\gamma _0\omega H$ is used for the velocity scale, and the time is still non-dimensionalised with 
$\omega$. Then 
$De$ and 
$Wi$ are defined as they are for the 0-D approach.
Assuming that the only non-zero velocity component is in the 
$x$-direction, and the flow is uniform in the 
$x$-direction, the resulting system of partial differential equations (PDEs) to be solved can be expressed in dimensionless form as
$$\begin{gather} Re \left(\frac{\partial u}{\partial t} \right)_i = \left(\frac{\partial \tau_{p, 12}}{\partial y}\right)_i + \beta \left(\frac{\partial^2 u}{\partial y^2} \right)_i, \end{gather}$$
$$\begin{gather}\left(\frac{\partial}{\partial t} A_{11} \right)_i = (\mathcal{F}_{11})_i, \end{gather}$$
$$\begin{gather}\left(\frac{\partial}{\partial t} A_{12} \right)_i = (\mathcal{F}_{12})_i, \end{gather}$$
$$\begin{gather}\left(\frac{\partial}{\partial t} A_{22} \right)_i = (\mathcal{F}_{22})_i, \end{gather}$$
$$\begin{gather}\left(\frac{\partial}{\partial t} A_{33} \right)_i = (\mathcal{F}_{33})_i, \end{gather}$$
where the Reynolds number 
$Re$ is defined as 
$Re \equiv \rho H^2\gamma _0\omega /(\eta _s + \eta _p)$, and the tensor 
$\mathcal {F}$ is the right-hand side of the constitutive model when expressed for the time derivative in conformation tensor form. First- and second-order derivatives are replaced with the discretised forms in (2.4) and (2.5), which turns the system of PDEs into a system of ODEs. For the momentum equation, 
$\tau _{p,12}$ is computed from 
$\boldsymbol{\mathsf{A}}$ using (1.13), then its gradient is discretised with (2.4). For the 1-D MOL modelling, we could not omit totally the solvent contribution to the stress due to stability issues. We therefore used 
$\beta = 1/1001$, which was found to be large enough to stabilise the simulations, but small enough so that the results were essentially insensitive to the value of 
$\beta$ in the ranges of 
$De$ and 
$Wi$ investigated. This is shown in the supplementary material available at https://doi.org/10.1017/jfm.2023.977. We also enforce true creeping flow so that inertia is neglected (i.e. the left-hand side of (2.6a) is zero).
For the spatial resolution, we used 
$N_y = 128$ , which proved sufficiently accurate to ensure that the results were independent of 
$N_y$. This is also shown in the supplementary material. At the bottom wall, the velocity was fixed at zero. The components 
$\boldsymbol{\mathsf{A}}$ were extrapolated linearly at the top and bottom boundaries. Simulations were initiated with 
$u = 0$ and 
$\boldsymbol{\mathsf{A}} = \boldsymbol{\mathsf{I}}$. We integrated the resulting system of equations in MATLAB using the adaptive-step ODE solver ode15s, which can solve systems of differential-algebraic equations (DAEs) using the mass matrix approach. We simulated the flow until the limit cycle was reached. Again, only data for the limit cycle are presented in § 3.
3. Results and discussion
For all of the results except those presented in § 3.5, the 0-D method is used (with 
$\beta = 0$) to obtain the solutions. In §§ 3.1 and 3.2, we investigate the model responses in LAOS with the parameter substitutions 
$\epsilon = 1/L^2$, 
$Wi_{{UCM}} = Wi_{{sPTT}} = Wi_{{FP}}/a$ and 
$De_{{UCM}} = De_{{sPTT}} = De_{{FP}}/a$. We present the results for various values of 
$De/a$ (
$De$) and 
$Wi/a$ (
$Wi$) for the FENE-P (sPTT or UCM) model. In § 3.3, we investigate the responses of ‘toy’ models to help to explain the observations from §§ 3.1 and 3.2. In § 3.4, we analyse the model responses using the sequence of physical processes methodology. Finally, in § 3.5, we use 1-D MOL modelling to assess whether the constitutive models are prone to shear banding in LAOS. Throughout much of this study, we present the results by showing the Lissajous–Bowditch curves. For the shear stress, these are displayed as plots of 
$\tau _{p,12}$ versus 
$\gamma$, and plots of 
$\tau _{p,12}$ versus 
$\dot {\gamma }$. The former is referred to as the elastic projection, and the latter is referred to as the viscous projection. The resulting patterns are often presented in Pipkin (or 
$De$ and 
$Wi$) space. For a more detailed overview of Lissajous–Bowditch plots, the reader is referred to the review by Hyun et al. (Reference Hyun, Wilhelm, Klein, Cho, Nam, Ahn, Lee, Ewoldt and McKinley2011).
3.1. Scaling of the Lissajous curves
Figures 1 and 2 show, respectively, the viscous and elastic projections of the Lissajous–Bowditch plots in the 
$De/a$ (
$De$)–
$Wi/a$ (
$Wi$) space for the UCM (black solid lines), FENE-P (yellow to red solid lines) and sPTT (cyan to blue dashed lines) models with varying values of 
$L^2 = 1/\epsilon$ for the FENE-P and sPTT models. The FENE-P and sPTT models deviate from the UCM model at high 
$Wi/a$ (
$Wi$) and particularly at low (high) values of 
$L^2$ (
$\epsilon$). For the majority of the plots (except the four in the upper right quadrant), the responses of the FENE-P and sPTT models to uniform oscillation are practically identical. However, for the upper right quadrant, the responses of the FENE-P and sPTT models differ from each other, particularly for the lower (higher) values of 
$L^2$ (
$\epsilon$), which matches the observations of Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022) for start-up shear flow.
Figure 1. Viscous Lissajous–Bowditch plots in 
$De/a$ (
$De$)–
$Wi/a$ (
$Wi$) space for the FENE-P (sPTT) model. Black curves represent the UCM response. Black numbers in each plot represent the maximum value of 
$\tau _{p,12}$ in the UCM response, since the 
$y$-axis is scaled differently in each plot.
Figure 2. Elastic Lissajous–Bowditch plots in 
$De/a$ (
$De$)–
$Wi/a$ (
$Wi$) space for the FENE-P (sPTT) model. Black curves represent the UCM response. Black numbers in each plot represent the maximum value of 
$\tau _{p,12}$ in the UCM response, since the 
$y$-axis is scaled differently in each plot.
In order to investigate how the responses of the models in LAOS scale with the model parameters, we look first at the way in which the model responses scale under steady simple shear flow (SSSF). For SSSF, the non-dimensional FENE-P model (1.7) yields the following solution for the (polymeric) shear stress:
\begin{equation} \frac{2}{(1-\beta)^2}\left(\frac{{Wi}}{aL}\right)^2 \tau_{p,12}^3 + \tau_{p, 12} = (1-\beta) \dot{\gamma}. \end{equation}
For constant 
$\beta$, the solution for 
$\tau _{p,12}$ in SSSF then evidently depends on the parameter 
$Wi / (aL)$. Oliveira & Pinho (Reference Oliveira and Pinho1999) discussed this scaling, but for the sPTT model response instead, when they derived analytical solutions for fully developed channel flow and showed that the solution scaled with the parameter 
$Wi\,\sqrt {\epsilon }$. With the aforementioned substitution of model parameters, the scaling parameters of both the FENE-P and sPTT models are the same for SSSF. The existence of the scaling parameter 
$Wi/(aL)$ (
$Wi\,\sqrt {\epsilon }$) for the FENE-P (sPTT) models in SSSF was also shown by Oliveira, Coelho & Pinho (Reference Oliveira, Coelho and Pinho2004) and Latreche et al. (Reference Latreche, Sari, Kezzar and Eid2021). A recent study by Yamani & McKinley (Reference Yamani and McKinley2022) showed, analytically, that the SSSF response of the FENE-P model scales with the dimensionless parameter 
$Wi / L$, which differs from the parameter used in this study, 
$Wi / (aL)$. However, in the version of the FENE-P model used in Yamani & McKinley (Reference Yamani and McKinley2022), the value of 
$a$ was assumed to be unity. The different versions of the FENE-P model that appear in the literature have been presented and discussed by Alves et al. (Reference Alves, Oliveira and Pinho2021) and Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022).
Figures 3 and 4 show the viscous and elastic Lissajous–Bowditch plots, respectively, in the 
$Wi/(aL)$ (
$Wi\,\sqrt {\epsilon }$)–
$De/a$ (
$De$) space for the FENE-P (sPTT) model. For the sPTT model, scaling the curves with 
$Wi\,\sqrt {\epsilon }$ causes the responses for the various values of 
$\epsilon$ to become universal (for 
$\beta = 0$). Whilst this is expected in SSSF due to the form of the analytical solution (3.1), it may not be immediately obvious why this is also the case in LAOS. However, this can be shown by considering the following system of equations for the sPTT model (acknowledging that 
$A_{22} = A_{33} = 1$ due to the fact that 
$F(A)$ is on the outside of the 
$(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ term in the constitutive model),
$$\begin{gather} \frac{\mathop{}\!\mathrm{d} {A}_{11}}{\mathop{}\!\mathrm{d} t} = 2 \left(\frac{Wi}{De}\right) {A}_{12} \cos(t) - \frac{1}{De}\,(1+\epsilon({A}_{11}-1))({A}_{11} - 1), \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{12}}{\mathop{}\!\mathrm{d} t} = \left(\frac{Wi}{De}\right) \cos(t) - \frac{1}{De}\,(1+\epsilon({A}_{11}-1))({A}_{12}), \end{gather}$$
$$\begin{gather}\tau_{p,12} = \frac{1-\beta}{Wi}\,A_{12}, \end{gather}$$
and introducing new variables 
$x = \epsilon (A_{11} - 1)$ and 
$y = A_{12}/Wi$, which gives
$$\begin{gather} De\,\frac{\mathop{}\!\mathrm{d}\kern 0.06em x}{\mathop{}\!\mathrm{d} t} = 2 (Wi\,\sqrt{\epsilon})^2 y \cos(t) - (1+x)x, \end{gather}$$
$$\begin{gather}De\,\frac{\mathop{}\!\mathrm{d} y}{\mathop{}\!\mathrm{d} t} = \cos(t)- (1+x)y, \end{gather}$$
$$\begin{gather}\tau_{p,12} = (1-\beta)y. \end{gather}$$
Figure 3. Viscous Lissajous–Bowditch plots in 
$Wi/(aL)$ (
$Wi\,\sqrt {\epsilon }$)–
$De/a$ (
$De$) space for the FENE-P (sPTT) model. Blue numbers in each plot represent the maximum value of 
$\tau _{p,12}$ in the sPTT response, since the 
$y$-axis is scaled differently in each plot. Plots in the black dashed box are shown at a larger scale in figure 5.
Figure 4. Elastic Lissajous–Bowditch plots in 
$Wi/(aL)$ (
$Wi\,\sqrt {\epsilon }$)–
$De/a$ (
$De$) space for the FENE-P (sPTT) model. Black curves represent the UCM response. Blue numbers in each plot represent the maximum value of 
$\tau _{p,12}$ in the sPTT response, since the 
$y$-axis is scaled differently in each plot.
Thus for constant values of 
$\beta$ (in our case 
$\beta = 0$), the system depends only on the parameters 
$De$ and 
$Wi\,\sqrt {\epsilon }$.
The FENE-P LAOS response does not scale universally with 
$Wi/(aL)$ under uniform oscillation, as its SSSF response does. For 
$De / a < 0. 1$, there is practically no difference in the curves for the various values of 
$L^2$ since the flow is approaching SSSF. For 
$De / a \geq 0.1$ and 
$Wi / (aL) \geq 1$, the difference in the FENE-P response with varying 
$L^2$ becomes significant. However, the FENE-P response does seem to become universal at constant values of 
$De/a$ and 
$Wi/(aL)$ for large enough values of 
$L^2$. This is highlighted in figure 5 by the fact that the responses for 
$L^2 = 100$ and 
$L^2 = 1000$ are practically identical. There are at least two potential reasons why the FENE-P response is not universal for constant values of 
$De/a$ and 
$Wi/(aL)$ at low values of 
$L^2$. One is that the functional form of 
$F(A)$ in the FENE-P model is different to that in the sPTT model (see (1.14)). Another is that the position of 
$F(A)$ in the conformation tensor form of the FENE-P model is different to that in the sPTT model (i.e. on the inside of the brackets in the recoil term rather than on the outside of the brackets). If the latter is the cause of the difference in the scaling behaviour between the sPTT and FENE-P models, then it should likely be the case that a universal solution exists for the FENE-CR model, which is given under LAOS as
$$\begin{gather} De\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - Wi\,(\boldsymbol{\mathsf{A}} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{\mathsf{A}}) =- F(A)_{{FP}}\,( \boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}} ), \end{gather}$$
$$\begin{gather}\boldsymbol{\tau\!}_p = \frac{1-\beta}{Wi}\,F(A)_{{FP}}\,(\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}}). \end{gather}$$
Noting that 
$A_{22} = A_{33} = 1$ in the FENE-CR model due to the difference in the position of 
$F(A)$ compared to the FENE-P model, the system of equations for the FENE-CR model under LAOS is given as
$$\begin{gather} \frac{\mathop{}\!\mathrm{d} {A}_{11}}{\mathop{}\!\mathrm{d} t} = 2 \left(\frac{Wi}{De}\right) {A}_{12} \cos(t) - \frac{1}{De} \left(\frac{L^2}{L^2 - A_{11} - 2}\right)(A_{11} - 1), \end{gather}$$
$$\begin{gather}\frac{\mathop{}\!\mathrm{d} {A}_{12}}{\mathop{}\!\mathrm{d} t} = \left(\frac{Wi}{De}\right) \cos(t) - \frac{1}{De} \left(\frac{L^2}{L^2 - A_{11} - 2} \right) A_{12}, \end{gather}$$
$$\begin{gather}\tau_{p,12} = \frac{1-\beta}{Wi} \left(\frac{L^2}{L^2 - A_{11} - 2} \right) A_{12}. \end{gather}$$
Introducing the new variable 
$x = (A_{11}-1)/L^2$, the extensibility function can be rewritten as 
$(1-x-3/L^2)^{-1}$, which, for 
$L^2 \gg 3$ becomes 
$(1-x)^{-1}$. The system of equations for the FENE-CR model for 
$L^2\gg 3$ can then be rewritten, using also 
$y = A_{12}/Wi$, as
$$\begin{gather} De\,\frac{\mathop{}\!\mathrm{d}\kern 0.06em x}{\mathop{}\!\mathrm{d} t} = 2 \left(\frac{Wi}{L}\right)^2 y \cos(t) - \left(\frac{1}{1-x}\right)x, \end{gather}$$
$$\begin{gather}De\,\frac{\mathop{}\!\mathrm{d} y}{\mathop{}\!\mathrm{d} t} = \cos(t)-\left(\frac{1}{1-x}\right)y, \end{gather}$$
$$\begin{gather}\tau_{p,12} = (1-\beta) \left(\frac{1}{1-x}\right)y. \end{gather}$$
Figure 5. Zoomed Lissajous–Bowditch plots (viscous projection) for (a) 
$De/a \ (De) = 1$ and 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon }) = 10$, and (b) 
$De/a \ (De) = 10$ and 
$Wi/(aL)\ (Wi\,\sqrt {\epsilon }) = 10$. See figure 3 for the legend.
Therefore, the FENE-CR model has universal solutions for constant values of 
$De$ and 
$Wi/L$ only in the case that 
$L^2 \gg 3$, and the breakdown of the universality can be caused solely by a change in the functional form of the extensibility function 
$F(A)$ without a change in its position in the constitutive model. This result alone cannot explain the scaling of the FENE-P response due to the fact that both the form of 
$F(A)$ and its position in the model (i.e. 
$F(A)\,(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ versus 
$F(A)\,\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}$) are different to those in the sPTT model. We will explore the effect of the position of 
$F(A)$ in the constitutive model on the scaling of the response further in § 3.3. We also note here that 
$(1-x)^{-1}$ can be expanded as 
$1 + x + O(x^2)$. Therefore, for 
$\epsilon = 1/L^2$, the evolution of 
$\boldsymbol{\mathsf{A}}$ for the FENE-CR model becomes mathematically identical to that for the sPTT model in the MAOS regime in the case 
$L^2 \gg 3$. In the MAOS regime, where the response is weakly nonlinear, terms of 
$O(x^2)$ can be neglected in the expansion of 
$F(A)$. There is, however, a difference in the stress response due to the presence of the extensibility function in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship in the FENE-CR model. We highlight this further in Appendix A.
3.2. Comparing the LAOS response of the FENE-P and sPTT models
One of the primary aims of this study is to assess the difference between the sPTT and FENE-P models in oscillatory shear flow. As mentioned, in SSSF, both models exhibit identical responses in terms of 
$\boldsymbol {\tau\!}_p$ with 
$\epsilon = 1/L^2$ and 
$Wi_{{sPTT}} = Wi_{{FP}} / a$, and differences in the model responses arise only due to transients in the flow.
In figures 3 and 4, it is observed that the sPTT and FENE-P models (with the aforementioned substitution of parameters) have, naturally, identical responses in oscillatory shear flow for 
$De / a \ (De) \rightarrow 0$, which corresponds to the system approaching SSSF. Both models also exhibit identical responses for 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon }) \rightarrow 0$; however, in this case, both models reduce to the UCM model. For large values of both 
$De/a \ (De)$ and 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon })$, there is a significant difference between the responses of the two models. One such difference is that the FENE-P response exhibits sharp overshoots in the shear stress, which are often observed in the responses of strongly nonlinear viscoelastic models, such as those derived specifically for worm-like micelles. These can be observed in particular in the plots for 
$De / a = 1$ and 
$Wi / (aL) = 10$, which are shown at a larger scale in figure 5(a). Similar stress overshoots were also observed in the FENE-P model response during start-up shear flow (Davoodi et al. Reference Davoodi, Zografos, Oliveira and Poole2022). In LAOS, the pronounced shear stress overshoots manifest as self-intersecting secondary loops in the viscous Lissajous curves. The criterion for the presence of the self-intersecting secondary loops is that the gradient of the decomposed elastic stress with respect to the strain is negative at 
$\gamma = 0$ (or 
$\dot {\gamma } = 1$), indicating that elasticity is being relieved through recoil faster than it is being accumulated through increased rates of deformation (Ewoldt & McKinley Reference Ewoldt and McKinley2010). These secondary loops are associated with strongly nonlinear viscoelastic responses and have been observed both in experimental LAOS data and in the responses of several viscoelastic constitutive models, as well as, as already mentioned, from simulations of network models.
In order to quantify the deviation of the FENE-P and sPTT models from the UCM response, we define a function 
$G(A)$ given by
\begin{equation} G(A) \equiv\begin{cases} \dfrac{F(A)_{{FP}}}{a} -1, & \text{FENE-P}, \\ F(A)_{{sPTT}}- 1, & \text{sPTT}. \end{cases} \end{equation} Note that with the relevant expression relating 
$\boldsymbol {\tau\!}_p$ to 
$\boldsymbol{\mathsf{A}}$ for each model (see (1.13)), it is the case that 
$F(A)_{{FP}} = F(\tau _p)_{{FP}}$ and 
$F(A)_{{sPTT}} = F(\tau _p)_{{sPTT}}$. We point this out just to clarify that as 
$G(A) \rightarrow 0$, both the conformation tensor and stress tensor forms of the models asymptote towards the UCM model.
Figure 6 shows 
$G(A)$ against the dimensionless strain rate for the corresponding Lissajous curves presented in figure 3. For both models, 
$G(A)$ is approximately 0 for values of 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon }) \leq 0.1$ at any value of 
$De/a\ (De)$. Consequently, the dimensionless shear stress response shown in the corresponding Lissajous curves is essentially that of the UCM model. For 
$De/a \ (De) = 0.01$, even when 
$G(A)$ increases and the model responses become increasingly nonlinear, for each model, the solution is universal for various values of 
$L^2$ and 
$\epsilon$ at constant 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon })$, which has been explained already by the fact that the steady-shear response of the dimensionless model contains only the scaling parameter 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon })$. Note here that the response of 
$G(A)$ is also the same for both the sPTT and FENE-P models, despite the fact that the functions 
$F(A)_{{FP}}/a$ and 
$F(A)_{{sPTT}}$ are not explicitly equivalent for 
$\epsilon = 1/L^2$, 
$Wi_{{sPTT}} = Wi_{{FP}} / a$ and 
$De_{{sPTT}} = De_{{FP}} / a$. Therefore, in SSSF, i.e. 
$De \rightarrow 0$, it is the case that 
$\mathrm {tr}(\boldsymbol{\mathsf{A}})_{{sPTT}}= (F(A)_{{FP}}/a) \times \mathrm {tr}(\boldsymbol{\mathsf{A}})_{{FP}}$, which is also implied from (1.13) if 
$(\boldsymbol {\tau\!}_p)_{{FP}} = (\boldsymbol {\tau\!}_p)_{{sPTT}}$. Thus in steady and homogeneous flows, 
$\mathrm {tr}(\boldsymbol{\mathsf{A}})$ for each model differs by a factor of 
$F(A)_{{FP}}/a$, highlighting the difference in the physical interpretation of the polymeric stress from 
$\boldsymbol{\mathsf{A}}$ in each model. For the higher values of 
$De/a \ (De)$ and 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon })$, it is observed again that the sPTT solution for 
$G(A)$ is universal for constant 
$De$ and 
$Wi\,\sqrt {\epsilon }$ with varying 
$\epsilon$, but the FENE-P solution for 
$G(A)$ is universal only for constant 
$De/a$ and 
$Wi/(aL)$ at large values of 
$L^2$. It is also observed that there is significant correspondence between the results in figures 6 and 3. Notably, the overshoots in 
$G(A)$ appear to correspond at least qualitatively with the overshoots in 
$\tau _p$. This will be discussed in more detail next.
Figure 6. Plots of 
$G(A)$ versus 
$\dot {\gamma }$ in 
$Wi/(aL)$ (
$Wi\,\sqrt {\epsilon }$)–
$De/a$ (
$De$) space for the FENE-P (sPTT) model. In each plot, the 
$y$-axis runs from 
$-1$ to 6.31. The 
$x$-axis runs from 
$-1.04$ to 1.04. The black solid line shows 
$y = 0$ on each plot.
For 
$De/a$ and 
$Wi/(aL) = 10$, we present in figure 7 three-dimensional (3-D) plots of 
$G(A)$, 
$A_{12}$, 
$A_{22}$ and 
$\tau _{p, 12}$ against 
$\gamma$ and 
$\dot {\gamma }$ for the FENE-P model response in LAOS (with 
$L^2 = 100$). We present this to highlight more clearly the mechanisms responsible for the distinct behaviour exhibited by the FENE-P model in transient flows. We split the curve into four regions (note that only one half of the curve is shown in terms of the range of 
$\dot {\gamma }$). It is important here to refer back to (2.3) and (1.13) to explain the evolution of the shear stress. The regions are specified in chronological order (i.e. the system moves in time from Region I to Region IV).
Figure 7. Three-dimensional plots showing (a) 
$G(A)$, (b) 
$A_{12}$, (c) 
$A_{22}$, (d) 
$\tau _{p, 12}$, in 
$\gamma$–
$\dot {\gamma }$ space for the FENE-P response (
$L^2 = 100$) for 
$De/a = 1$ and 
$Wi/(aL) = 10$. Only positive 
$\dot {\gamma }$ are displayed for convenience. Regions denoted with roman numerals are discussed in the text.
In Region I, 
$A_{12}$ is negative and the rate of deformation is increasing, meaning that both the growth of 
$A_{12}$ due to deformation and the elastic recoil are acting in the same direction, causing a positive rate of change in time of 
$A_{12}$. The rate of change of 
$A_{12}$ is governed by a balance between growth due to deformation and reduction due to elastic recoil. In the sPTT model, 
$A_{22} = 1$, so the growth of 
$A_{12}$ due to increasing rates of deformation is proportional to the strain rate. In the FENE-P model, however, 
$A_{22}$ is time-dependent for large values of 
$Wi/(aL)$, so the growth of 
$A_{12}$ due to increasing deformation rates is nonlinear (this will be discussed in more detail in § 3.3). Whilst 
$G(A)$ is decreasing in Region I, which reduces the degree of the elastic recoil, the growth of 
$A_{22}$ is relatively large, so ultimately 
$A_{12}$ grows nonlinearly in Region I. Despite the fact that 
$A_{12}$ grows in Region I, 
$\tau _{p,12}$ grows only slightly (and seemingly linearly) due to the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship and the fact that 
$F(A)$ is decreasing.
Region II is characterised by a sharp increase in 
$G(A)$ as 
$\mathrm {tr}(\boldsymbol{\mathsf{A}}) \rightarrow L^2$, which causes 
$\tau _{p,12}$ to increase rapidly but also limits somewhat the growth of 
$A_{12}$ since the elastic recoil increases. In Region II, 
$A_{22}$ also goes through a maximum and begins to decrease, which causes a reduction in the growth of 
$A_{12}$ with increasing rates of deformation, exacerbating the recoil effect. In Region III, 
$G(A)$ is large enough that the elastic recoil now exceeds the growth of 
$A_{12}$ due to the increasing rate of deformation, which drives a negative rate of change of 
$A_{12}$ for increasing strain rates. In this region, 
$G(A)$ also decreases as 
$\mathrm {tr}(\boldsymbol{\mathsf{A}})$ decreases, so 
$\tau _{p,12}$ decreases rapidly due to 
$A_{12}$ and 
$G(A)$ both decreasing (and the fact that 
$F(A)$ is present in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship).
At the start of Region IV, 
$A_{12}$ remains fairly constant, which is a consequence of the balance of the elastic recoil with the building of elasticity due to the rate of deformation. This might be thought of as being representative of the system approaching steady state. Further into Region IV, there is a decrease in 
$G(A)$, which drives a reduction in the shear stress. However, as the deformation rate decreases, there is no significant decrease in 
$A_{12}$ due to the fact that 
$A_{22}$ is increasing and 
$G(A)$ is decreasing, which acts ultimately to slow the recoil of 
$A_{12}$ when the rate of deformation is decreased.
Considering the previous analysis, we then present 3-D plots of 
$G(A)$, 
$A_{12}$, 
$A_{22}$ and 
$\tau _{p, 12}$ against 
$\gamma$ and 
$\dot {\gamma }$ for the FENE-P (
$L^2 = 100$) and sPTT model responses for 
$De/a \ (De) = 1$ and 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon }) = 10$ in figure 8. Note that the sharp overshoots in 
$G(A)$ are not observed for the sPTT model response, and also note the significant differences in the 
$A_{22}$ responses of each model. Here, 
$A_{12}$ is generally much larger for the sPTT model than for the FENE-P model, but the values of 
$\tau _{p, 12}$ are fairly similar for large parts of the oscillation due to (1.13). This analysis unravels exactly where the differences arise between the sPTT and FENE-P models in unsteady shear flows. It should be noted also that the unsteadiness here can be Eulerian or Lagrangian in nature, since the only difference between the two models written in stress tensor form (1.7) and (1.8) is a Lagrangian derivative term on the right-hand side of the FENE-P model. This is not so easily observed when the models are expressed for the conformation tensor. This has significant consequences when these constitutive models are used to model flows in complex geometries that at first might seem steady-state due to the Eulerian steadiness, but might be Lagrangian unsteady (or inhomogeneous).
Figure 8. Three-dimensional plots showing (a) 
$G(A)$, (b) 
$A_{12}$, (c) 
$A_{22}$, (d) 
$\tau _{p, 12}$ in 
$\gamma$–
$\dot {\gamma }$ space for the FENE-P response (
$L^2 = 100$) and the sPTT response for 
$De/a \ (De) = 1$ and 
$Wi/(aL) \ (Wi\,\sqrt {\epsilon }) = 10$.
3.3. Toy sPTT/FENE-P models
In this subsection, we focus on identifying and investigating more closely the differences in the LAOS responses of the sPTT and FENE-P models, and particularly the origins of the pronounced stress overshoots in the FENE-P response. To do this, we define ‘toy’ models by manipulating slightly the standard sPTT and FENE-P models. In this subsection, we do not focus on the aforementioned substitution of parameters to equate the FENE-P and sPTT responses for steady and homogeneous flows, as we just compare qualitatively the responses of each toy model. We employ only the 0-D modelling approach to obtain the solutions, and we also use again 
$\beta = 0$ for all results.
For the toy sPTT model, we start with the generic network model given in (1.15) and adjust the rates of micro-structural destruction 
$D(A)$ and creation 
$C(A)$ as
$$\begin{gather} D(A) = 1 + \alpha \epsilon [\mathrm{tr}(\boldsymbol{\mathsf{A}})-3], \end{gather}$$
$$\begin{gather}C(A) = 1 + (1-\alpha) \epsilon [\mathrm{tr}(\boldsymbol{\mathsf{A}})-3]. \end{gather}$$
Thus for 
$\alpha = 0.5$, we recover the sPTT model (with the value for 
$\epsilon$ halved), and for 
$\alpha = 1$, we recover a model with a similar form to the FENE-P model but without 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship. Therefore, in a way, 
$\alpha$ controls the position of 
$F(A)$ in the recoil term of the constitutive model. The LAOS response for the toy sPTT model with 
$\epsilon = 1/100$ at 
$De = 0.5$ and 
$Wi = 200$ is presented for various values of 
$\alpha$ in figure 9. As 
$\alpha$ is increased, the stress overshoots, and self-intersecting secondary loops become significantly more pronounced, as is expected for systems that exhibit large rates of micro-structural destruction (Davoodi et al. Reference Davoodi, Zografos, Oliveira and Poole2022; Vargas et al. Reference Vargas, Gómez-López, Escandón, Mil-Martínez and Phillips2023), and the response begins to appear qualitatively similar to the FENE-P response. This suggests that the primary reason for the FENE-P response exhibiting large stress overshoots in transient flows is that the extensibility function is multiplied by 
$\boldsymbol{\mathsf{A}}$ instead of 
$(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$, and it is not due to the fact that the extensibility function appears in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship, or due to the difference in the natures of 
$F(A)_{{FP}}$ and 
$F(A)_{{sPTT}}$. We also show this more explicitly using a toy FENE-P model where the evolution equation for 
$\boldsymbol{\mathsf{A}}$ is unchanged (given by (1.11)), but 
$\boldsymbol {\tau\!}_p$ is given by
\begin{equation} \boldsymbol{\tau\!}_p = \frac{a(1-\beta)}{Wi} \left(\left( \frac{F(A)_{{FP}}}{a} \right)^b \boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}}\right)\!, \end{equation}
Figure 9. Viscous Lissajous curves for the toy sPTT model (
$\epsilon = 1/100$) with various values of 
$\alpha$ between 0.5 (for which 
$C(A) = D(A)$) and 1 (for which 
$C(A) = 0$). Here, 
$De = 0.5$ and 
$Wi = 200$, 
$y$-axis limits are shown by the numbers adjacent to the ends of the axes, and 
$x$-axes run from 
$-1.04$ to 1.04.
where 
$0 \leq b \leq 1$. When 
$b = 1$, the original FENE-P model is obtained, and when 
$b = 0$, the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship reverts to the original form given by the Kramers relation (Kramers Reference Kramers1944), and that used for the sPTT model. The viscous Lissajous curves for the toy FENE-P model with 
$L^2 = 100$ at 
$De = 1$ and 
$Wi = 100$ with varying values of 
$b$ are shown in figure 10. The stress overshoots, and self-intersecting loops are observed for all values of 
$b$, indicating that the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship does not explicitly cause pronounced shear stress overshoots in the transient response. As 
$b \rightarrow 0$, the region directly after the stress overshoot becomes significantly flatter, and the Lissajous curves resemble more those of constitutive models derived for worm-like micellar systems. We also note for both figures 10 and 9 that the scale of the 
$y$-axis in each plot is not fixed. The parameters 
$\alpha$ and 
$b$ do affect significantly the maximum stresses reached in the oscillation. In this sense, the stress overshoots that we are discussing here are relative.
Figure 10. Viscous Lissajous curves for the toy FENE-P model (
$L^2 = 100$) with various values of 
$b$ between 1 and 0. Here, 
$De = 1$ and 
$Wi = 100$, 
$y$-axis limits are shown by the numbers adjacent to the ends of the axes, and 
$x$-axes run from 
$-1.04$ to 1.04.
In the standard form of the sPTT model (
$\alpha = 0.5$), the growth term for 
$A_{12}$ due to the rate of deformation is equal to 
$(Wi/De) \cos (t)$ since 
$A_{22} = 1$. As 
$\alpha \rightarrow 1$, there is significant deviation of 
$A_{22}$ from unity during the oscillation, hence deviation of the growth term from a pure cosine wave. Since the time rate of change of 
$A_{12}$ is governed by a balance between the growth and the elastic recoil, the deviation of 
$A_{22}$ from unity contributes significantly to the occurrence of pronounced stress overshoots and self-intersecting secondary loops in LAOS. Ultimately, the origin of this behaviour lies in the use of the upper-convected time derivative combined with the specific positioning of 
$F(A)$ within the constitutive model (i.e. 
$F(A)\,(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ versus 
$F(A)\,\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}$). Figure 11 shows 
$A_{22} \cos (t)$ during one oscillation for the toy sPTT model, with 
$\epsilon = 1/100$ at 
$Wi = 200$ and 
$De = 0.5$ for various values of 
$\alpha$. Clear overshoots of 
$A_{22}\cos (t)$ are observed at approximately 
$t/t_p = 0.3$ and 
$0.8$ as 
$\alpha \rightarrow 1$.
Figure 11. Plots of 
$A_{22} \cos (t)$ versus time for the toy sPTT model (
$\epsilon = 1/100$) for 
$De = 0.5$ and 
$Wi = 200$ during one oscillation with varying 
$\alpha$. The inset highlights the point of the overshoot in 
$A_{22}$ when 
$\alpha \rightarrow 1$.
To highlight more clearly the role of 
$A_{22}$ in the generation of the pronounced stress overshoots and self-intersecting secondary loops, we define the growth and recoil terms for the evolution of 
$A_{12}$ as
$$\begin{gather} Q_g = \left(\frac{Wi}{De}\,A_{22}\right) \cos(t), \end{gather}$$
$$\begin{gather}Q_r = \left(\frac{1}{De}\right)F(A)\,A_{12}. \end{gather}$$
Note that these are just, respectively, the first and second terms on the right-hand sides of (2.2a) and (2.3b). Figure 12 shows 
$Q_g$, 
$Q_r$ and 
$A_{12}$ for the toy sPTT model with 
$\epsilon = 1/100$ when 
$\alpha = 0.5$ and 
$1$, during a quarter of the oscillation period. During this quarter of the period, 
$\dot {\gamma }$ is increasing from 0 to 1. Naturally, stress overshoots are observed in this region for the case when 
$Q_r > Q_g$, which is seemingly the case for the toy sPTT model when 
$\alpha \rightarrow 1$ and thus 
$C(A)=1$. Whilst it is difficult to assess explicitly the role that 
$A_{22}$ plays in the evolution of 
$Q_g$ and 
$Q_r$ due to the strong coupling in the equations, it is clearly the case that the sharp decrease in 
$Q_g$, which can be caused only by a change in 
$A_{22}$ for fixed values of 
$De$ and 
$Wi$, occurs before any observed decrease in 
$Q_r$, which allows for a region where 
$Q_r$ is significantly larger than 
$Q_g$. In figure 12(b), the point of intersection of 
$Q_g(\dot {\gamma })$ and 
$Q_r(\dot {\gamma })$ of course defines the exact position for the maximum of 
$A_{12}$ (and hence 
$\tau _{p, 12}$) associated with the stress overshoot. As 
$\dot {\gamma } \rightarrow 1$, there is a large region where 
$Q_g \approx Q_r$, which as mentioned might represent the system approaching steady shear.
Figure 12. Plots of 
$Q_g$, 
$Q_r$, and 
$A_{12}$ versus 
$\dot {\gamma }$ in one quarter of an oscillation for the toy sPTT model (
$\epsilon = 1/100$) for 
$De = 0.5$ and 
$Wi = 200$, with (a) 
$\alpha = 0.5$ and (b) 
$\alpha = 1$. The dashed vertical line shows the crossover point of 
$Q_g$ and 
$Q_r$, and hence also the maxima of 
$A_{12}$.
Figure 13 shows 
$Q_g$, 
$Q_r$ and 
$A_{12}$ for the FENE-P model (or toy FENE-P model with 
$b=1$) with 
$L^2 = 100$ at 
$De = 1$ and 
$Wi = 100$. Again, the decrease in 
$Q_g$, caused by the time-dependence of 
$A_{22}$, is significant and is primarily responsible for the generation of the pronounced stress overshoots. Also, 
$Q_r$ grows much more sharply for the FENE-P model than for the toy sPTT model due to the difference between 
$F(A)_{{FP}}$ and 
$F(A)_{{sPTT}}$. The overshoots in the FENE-P response are likely exacerbated somewhat by this. We should note that if 
$F(A)_{{sPTT}}$ is replaced by 
$F(A)_{{FP}}$ in the original sPTT model, then stress overshoots can still occur even when 
$A_{22} = 1$ solely due to the increased nonlinearity of 
$Q_r$. However, these overshoots are significantly smaller than those observed when 
$Q_g$ is nonlinear. This replacement of 
$F(A)_{{sPTT}}$ by 
$F(A)_{{FP}}$ in the sPTT model essentially corresponds to a toy FENE-CR model, which we explore in Appendix A.
Figure 13. 
$Q_g$, 
$Q_r$, and 
$A_{12}$ versus 
$\dot {\gamma }$ in one quarter of an oscillation for the FENE-P model (
$L^2 = 100$) for 
$De = 1$ and 
$Wi = 100$. Dashed vertical line shows the cross-over point of 
$Q_g$ and 
$Q_r$, and hence also the maxima of 
$A_{12}$.
As discussed in § 3.1, the standard sPTT model has universal solutions both in steady shear and in LAOS (for constant values of 
$De$) for constant values of 
$Wi\,\sqrt {\epsilon }$. However, the FENE-P model LAOS response has universal solutions only for constant values of 
$Wi/(aL)$ for large enough values of 
$L^2$. In § 3.1, we show that replacing 
$F(A)_{{sPTT}}$ with 
$F(A)_{{FP}}$ in the sPTT model causes a breakdown of this universality, and therefore that the existence of universal solutions is dependent on the specific form of 
$F(A)$ used. Here, we use the toy models to also investigate the effect of the positioning of 
$F(A)$ on the scaling of the LAOS response (i.e. 
$F(A)\,(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ versus 
$F(A)\,\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}$). In figure 14, we show the LAOS response of the toy sPTT model with 
$\alpha = 1$ (such that 
$C(A) = 1$) for 
$De = 0.5$ and 
$Wi\,\sqrt {\epsilon } = 20$ for various values of 
$\epsilon$. The solution seems to be universal only for small enough values of 
$\epsilon$, which is strongly reminiscent of the FENE-P model response presented in § 3.1. Note that 
$\alpha = 1$ makes the toy sPTT model appear in a similar form to the FENE-P model in the network model framework. As highlighted by Tanner (Reference Tanner2006) and Davoodi et al. (Reference Davoodi, Zografos, Oliveira and Poole2022), the generic network model (1.15) can be written for steady and homogeneous flows in stress tensor form as
\begin{equation} D(\tau_p)\,\boldsymbol{\tau\!}_p + Wi\,\overset{\kern0em\triangledown}{\boldsymbol{\tau\!}}_p =2 (1-\beta)\boldsymbol{\mathsf{D}} - \frac{D(\tau_p)\,(1-\beta)}{Wi} \left(1 - \frac{C(\tau_p)}{D(\tau_p)}\right) \boldsymbol{\mathsf{I}}. \end{equation}
Note here that both the time and the velocity gradient are made dimensionless by 
$\dot {\gamma }$, so the entire upper-convected time derivative is multiplied by 
$Wi$. For 
$D(A) \neq C(A)$, the last term on the right-hand side of (3.12) is non-zero. For the toy sPTT model with 
$\alpha = 1$, 
$D(\tau _p) = F(\tau _p)_{{sPTT}}$ and 
$C(\tau _p) = 1$. In this case, the solution of (3.12) in SSSF is given by the system
$$\begin{gather} \left[1 + \frac{Wi\,\epsilon}{1-\beta}\,(\tau_{p,11}+\tau_{p,22}+\tau_{p,33})\right] \tau_{p,11} -2\,Wi\,\dot{\gamma}\tau_{p,12}=-\epsilon(\tau_{p,11}+\tau_{p,22}+\tau_{p,33}), \end{gather}$$
$$\begin{gather}\left[1 + \frac{Wi\,\epsilon}{1-\beta}\,(\tau_{p,11}+\tau_{p,22}+\tau_{p,33})\right] \tau_{p,12} -Wi\,\dot{\gamma}\tau_{p,22}= (1-\beta)\dot{\gamma}, \end{gather}$$
$$\begin{gather}\left[1 + \frac{Wi\,\epsilon}{1-\beta}\,(\tau_{p,11}+\tau_{p,22}+\tau_{p,33})\right] \tau_{p,22} =-\epsilon(\tau_{p,11}+\tau_{p,22}+\tau_{p,33}), \end{gather}$$
$$\begin{gather}\left[1 + \frac{Wi\,\epsilon}{1-\beta}\,(\tau_{p,11}+\tau_{p,22}+\tau_{p,33})\right] \tau_{p,33} =-\epsilon(\tau_{p,11}+\tau_{p,22}+\tau_{p,33}). \end{gather}$$
Figure 14. Viscous Lissajous curves for the toy sPTT model for 
$\alpha = 1$, 
$De = 0.5$ and 
$Wi\,\sqrt {\epsilon } = 20$, with varying 
$\epsilon$.
For finite values of 
$Wi\,\epsilon$ (or 
$Wi\,\sqrt {\epsilon }$), but in the limit that 
$\epsilon \rightarrow 0$, the term 
$-\epsilon (\tau _{p,11}+\tau _{p,22}+\tau _{p,33})$ on the right-hand side of the diagonal components approaches zero, in which case 
${\tau _{p, 22} = \tau _{p, 33} \rightarrow 0}$ and (3.13) reduces to
$$\begin{gather} \left[1 + \frac{Wi\,\epsilon}{(1-\beta)}\,\tau_{p,11}\right] \tau_{p,11}-2\,Wi\,\dot{\gamma}\tau_{p,12}=0, \end{gather}$$
$$\begin{gather}\left[1 + \frac{Wi\,\epsilon}{(1-\beta)}\,\tau_{p,11}\right] \tau_{p,12} = (1-\beta)\dot{\gamma}, \end{gather}$$ which is the solution for the original sPTT model. As discussed in § 3.1, the solution to (3.14) depends only on the parameter 
$Wi\,\sqrt {\epsilon }$ (something which is also seen by introducing the new variable 
$x = Wi\,\epsilon \tau _{p,11}/(1-\beta )$ into (3.14)). It highlights that the breakdown of the universal scaling with 
$Wi\,\sqrt {\epsilon }$ can be caused, at least in SSSF, by setting 
$D \neq C$ even when 
$F(A)_{{sPTT}}$ is used for the extensibility function. Therefore, considering also the scaling analysis for the FENE-CR model presented in § 3.1, the breakdown of the universal scaling for low values of 
$L^2$ in the FENE-P LAOS response is likely a combined effect of the functional form of 
$F(A)$ and its position in the conformation tensor form of the constitutive model (i.e. 
$F(A)\,(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ versus 
$F(A)\,\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}}$). We note briefly that for the FENE-P model, the extra term on the right-hand side of (3.12) is multiplied by the Lagrangian derivative of 
$1/F(\tau _p)$ due to the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship, which is why this term does not affect the SSSF solution for the FENE-P model.
For clarity, we will summarise this subsection briefly. When 
$D(A) \neq C(A)$ in the network model framework, temporal changes in 
$A_{22}$ cause both nonlinear growth and nonlinear recoil of 
$A_{12}$ simultaneously. This causes a region where 
$Q_r$ is significantly larger than 
$Q_g$, which manifests as pronounced stress overshoots in the Lissajous curves. The presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship does not seem to have much effect on the relative stress overshoots. The sPTT model response scales universally due to both the specific functional form of 
$F(A)$ used in the model and the fact that 
$F(A)$ is on the outside of the brackets in the recoil term (or 
$D(A) = C(A)$ in the network model framework). The FENE-P response does not scale universally in LAOS due to the specific form of 
$F(A)$ used and the fact that 
$F(A)$ is on the inside of the brackets in the recoil term (or 
$D(A) \neq C(A)$ in the network model framework).
3.4. Sequence of physical processes
Previously, we have discussed the responses of the FENE-P and sPTT models predominantly in terms of the nature of the underlying ODEs being solved. Here, we analyse the responses of the models in terms of the physical phenomena being represented by the models. To do this, we employ the sequence of physical processes (SPP) analysis (Rogers et al. Reference Rogers, Erwin, Vlassopoulos and Cloitre2011), which will now be outlined briefly.
Whilst the moduli 
$G'$ and 
$G''$ are time-independent during SAOS, and their values can be found from simple Fourier analysis, these moduli are transient during LAOS for a nonlinear viscoelastic material or constitutive model. The SPP framework identifies, at each instant, these transient moduli, denoted as 
$G'_t(t)$ and 
$G''_t(t)$, by utilising the Frenet–Serret theorem. In this approach, each point in the oscillation is given by a position vector in the strain-rate–stress space 
${\boldsymbol {P}(t)=\langle \gamma (t), \dot {\gamma }(t), \tau _{p,12}(t)\rangle }$, noting that 
$\beta = 0$ in this study so 
$\tau_{p,12} = \tau_{12}$. Three additional vectors are then defined at each position on the Lissajous curve, which are used to identify the transient moduli. The tangent vector 
${\boldsymbol {T}(t) = \dot {\boldsymbol {P}}(t)/|\dot {\boldsymbol {P}}(t)|}$ points in the direction of the instantaneous trajectory (the overdot here denoting total differentiation with respect to time). The normal vector 
$\boldsymbol {N}(t) = \dot {\boldsymbol {T}}(t)/|\dot {\boldsymbol {T}}(t)|$ points to the centre of the curvature of the instantaneous trajectory. Finally, the binormal vector 
$\boldsymbol {B}(t) = \boldsymbol {T}(t)\times \boldsymbol {N}(t)$ points in the direction normal to the plane in which the instantaneous trajectory sits (the osculating plane). The transient moduli are then defined (noting that 
$\omega =1$ in our case due to the equations being solved in non-dimensional form) using the components of 
$\boldsymbol {B}$ as
$$\begin{gather} G'_t(t) = \frac{\mathop{}\!\mathrm{d} \tau_{p,12}}{\mathop{}\!\mathrm{d} \gamma} =- \frac{B_{\gamma}(t)}{B_{\tau_{p,12}}(t)}, \end{gather}$$
$$\begin{gather}G''_t(t) = \frac{\mathop{}\!\mathrm{d} \tau_{p,12}}{\mathop{}\!\mathrm{d} \dot{\gamma}} =- \frac{B_{\dot{\gamma}}(t)}{B_{\tau_{p,12}}(t)}. \end{gather}$$
The stress response, for 
$\beta = 0$, is then reconstructed as
\begin{equation} \tau_{p,12}(t) = G'_t\gamma + G''_t\dot{\gamma} + \tau_{p,12}^d, \end{equation}
where 
$\tau _{p,12}^d$ is the displacement stress. For a more detailed explanation of the SPP framework, the reader is referred to Rogers et al. (Reference Rogers, Erwin, Vlassopoulos and Cloitre2011), Rogers (Reference Rogers2012, Reference Rogers2017) and Lee & Rogers (Reference Lee and Rogers2017).
The time-dependent behaviour of 
$G'_t$ and 
$G''_t$ during the (half) oscillation informs us of underlying physical phenomena occurring in the stress response. Increasing values of 
$G'_t$ represent intra-cycle strain-hardening, whilst decreasing values of 
$G'_t$ represent intra-cycle strain-softening. Similarly, for the viscous modulus, increasing values of 
$G''_t$ represent intra-cycle shear-thickening, whereas decreasing values of 
$G''_t$ represent intra-cycle shear-thinning. Negative instantaneous values of 
$G'_t$ can be thought of as representing elastic recoil, and negative instantaneous values of 
$G''_t$ can be thought of as representing viscous backflow (Rogers Reference Rogers2017; Choi, Nettesheim & Rogers Reference Choi, Nettesheim and Rogers2019; Donley, Bantawa & Gado Reference Donley, Bantawa and Gado2022). In this subsection, we perform the SPP analysis for the sPTT (FENE-P) responses for a fixed 
$De \ (De/a) = 0.5$ with 
$L^2 = 1/\epsilon = 100$. We also investigate the responses of the toy models outlines in § 3.4. The SPP freeware (https://publish.illinois.edu/rogerssoftmatter/freeware/) for MATLAB was used for the analysis, which was kindly provided to us by the developers. Central differencing was used for differentiation of the stress response, and a single limit cycle period with 401 time points was used for the analysis.
Figure 15 shows the 3-D Lissajous curves for the sPTT and FENE-P models as well as the Cole–Cole plots (
$G''_t$ versus 
$G'_t$). For the sPTT model, the response manifests as deltoids in the Cole–Cole plots, which are commonly observed experimentally for a range of viscoelastic materials, including doughs (Park & Rogers Reference Park and Rogers2020; Erturk, Rogers & Kokini Reference Erturk, Rogers and Kokini2022). As 
$Wi\,\sqrt {\epsilon }$ is increased, the deltoids increase in size, which physically might represent an increase in the degree of micro-structural change during the oscillation (Park & Rogers Reference Park and Rogers2020). For the FENE-P response, as 
$Wi$ is increased (note that 
$L^2$ is fixed at 
$100$, and we do not assume that the solution truly scales with 
$Wi/(aL)$), the Cole–Cole plots resemble instead arrowhead shapes that are significantly larger in size than the deltoids observed in the sPTT response (see the scales of the axes).
Figure 15. The SPP analysis for the sPTT and FENE-P model responses at 
$De \ (De/a) = 0.5$ and 
$L^2 = 1/\epsilon = 100$. (a,c) The 3-D Lissajous curves for the sPTT and FENE-P models, respectively. (b,d) The respective Cole–Cole plots.
Figure 16 highlights more clearly the temporal evolution of 
$G'_t$ and 
$G''_t$ during the (half) oscillation for the sPTT and FENE-P models for the largest value of 
$Wi$ shown in figure 15. Figures 16(a,c) show the 3-D Lissajous curves for the sPTT and FENE-P models, respectively, whilst figures 16(b,d) show the respective Cole–Cole plots. The colour bar indicates the normalised time between the start point 
$t_0$ (chosen arbitrarily as 
$\gamma = -1$ and 
$\dot {\gamma }=0$) and the end point 
$t_0 +{\rm \pi}$ (i.e. half a period after 
$t_0$). For both the sPTT and FENE-P model responses, the majority of the temporal change in 
$G'_t$ and 
$G''_t$ takes place between 
$t_0$ and 
$(t-t_0)/{\rm \pi} \approx 0.4$, which corresponds to the region directly before and after the stress overshoot in the FENE-P model. In the sPTT model response, initially (i.e. at 
$t_0$) both 
$G'_t$ and 
$G''_t$ are positive, with 
$G'_t$ being slightly larger than 
$G''_t$. Then in the period between 
$t_0$ and 
$(t-t_0)/{\rm \pi} \approx 0.4$, there is first a decrease in 
$G'_t$ with an increase in 
$G''_t$ (which corresponds to strain-softening and shear-thickening), followed by an increase in 
$G'_t$ with a decrease in 
$G''_t$ (which corresponds to strain-stiffening and shear-thinning). For a significant portion of the region between between 
$t_0$ and 
$(t-t_0)/{\rm \pi} \approx 0.4$, 
$G'_t$ is negative, indicating the presence of elastic recoil despite there being practically no self-intersecting loops in the viscous Lissajous curve. In the FENE-P response, both 
$G'_t$ and 
$G''_t$ are smaller at 
$t_0$ than for the sPTT model, and 
$G''_t$ is slightly larger than 
$G'_t$. In the period between 
$t_0$ and 
$(t-t_0)/{\rm \pi} \approx 0.4$, the FENE-P response is drastically different to the sPTT response. Initially, there is a sharp increase in 
$G'_t$ with a decrease in 
$G''_t$. This decrease in 
$G''_t$ is sharp enough that it becomes negative in this region before the stress overshoot. Then, similarly to the sPTT response, a decrease in 
$G'_t$ with an increase in 
$G''_t$ is observed. This behaviour is, however, more extreme for the FENE-P response than for the sPTT response. Then an increase in 
$G'_t$ with a decrease in 
$G''_t$ is observed as the trajectory passes through the sharp stress overshoot. The decrease in 
$G''_t$ in this region for the FENE-P response is so large that 
$G''_t$ is negative after the stress overshoot. In supplementary movies 1 and 2, we show the evolution of the Frenet–Serret frame along the Lissajous curves displayed in figure 16, along with the projections of 
$\boldsymbol {B}$ in the 
$\dot {\gamma }$–
$\tau _{p,12}$ plane (highlighting the sign of 
$G''_t$) and the current position in the Cole–Cole plots.
Figure 16. The SPP analysis for the (a,b) sPTT and (c,d) FENE-P models for 
$De \ (De/a) = 0.5$ and 
$Wi\,\sqrt {\epsilon } \ (Wi/(aL)) = 2$, with 
$\epsilon = 1/L^2 = 1/100$. The colour bar indicates the normalised time in one half of the oscillation from the point indicated by 
$t_0$. (a,c) The 3-D Lissajous curves. (b,d) The Cole–Cole plots.
We now use the toy models outlined in § 3.3 to further highlight the link between the nature of the constitutive model and its behaviour in terms of the SPP analysis. Figures 17(a,b) show the Cole–Cole plots obtained from the SPP framework for the toy FENE-P and toy sPTT models, respectively. Note that 
$De = 0.5$ and 
$Wi = 200$ for the toy sPTT model, whilst 
$De = 1$ and 
$Wi = 100$ for the toy FENE-P model. This corresponds to the conditions for the Lissajous curves shown in figures 9 and 10. For the toy sPTT model, increasing the value of 
$\alpha$ has a seemingly minor effect on the Cole–Cole plots. The general shape remains fairly constant; however, a region of negative 
$G''_t$ develops that corresponds to the developing stress overshoots (seen in figure 9). The link between this negative region of 
$G''_t$ and the stress overshoot is seen clearly in supplementary movie 1 for the FENE-P model response. Essentially, after the stress overshoot, the normal vector 
$\boldsymbol {N}$ points in the positive stress direction as the recoil fades, which points the binormal vector 
$\boldsymbol {B}$ towards the negative strain rate direction. For the toy FENE-P model, the value of 
$b$ has a significant effect on the Cole–Cole plot. For 
$b \rightarrow 1$, the extremities of 
$G'_t$ and 
$G''_t$ become significantly large in magnitude. This indicates that the SPP analysis is highly sensitive to the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship (which relates physically to explicit finite extensibility of the polymer chains) despite the fact that the primary features of the Lissajous curves do not, at least qualitatively, appear to change significantly (see figure 10). More quantitatively, the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship appears to give rise to a region of backflow (i.e. negative 
$G''_t$) before the onset of the stress overshoot. This behaviour cannot be explained easily in terms of the evolution of 
$\boldsymbol{\mathsf{A}}$, as the region of backflow after the stress overshoot can (i.e. 
$Q_r < Q_g$). The region of backflow before the stress overshoot is observed only in the toy FENE-P and toy FENE-CR (see Appendix A) models when 
$b>0$. In supplementary movie 1 for the FENE-P response, it is evident that this region of backflow in the FENE-P response arises due to the exceptionally sharp increase in the stress before the overshoot, which is explained by the fact that the stress grows nonlinearly with 
$\boldsymbol{\mathsf{A}}$ according to (1.13). Figure 18 shows the Cole–Cole plots for the toy sPTT model with 
$\alpha = 1$ and the toy FENE-P model with 
$b = 0$, such that the only difference between the models is the functional form of 
$F(A)$. The qualitative features of the Cole–Cole plots are very similar. The sizes of the deltoids are similar, both responses exhibit backflow after the stress overshoot (owing to the position of 
$F(A)$ in the model and the transient nature of 
$A_{22}$), and neither of the responses exhibits backflow before the stress overshoot (owing to the absence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship). In this sense, the SPP framework can be used to identify quickly the presence of finite extensibility effects in LAOS responses, and differentiate them from other nonlinear effects such as micro-structural destruction.
Figure 17. Cole–Cole plots for (a) the toy FENE-P model at 
$De = 1$, 
$Wi = 100$, 
$L^2 = 100$, and (b) the toy sPTT model at 
$De = 0.5$, 
$Wi = 200$, 
$\epsilon = 1/100$. The corresponding Lissajous curves can be seen in figures 10 and 9.
Figure 18. Cole–Cole plots comparing the toy sPTT model response with 
$\alpha = 1$ and the toy FENE-P model response with 
$b=0$.
3.5. One-dimensional modelling of LAOS
In the previous subsections, it is assumed that the shear rate is uniform in space. For constitutive models that have non-monotonic stress–shear-rate relationships – such as the Giesekus model (Giesekus Reference Giesekus1982) or the Rolie-Poly model (Likhtman & Graham Reference Likhtman and Graham2003) – the material in a shear flow can shear band such that two distinct regions of shear rate exist in the flow. Neither the FENE-P nor sPTT model will shear band in SSSF since their underlying stress–shear-rate curves are both monotonic. However, as we have highlighted already, during LAOS, the FENE-P model behaves in a much more nonlinear manner than during steady shear, and aspects of the response of the model, such as the presence of strong stress overshoots and self-intersecting secondary loops, resemble those of models that do exhibit shear banding. Moreover, recent studies have highlighted that shear banding can occur due to stress overshoots in transient flows, even in models that have monotonic underlying stress–shear-rate constitutive curves (Adams & Olmsted Reference Adams and Olmsted2009; Moorcroft, Cates & Fielding Reference Moorcroft, Cates and Fielding2011; Moorcroft & Fielding Reference Moorcroft and Fielding2013; Carter, Girkin & Fielding Reference Carter, Girkin and Fielding2016). It is therefore sensible to check whether the FENE-P model is capable of shear banding in LAOS.
Using the 1-D modelling approach, and enforcing true creeping flow, the stress response at the top boundary will match the results obtained in the previous subsections if there is no shear banding. However, if shear banding occurs, then we will be able to observe the heterogeneous velocity gradient in the gap, and the stress response at the top wall will not match with the previous results. In order to avoid the known problems of stress selection (Lu, Olmsted & Ball Reference Lu, Olmsted and Ball2000; Olmsted Reference Olmsted2008) and discontinuous strain rates in shear banding, we add a small diffusive term (
$\kappa \,\nabla ^2 \boldsymbol{\mathsf{A}}$) to the right-hand side of the constitutive models during the simulations. The value of 
$\kappa$ was fixed at 
$10^{-9}$. Note that 
$\kappa$ here is dimensionless and is given by 
$\kappa = \lambda D/H^2$, where 
$D$ is the diffusion coefficient in m
$^2$ s
$^{-1}$. In the supplementary material, we show that this methodology is capable of capturing shear banding in LAOS using the Rolie-Poly (ROuse LInear Entangled POLYmers) model.
Figures 19 and 20 show the viscous Lissajous–Bowditch plots for the sPTT and FENE-P with 
$L^2 1/\epsilon = 100$ models, respectively, where the stress response has been computed with the 0-D approximation (black dashed lines) and 1-D simulation (solid lines). For both models, both the 0-D and 1-D methods give practically identical responses, indicating that the 0-D approximation is adequate for describing the model responses properly, and that no shear banding is occurring in the 1-D simulations. Since the methodology used for the 1-D simulations is capable of predicting shear banding, it is reasonable to assume that the FENE-P model does not shear band in LAOS, despite the fact that the model response is significantly more nonlinear in LAOS than it is in SSSF, at least in the ranges of 
$De/a$ and 
$Wi/(aL)$ investigated.
Figure 19. Lissajous–Bowditch plots (viscous projection) for the sPTT model. Blue solid lines show the 0-D approximation solution from § 3.1; black dashed lines show the results from the 1-D simulations, where the stress is computed at the top (moving) boundary.
Figure 20. Lissajous–Bowditch plots (viscous projection) for the FENE-P model (
$L^2 = 100$). Red solid lines show the 0-D approximation solution from § 3.1; black dashed lines show the results from the 1-D MOL simulations, where the stress is computed at the top (moving) boundary.
4. Conclusions
We have compared the response of the sPTT and FENE-P constitutive models in LAOS when the parameters in the models are chosen such that the models are mathematically identical for steady and homogeneous flows. The results show that the FENE-P model behaves in a significantly more nonlinear manner than the sPTT model in LAOS, as it does in other transient flows such as start-up flow. The FENE-P model exhibits strong stress overshoots and large self-intersecting secondary loops in the viscous Lissajous curves, whereas the sPTT model does not. The stress overshoots and self-intersecting secondary loops arise from the FENE-P model due primarily to the extensibility function being multiplied with 
$\boldsymbol{\mathsf{A}}$ instead of 
$(\boldsymbol{\mathsf{A}}-\boldsymbol{\mathsf{I}})$ in the conformation tensor form of the constitutive model, which can be thought of in terms of network theory as the system exhibiting faster rates of micro-structural destruction than creation. This drives a change in 
$A_{22}$ and 
$A_{33}$ during the oscillation that, due to the use of the upper-convected derivative, causes the elastic recoil of 
$A_{12}$ to exceed the growth of 
$A_{12}$ for increasing rates of deformation, which ultimately leads to the pronounced stress overshoots. It is important to note that the differences between these two models arise for both Eulerian unsteady flows (such as LAOS) and Lagrangian unsteady (or inhomogeneous) flows, so the differences between the model responses such as the stress overshoots will also occur in Eulerian steady flows that are Lagrangian unsteady due to, perhaps, expansions and contractions in a complex geometrical domain. Such complex geometries are often encountered in many industrial flows and processes. In fluid regions with strong accelerations, one can expect much sharper changes in the polymeric stress with the FENE-P model than the sPTT model.
Although it has been shown previously analytically for the FENE-P and sPTT models that the stress–strain-rate curves scale with 
$Wi/(aL)$ and 
$Wi\,\sqrt {\epsilon }$, respectively, for simple steady shear, we have been able to show with our numerical results that the sPTT LAOS response also scales with 
$Wi\,\sqrt {\epsilon }$, but the FENE-P response in LAOS scales only with 
$Wi/(aL)$ at large enough values of 
$L^2$. We have also been able to explain this analytically. This is shown to be caused by the specific functional form of the extensibility function as well as its position in the FENE-P constitutive model (inside the brackets in the recoil term rather than on the outside).
Using the sequence of physical processes framework, we highlight the differences in the model responses in terms of the transient moduli 
$G'_t$ and 
$G''_t$. The Cole–Cole plots show that the FENE-P model exhibits significantly more complex rheological behaviour during the oscillation. A key result obtained from the SPP analysis is that the FENE models (both FENE-P and FENE-CR) exhibit backflow (i.e. negative 
$G''_t$) both before and after the stress overshoot. The region of backflow before the stress overshoot is shown to be linked to the presence of the extensibility function in the stress–conformation-tensor relationship. This highlights that the SPP framework can be particularly useful for identifying the correct form of constitutive models for viscoelastic materials from LAOS data.
Although the FENE-P model in LAOS behaves more like models that exhibit shear banding, such as Giesekus and Rolie-Poly, than it does in steady shear, the FENE-P model does not appear to be capable of shear banding in LAOS. This was investigated by solving both the momentum equation and the constitutive model in a 1-D gap of fluid using the method of lines technique. For all cases with the sPTT and FENE-P models, the velocity gradient remained constant across the gap, and the Lissajous–Bowditch plots were almost identical when the 1-D and 0-D solutions were compared.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2023.977.
Acknowledgements
The authors acknowledge the developers of the SPP software for kindly providing us with a copy of the software. We would also like to thank all of the reviewers for their feedback and comments. We thank Reviewer 3 specifically for showing us how to prove the universal scaling of the sPTT model LAOS response in § 3.1.
Funding
The authors acknowledge the financial support of the Center in Advanced Fluid Engineering for Digital Manufacturing, UK (CAFE4DM) project (grant no. EP/R00482X/1)
Declaration of interests
The authors report no conflict of interest.
Appendix A. Comparing the FENE-CR and sPTT models
In § 3.1, we highlight that the FENE-CR model scales universally with 
$Wi/L$ only for 
$L^2 \gg 3$. We also point out that the evolutions of 
$\boldsymbol{\mathsf{A}}$ for the sPTT and FENE-CR models become identical for 
$L^2 \gg 3$ and in the case that terms of 
$O(x^2)$ can be neglected in the expansion of 
$F(A)_{{FP}}$, where 
$x = (A_{11}-1)/L^2$, which corresponds to the MAOS regime. This is also highlighted in figure 21. We introduce here a toy FENE-CR model that, under LAOS flow, is defined as
$$\begin{gather} De\,\frac{\partial}{\partial t} \boldsymbol{\mathsf{A}} - Wi\,(\boldsymbol{\mathsf{A}} \boldsymbol{\cdot}\boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\mathrm{T}} \boldsymbol{\cdot} \boldsymbol{\mathsf{A}} - \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{\mathsf{A}}) =- F(A)_{{FP}}\, (\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}} ), \end{gather}$$
$$\begin{gather}\boldsymbol{\tau\!}_p = \frac{1-\beta}{Wi}\,(F(A)_{{FP}})^b(\boldsymbol{\mathsf{A}} - \boldsymbol{\mathsf{I}}), \end{gather}$$
where 
$0 \leq b \leq 1$. We remind the reader that 
$F(A)_{{FP}} = L^2/(L^2 - \mathrm {tr}(\boldsymbol{\mathsf{A}}))$. When 
$b = 0$, the only difference between the toy FENE-CR and sPTT models is the specific form of 
$F(A)$ used (i.e. 
$1+\epsilon (\mathrm {tr}(\boldsymbol{\mathsf{A}})-3)$ versus 
$L^2/(L^2 - \mathrm {tr}(\boldsymbol{\mathsf{A}}))$). In this case, the stress response of the toy FENE-CR model will now also become equivalent to that of the sPTT model in the MAOS regime. We show this in figure 22. For 
$Wi\,\sqrt {\epsilon } \ (Wi/L) = 0.2$, the responses of both models are nonlinear, due to the fact that 
$F(A)$ becomes transient and larger than unity during the oscillation. However, since 
$x$ is small during the oscillation and 
$L^2 \gg 3$, 
$F(A)_{{FP}}$ is essentially linear and equivalent to 
$F(A)_{{sPTT}}$, so the stress responses of both models are practically identical. For larger values of 
$Wi\,\sqrt {\epsilon } \ (Wi/L)$, 
$F(A)_{{FP}}$ starts to deviate from 
$F(A)_{{sPTT}}$ at points in the oscillation due to the increased values of 
$A_{11}$, so the stress responses begin to differ. As is likely expected, the toy FENE-CR model response appears more nonlinear than the sPTT response, owing to the increased nonlinearity in 
$F(A)$. The responses are, however, at least qualitatively similar even for moderate values of 
$Wi\, \sqrt {\epsilon } \ (Wi/L)$.
Figure 21. Extensibility function 
$F(A)$ versus 
$x$, where 
$x = \epsilon (A_{11}-1)$ for the sPTT model, and 
$x = (A_{11}-1)/L^2$ for the FENE-CR model. The black dashed line shows the expression for 
$F(A)_{{FP}}$ in the case that 
$L^2 \gg 3$. Yellow to red lines show 
$F(A)$ for the FENE-CR model, where 
$A_{22} = A_{33} = 1$, so 
$F(A) = (1-x-3/L^2)^{-1}$.
Figure 22. Viscous Lissajous curves (
$\tau _{p,12}$ versus 
$\dot {\gamma }$) for the sPTT model (toy FENE-CR model) for 
$De = 0.2$ and varying 
$Wi\,\sqrt {\epsilon } \ (Wi/L)$. Insets show the ranges of 
$F(A)$ for each model during the oscillation. Here, 
$b = 0$ and 
$L^2 = 10^5$ for the toy FENE-CR model.
Figure 23 shows the viscous Lissajous curves for the toy FENE-CR model with varying values of 
$b$ at 
$De = 1$, 
$Wi = 100$ and 
$L^2 = 100$. For all values of 
$b$ in the toy FENE-CR model, despite the fact that 
$Q_g$ is linear, a small stress overshoot is observed due to the nonlinearity of 
$F(A)_{{FP}}$, and hence nonlinearity of 
$Q_r$. This highlights again that the pronounced stress overshoots, which are observed in the FENE-P, toy FENE-P and toy sPTT model responses, are closely linked with the nonlinearity of 
$Q_g$ and thus the transient nature of 
$A_{22}$. This becomes particularly clear when figures 23(c) and 10(f) are compared. Note when comparing these two figures that 
$F(A)$ is moving from the inside to the outside of the brackets in the recoil term, and 
$F(A)$ does not appear in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship in either case.
Figure 23. Viscous Lissajous curves (
$\tau _{p,12}$ versus 
$\dot {\gamma }$) for the toy FENE-CR model with (a) 
$b=1$, (b) 
$b=0.5$, and (c) 
$b=0$. Note that small stress overshoots are observed in all cases.
Figure 24 shows the Cole–Cole plots for the toy FENE-CR model responses displayed in figure 23. For 
$b > 0$, a region of backflow is observed before the small stress overshoot, which is also observed in the toy FENE-P model response for 
$b > 0$. This further highlights that the SPP framework can identify clearly the presence of 
$F(A)$ in the 
$\boldsymbol {\tau\!}_p$–
$\boldsymbol{\mathsf{A}}$ relationship, and thus help to identify a suitable form of constitutive model from LAOS data.
Figure 24. Cole–Cole plots (
$G''_t$ versus 
$G'_t$) for the toy FENE-CR model with (a) 
$b=1$, (b) 
$b=0.5$, and (c) 
$b=0$. Note that for 
$b=0$, only one region of viscous backflow (i.e. 
$G''_t < 0$) is observed, whereas two regions of backflow are observed for 
$b=0.5$ and 
$b=1$.
Open access
