A comparison between the FENE-P and sPTT constitutive models in Large Amplitude Oscillatory Shear (LAOS)

The FENE-P and sPTT viscoelastic models are widely used for modelling of complex fluids. Although they are derived from distinct micro-structural theories, these models can become mathematically identical in steady and homogeneous flows with a particular choice of the value of the model parameters. However, even with this choice of parameter values, the model responses are known to differ from each other in transient flows. In this work, we investigate the responses of the FENE-P and sPTT constitutive models in Large Amplitude Oscillatory Shear (LAOS). In steady-shear, the shear stress scales with the non-dimensional group $Wi/(aL) \ (Wi\sqrt{\epsilon})$ for the FENE-P (sPTT) model, where $Wi$ is the Weissenberg number, $L^2$ is the limit of extensibility in the FENE-P model ($a$ being $L^2/(L^2-3)$) and $\epsilon$ is the extensibility parameter in the sPTT model. Our numerical and analytical results show that, in LAOS, the FENE-P model only shows this universality for large values of $L^2$ whereas the sPTT model shows this universality for all values of$\epsilon$. In the strongly non-linear region, there is a drastic difference between the responses of the two models, with the FENE-P model exhibiting strong shear stress overshoots which manifest as self-intersecting secondary loops in the viscous Lissajous curves. We quantify the non-linearity exhibited by each constitutive model using the Sequence of Physical Processes framework. Despite the high degree of non-linearity exhibited by the FENE-P model, we also show using fully non-linear 1D simulations that it does not shear band in LAOS within the range of conditions studied.


Introduction
Viscoelastic fluids are ubiquitous in many industrial sectors including fast moving consumer goods, food, and health-care, among others, and it is of significant importance that we are able to correctly model these flows in a wide range of geometries.The Oldroyd-B model (Oldroyd 1950) is given as where   is the polymeric stress,  is the viscoelastic relaxation time,   is the polymeric viscosity and D is the rate-of-strain tensor given by D = 1/2(∇ + ∇ T ).
▽   denotes the upper-convected time derivative of the polymeric stress tensor, which is given as ▽   = D  /D −   • ∇ − ∇ T •   .For many viscoelastic models, including the Oldroyd-B model, the total stress  appearing in the momentum equation is related to the polymeric stress by  =  − I =   + 2  D − I, where  is the extra-stress tensor,  is the pressure and I is the identity tensor.The solvent contribution to the stress is expressed as a Newtonian fluid with a viscosity   .
Whilst the simplicity of the Oldroyd-B model makes it particularly useful for solving problems analytically (Boyko & Stone 2022;Ghosh et al. 2021;Norouzi et al. 2018;Qi & Xu 2007;Rajagopal & Bhatnagar 1995;Zhao et al. 2013) and testing and validating computational codes (Duarte et al. 2008;Habla et al. 2014;Mompean & Deville 1997), it has a number of wellknown short-comings.Likely the most well-known short-coming 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 et al. 1987).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 FENE-P (Finitely Extensible Non-linear Elastic with Peterlin closure) (Bird et al. 1980) and the sPTT (simplified Phan-Thien-Tanner) (Phan-Thien & Tanner 1977) models.
The original FENE model (Warner 1972) 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 each polymer molecule.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. 1980;Keunings 1997).The springs are finitely extensible since the elastic stress increases non-linearly 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 inter-connected 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 simultaneously created and destroyed, 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  1 =  11 −  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 thickening 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  2 (1/) for the FENE-P (sPTT) model, where  2 and  represent the respective extensibility parameters in the FENE-P and sPTT models.
The FENE-P constitutive model is given in stress tensor form as where   ≡ tr(  ), or equivalently   .It should also be noted that the sPTT model (Equation (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 1978) or even generalised (Ferrás et al. 2019) terms, which are believed to help the model perform better under strong deformations.In this study, we will only use the sPTT model with the linear function (Equation (1.5)), we shall 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 Equations (1.4) and (1.5), it is observed that with the parameter substitutions  = 1/ 2 and  sPTT =  FP /, the FENE-P and sPTT-Linear models become mathematically identical for steady and homogeneous flows.The equivalence of these two models was first noted in the study of Cruz et al. (2005), who derived analytical solutions for fully-developed pipe and channel flows with the FENE-P and sPTT models.Latreche et al. (2021) 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. (2022) then investigated the FENE-P and sPTT models for a range of steady and homogeneous flows, as well as unsteady and in-homogeneous 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. (2022) was the cross-slot.For viscoelastic flows in the cross-slot, the elastic stresses causes a symmetry-breaking instability to occur at a critical , which has previously been well studied and characterised (Afonso et al. 2010;Cruz et al. 2014;Davoodi et al. 2019Davoodi et al. , 2021;;Haward et al. 2012;Poole et al. 2007;Rocha et al. 2009;Xi & Graham 2009).Davoodi et al. (2022) observed that that the critical value of  for the onset of the asymmetry is lower for the FENE-P model than for the sPTT model when relatively low (high) values of  2 () 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 which might induce Lagrangian unsteadiness, even for an Eulerian steady flow.Varchanis et al. (2022) recently highlighted that even the Oldroyd-B model exhibits complex rheological behaviour in Lagrangian transient flows, which has significant consequences on, for example, the understanding of pressure drop measurements across a contraction.It is therefore of significant importance to compare and understand how these non-linear models behave in transient flows.
An ideal way of probing the transient non-linear response of viscoelastic materials and models, and in particular classifying complex fluids (Hyun et al. 2002), is with Large Amplitude Oscillatory Shear (LAOS), which has become a widely used technique for characterising non-linear viscoelasticity experimentally (Hyun et al. 2011;Leblanc 2008;Sun et al. 2011;Szopinski & Luinstra 2016), theoretically (Bae & Cho 2017;Gurnon & Wagner 2012;Kammer & Castañeda 2020;Khair 2016), and numerically (Cordasco & Bagchi 2016;D'Avino et al. 2013;Ewoldt & McKinley 2010).In Small Amplitude Oscillatory Shear (SAOS) the shear stress response of a material or constitutive model is approximately linear and given by  ,12 =  0 [ ′ sin() +  ′′ cos()] where  0 and  are the amplitude and angular frequency of the oscillation respectively. ′ and  ′′ 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,  is very often estimated as the inverse of the frequency at which  ′ and  ′′ cross over in a frequency sweep.However, as  0 increases, flowinduced micro-structural changes take place during the oscillation (Gilbert & Giacomin 2016), 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 accurately reconstructed using a single mode of  ′ and  ′′ .Multiple frameworks have been developed for quantitative analysis of the non-linear stress response obtained from LAOS, namely these are Fourier Tranform Rheology (Wilhelm et al. 1998), Stress Decomposition (Cho et al. 2005) with Chebyshev analysis (Ewoldt et al. 2008), and a Sequence of Physical Processes (Rogers et al. 2011).LAOS is considered to be especially useful for the purpose of fitting constitutive models to experimental data (Bae & Cho 2015).Calin et al. (2010) used LAOS to fit the spectrum of the tensorial mobility parameter of the Giesekus model (Giesekus 1982) to experimental data using an iterative numerical solution.Gurnon & Wagner (2012) then derived an asymptotic solution for the Giesekus model in oscillatory shear, which they use to easily fit 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. 2014), the co-rotational Maxwell Model (Giacomin et al. 2015), and the White-Metzner model (Merger et al. 2016), among others.Hyun et al. (2007) 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. (2011) and Kamkar et al. (2022).
For a purely oscillatory shear flow where the strain rate is uniform in space, the strain () and strain rate () are given by () =  0 sin() and () =  0  cos(), respectively.We define here the non-dimensional polymeric stress  *  =   /( 0 [  +   ]), the non-dimensional velocity gradient (∇) * = ∇/( 0 ), and the non-dimensional time  * = .We also define the Weissenberg number as  =  0  and the Deborah number as  = .As usual, we define the dimensionless parameter  as the ratio of the solvent viscosity to the total viscosity (polymeric viscosity plus solvent viscosity), such that  =   /(  +   ).Using these definitions to nondimensionalise the FENE-P (Equation (1.3)) and sPTT (Equation (1.5)) models, and dropping the asterisks upon non-dimensionalisation, we have respectively where For all of the models discussed in this study, including the Oldroyd-B model, the extra-stress tensor is given as  =   + 2D.In dimensionless form, it is upon substitution of  = 1/ 2 and  sPTT =  FP / that the FENE-P and sPTT models become mathematically identical for steady ( = 0) and homogeneous ( With regards to the definitions of  and  for LAOS, Kamani et al. (2023) recently highlighted that using a time-independent value of  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 , 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  and  naturally appear 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  and  would naturally appear from the equations.In the FENE-P and sPTT models, one might also think of an "effective" relaxation time based on  and  (  ).Therefore, we note that, depending on the model or material in question, one may start to question the correct choice of definition for  and  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 only use the time-independent values for  and  defined previously.
In the limit that  (/) → 0 and  (/) → 0, the sPTT (FENE-P) models, as well as the Oldroyd-B model, reduce to that of a Newtonian fluid.Note that lim →0 (1/ (  ) FP ) = 1/ and D(1/)/D = 0.In the limit that  → 0, the response of each model reduces to its respective steady state response.In the case that  (  ) sPTT → 1 for the sPTT model, or  (  ) FP / → 1 for the FENE-P model, the Oldroyd-B model is obtained.Note that for the FENE-P model  (  ) FP / → 1 is equivalent to  (  ) FP →  and therefore / (1/ (  ) FP ) → 0 and  • ∇(1/ (  ) FP ) → 0 and so the last term on the right hand side of Equation (1.7) vanishes in this limit.
Viscoelastic constitutive models can also be written for the conformation tensor, A. For dumbbell models such as the FENE-P model, A can be given as A = ⟨⟩/ 2  where  is the end-to-end vector of an individual dumbbell (the angled brackets represents the ensemble average) and  2  is the square of the magnitude at equilibrium given as  2  = ⟨ • ⟩  /3 (Alves et al. 2021).In general, for other viscoelastic models,  might represent the end-to-end vector of polymer chains or sub-chains (Hoyle & Fielding 2016), rather than the dumbbell vector specifically.The dimensionless FENE-P model is given in conformation tensor form as follows where which can also be re-written as The sPTT model is given in conformation tensor form as where  ( ) = [1 +  (tr(A) − 3)] (1.12)Note that  ≡ tr(A).  is then recovered from the solutions of Equations (1.10) to (1.12) as follows ( ) is therefore also be defined on a per-model basis as As highlighted by Davoodi et al. (2022), the evolution equation for A in network theory models follows a general form given by where  ( ) and  ( ) represent, respectively, the rates of destruction and creation of microstructures.For the sPTT model,  ( ) =  ( ) =  ( ) sPTT .It is therefore observed, given Equation (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. 2022).
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. 2002).Physically, each classification is believed to correspond to a particular type of underlying micro-structural interaction.Sim et al. (2003) 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 (2018) 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  = 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  ′ .With increased oscillation frequency, the FENE-P response changed to a type III response where  ′′ 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 of standard FENE dumbbells and FENE-type networks in LAOS using a technique known as the Brownian Configuration Field Method (BCFM) (Gómez-López et al. 2019;Vargas et al. 2023).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 microstructures.As already mentioned, in the context of the PTT model framework (Equation (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 2010).Ng et al. (2011) 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 which is essentially a blend between  ( ) sPTT at low stretching and  ( ) FP at high stretching.They also include  ( ) FP in the   -A relationship and so the constitutive model represents a FENE-type network model.The model was able to at least qualitatively predict 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  ( ) FP ) and they introduce a modified function, which diverges to infinity before the FENE limit is approached to empirically temper the magnitude of the stress overshoots.Keunings (1997) 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 (2020), who showed that with increasing  and  the sPTT response clearly deviates away from the linear UCM/Oldroyd-B response.However, in this study, no quantitative analysis of the generated waveforms was conducted.
Despite the fact 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.

0D modelling
The majority of the results in this study are obtained by solving constitutive equations assuming that A,   , and  are uniform in space.We denote this approach to solving the equations as the 0D method.This methodology will now be detailed.
For an ideal oscillatory shear flow, the dimensionless constitutive model can be solved using Since the FENE-P model in stress tensor form (Equation (1.3)) cannot be easily expressed 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 only work with dimensionless variables and we re-confirm that the asterisks denoting the dimensionless variables have been dropped for brevity.
The following system of ODEs are obtained for the time evolution of A according to the FENE-P model and according to the sPTT model where   is recovered from A with Equation (1.13).With an initial condition of A = I (ie.  = 0), d 22 /d = d 33 /d = 0 at all times for the sPTT model, and so  22 and  33 remain fixed at unity.However, for the FENE-P model, under large oscillatory deformations  22 and  33 will become time-dependent and lower than unity, although it is still the case that  22 =  33 .We note here that, for the FENE-P model,  22 =  33 = / ( ) FP in steady-state conditions, and so  ,22 =  ,33 = 0 according to Equation (1.13). 22 and  33 therefore 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  ,22 and  ,33 are also nonzero.For all 0D simulations, we omit the solvent contribution to the stress by setting  = 0 so that we only study the response of the viscoelastic constitutive model itself.Therefore, the Oldroyd-B model is here-on-in denoted as the Upper-Convected Maxwell (UCM) model.For the FENE-P and sPTT models, we performed simulations for five values of  2 = 1/ = [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 Section 3, data is only plotted for the final oscillation when the system is steady-periodic (i.e. the limit cycle).

1D modelling
We also use a 1D modelling approach by solving both the momentum equation and the constitutive model in a 1D 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 1D approach, we use the Method Of Lines (MOL) technique, in which spacial 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 -direction.The first and second order spacial derivatives of a scalar variable  are discretised using a fourth-order finite difference scheme, respectively, as where the index  denotes the node number in a uniformly discretised domain with   elements.Δ is the distance between neighbouring cells, given as Δ =   −  −1 .The (dimensional) velocity at the top boundary    () is varied according to    =  0 cos() where  is the gap height. 0  represents the strain rate amplitude.For non-dimensionalisation,  is used for the length scale,  0  is used for the velocity scale, and the time is still non-dimensionalised with . and  are then defined as they are for the 0D approach.
Assuming that the only nonzero velocity component is in the -direction, and the flow is uniform in the -direction, the resulting system of Partial Differential Equations (PDEs) to be solved can be expressed in dimensionless form as where the Reynolds number Re is defined as Re =  2  0 /(  +  ), and the tensor 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 Equations (2.4) and (2.5), which turns the system of PDEs into a system of ODEs.For the momentum equation, ,12 is computed from A using Equation (1.13), and then its gradient is discretised with Equation (2.4).For the 1D MOL modelling, we could not totally omit the solvent contribution to the stress due to stability issues.We therefore used  = 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  in the range of  and  investigated.This is shown in the Supplementary Material.We also enforce true creeping flow so that inertia is neglected (i.e. the left hand side of Equation (2.6a) is zero).
For the spacial resolution, we used   = 128 , which proved sufficiently accurate to ensure that the results were independent of   .This is also shown in the Supplementary Material.At the bottom wall the velocity was fixed at zero.The components A were linearly extrapolated at the top and bottom boundaries.Simulations were initiated with  = 0 and A = I.We integrated the resulting system of equations in MATLAB using the adaptive-step ODE solver ode15s, which can solve systems of DAEs using the Mass Matrix approach.We simulated the flow until the limit cycle was reached.Again, only data for the limit cycle is presented in Section 3.

Results and discussion
For all of the results except those presented in Section 3.5, the 0D method is used (with  = 0) to obtain the solutions.In Sections 3.1 and 3.2, we investigate the model responses in LAOS with the parameter substitutions  = 1/ 2 ,  UCM =  sPTT =  FP /, and  UCM =  sPTT =  FP /.We present the results for various values of / () and / () for the FENE-P (sPTT or UCM) model.In Section 3.3, we investigate the responses of "toy" models to help explain the observations from Sections 3.1 and 3.2.In Section 3.4, we analyse the model responses using the Sequence of Physical Processes methodology.Finally, in Section 3.5 we use the 1D 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  ,12 versus  and plots of  ,12 versus .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  and ) space.For a more detailed overview of Lissajous-Bowditch plots, the reader is referred to the review by Hyun et al. (2011).

Scaling of the Lissajous curves
Figures 1 and 2 show, respectively, the viscous and elastic projections of the Lissajous-Bowditch plots in the / () − / () 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  2 = 1/ for the FENE-P and sPTT models.The FENE-P and sPTT models deviate from the UCM model at high / () and particularly at low (high) values of  2 ().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  2 (), which matches the observations of Davoodi et al. (2022) for start-up shear flow.
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 (Equation (1.7)) yields the following solution for the (polymeric) shear stress For constant , the solution for  ,12 in SSSF then evidently depends on the parameter /().
Oliveira & Pinho (1999) 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  √ .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 /() ( √ ) for the FENE-P (sPTT) models in SSSF was also shown by Oliveira et al. (2004) and Latreche et al. (2021).A recent study by Yamani & McKinley (2022) showed, analytically, that the SSSF response of the FENE-P model scales with the dimensionless parameter /, which differs from the parameter used in this study, /().However, in the version of FENE-P model used in the study of Yamani & McKinley (2022) the value of  was  ,12 vs   causes the responses for the various values of  to become universal (for constant ).Whilst this is expected in SSSF due to the form of the analytical solution (Equation (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  22 =  33 = 1 due to the fact  ( ) is on the outside of the (A − I) term in the constitutive model) and introducing new variables  =  (  11 − 1) and  =  12 /, which gives  ,12 vs Thus, for constant values of  (in our case  = 0), the system only depends on the parameters  and  √ .The FENE-P LAOS response does not scale universally with /() under uniform oscillation, as its SSSF response does.For / < 0.1 there is practically no difference in the curves for the various values of  2 since the flow is approaching SSSF.For / ⩾ 0.1 and /() ⩾ 1, the difference in the FENE-P response with varying  2 becomes significant.However, the FENE-P response does seem to become universal at constant values of / and /() for large enough values of  2 .This is highlighted in Figure 5 by the fact that the responses for  2 = 100 and  2 = 1000 are practically identical.There are at least two potential reasons that the FENE-P response is not universal for constant values of / and /() at low values of  2 .One is that the functional form of  ( ) in the FENE-P model is different to that in the sPTT model (see Equation (1.14)).Another is that the position of  ( ) in the conformation tensor form of  ,12 vs  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 for 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 Introducing the new variable  = (  11 − 1)/ 2 , the extensibility function can be rewritten as The system of equations for the FENE-CR model for  2 ≫ 3 can then be rewritten, using also  =  12 /, as Therefore, the FENE-CR model has universal solutions for constant values of  and / only in the case that  2 ≫ 3, and the breakdown of the universality can be caused solely due to a change in the functional form of the extensibility function  ( ) 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  ( ) and its position in the model (i.e. ( ) (A − I) versus ( ( )A − I)) are different to those in the sPTT model.We will explore the effect of the position of  ( ) in the constitutive model on the scaling of the response further in Section 3.3.We also note here that (1 − ) −1 can be expanded as 1 +  + O ( 2 ).Therefore, for  = 1/ 2 , the evolution of A for the FENE-CR model becomes mathematically identical to that for the sPTT model in the MAOS regime in the case that  2 ≫ 3.In the MAOS regime, where the response is weakly non-linear, terms of O ( 2 ) can be neglected in the expansion of  ( ).There is, however, a difference in the stress response due to the presence of the extensibility function in the   -A relationship in the FENE-CR model.We highlight this further in Appendix A.

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   with  = 1/ 2 and  sPTT =  FP /, 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 / () → 0, which corresponds to the system approaching SSSF.Both models also exhibit identical responses for /() ( √ ) → 0, however in this case both models reduce to the UCM model.For large values of both / () and /() ( √ ), 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 non-linear viscoelastic models, such as those derived specifically for wormlike micelles.These can be observed in particular in the plots for / = 1 and /() = 10, which are  ( ) vs  shown at a larger scale in Figure 5a.Similar stress overshoots were also observed in the FENE-P model response during start-up shear flow (Davoodi et al. 2022).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  = 0 (or  = 1), indicating that elasticity is being relieved through recoil faster than it is being accumulated through increased rates of deformation (Ewoldt & McKinley 2010).These secondary loops are associated with strongly non-linear 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  ( ) given by Note that with the relevant expression relating   to A for each model (Equation (1.13)) it is the case that  ( ) FP =  (  ) FP and  ( ) sPTT =  (  ) sPTT .We point this out just to clarify that as  ( ) → 0, both the conformation tensor and stress tensor forms of the models asymptote towards the UCM model.
Figure 6 shows  ( ) against the dimensionless strain rate for the corresponding Lissajous curves presented in Figure 3.For both models,  ( ) is approximately 0 for values of /() ( √ ) ⩽ 0.1 at any value of / ().Consequently, the dimensionless shear stress response shown in the corresponding Lissajous curves is essentially that of the UCM model.For / () = 0.01, even when  ( ) increases and the model responses become increasingly non-linear, for each model, the solution is universal for various values of  2 and  at constant /() ( √ ), which has been explained already by the fact that the steady-shear response of the dimensionless model contains only the scaling parameter /() ( √ ).Note here that the response of  ( ) is also the same for both the sPTT and FENE-P models, despite the fact that the functions  ( ) FP / and  ( ) sPTT are not explicitly equivalent for  = 1/ 2 ,  sPTT =  FP /, and  sPTT =  FP /.Therefore, in SSSF, i.e.,  → 0, it is the case that tr(A) sPTT = ( ( ) FP /) • tr(A) FP , which is also implied from Equations (1.13) if (  ) FP = (  ) sPTT .Thus, in steady and homogeneous flows, tr(A) for each model differs by a factor of  ( ) FP /, highlighting the difference in the physical interpretation of the polymeric stress from A in each model.For the higher values of / () and /() ( √ ), it is observed again that the sPTT solution for  ( ) is universal for constant  and

√
with varying , but the FENE-P solution for  ( ) is only universal for constant / and /() at large values of  2 .It is also observed that there is significant correspondence between the results in Figures 6 and 3. Notably, the overshoots in  ( ) appear to correspond at least qualitatively with the overshoots in   .This will be discussed in more detail next.
For / and /() = 10, we present in Figure 7 3D plots of  ( ),  12 ,  22 and  ,12 against  and  for the FENE-P model response in LAOS (with  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 4 regions (note that only one half of the curve is shown in terms of the range of ).It is important here to refer back to Equations (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).
In Region I,  12 is negative and the rate of deformation is increasing, meaning that both the growth of  12 due to deformation and the elastic recoil are acting in the same direction, causing a positive rate of change in time of  12 .The rate of change of  12 is governed by a balance between growth due to deformation and reduction due to elastic recoil.In the sPTT model,  22 = 1 and so the growth of  12 due to increasing rates of deformation is proportional to the strain rate.In the FENE-P model however,  22 is time-dependent for large values of /() and so the growth of  12 due to increasing deformation rates is non-linear (this will be discussed in more detail in Section 3.3).Whilst  ( ) is decreasing in Region I, which reduces the degree of the elastic recoil, the growth of  22 is relatively large and so ultimately  12 grows non-linearly in Region I. Despite the fact  12 grows in Region I,  ,12 grows only slightly (and seemingly linearly) due to the presence of  ( ) in the   -A relationship and the fact that  ( ) is decreasing.
Region II is characterised by a sharp increase in  ( ) as tr(A) →  2 , which causes  ,12 to increase rapidly but also limits somewhat the growth of  12 since the elastic recoil increases. 22 also goes through a maxima in Region II and begins to decrease, which causes a reduction in the growth of  12 with increasing rates of deformation, exacerbating the recoil effect.In Region III,  ( ) is large enough that the elastic recoil now exceeds the growth of  12 due to the increasing rate of deformation, which drives a negative rate of change of  12 for increasing strain rates.In this region,  ( ) also decreases as tr(A) decreases, and so  ,12 decreases rapidly due to both decreasing  12 and  ( ) (and the fact  ( ) is present in the   -A relationship).
At the start of Region IV,  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  ( ), which drives a reduction in the shear stress.However, as the deformation rate decreases, there is no significant decrease in  12 due to the fact that  22 is increasing and  ( ) is decreasing, which acts ultimately to slow the recoil of  12 when the rate of deformation is decreased.
Considering the previous analysis, we then present 3D plots of  ( ),  12 ,  22 and  ,12 against  and  for the FENE-P ( 2 = 100) and sPTT model response for / () = 1 and /() ( √ ) = 10 in Figure 8.Note that the sharp overshoots in  ( ) are not observed for the sPTT model response and also note the significant differences in the  22 responses of each model. 12 is generally much larger for the sPTT model than for the FENE-P model, but the values of  ,12 are fairly similar for large parts of the oscillation due to Equation (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 (Equations (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 which at first might seem steady-state due to the Eulerian steadiness, but might be Lagrangian unsteady (or inhomogeneous).

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 only employ the 0D modelling approach to obtain the solutions, and we also use again  = 0 for all results.
For the toy sPTT model, we start with the generic network model given in Equation (1.15) and adjust the rates of micro-structural destruction  ( ) and creation  ( ) as Thus, for  = 0.5 we recover the sPTT model (with the value for  halved), and for  = 1 we recover a model with a similar form to the FENE-P model but without  ( ) in the   -A relationship. therefore, in a way, controls the position of  ( ) in the recoil term of the constitutive model.The LAOS response for the toy sPTT model with  = 1/100 at  = 0.5 and  = 200 is presented for various values of  in Figure 9.As  is increased, the stress overshoots and self-intersecting secondary loops become significantly more pronounced, as is expected for systems which exhibit large rates of micro-structural destruction (Davoodi et al. 2022;Vargas et al. 2023), 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 A instead of (A − I), and not due to the fact that the extensibility function appears in the   -A relationship, or due the difference in the natures of  ( ) FP and  ( ) sPTT .We also show this more explicitly using a toy FENE-P model where the evolution equation for A is unchanged (given by Equation (1.11)), but   is given by where 0 ⩽  ⩽ 1.When  = 1, the original FENE-P model is obtained, and when  = 0, the   -A relationship reverts to the original form given by the Kramers' relation (Kramers 1944), and that used for the sPTT model.The viscous Lissajous curves for the toy FENE-P model with  2 = 100 at  = 1 and  = 100 with varying values of  are shown in Figure 10.The stress overshoots and self-intersecting loops are observed for all values of , indicating that the presence of  ( ) in the   -A relationship does not explicitly cause pronounced shear stress overshoots in the transient response.As  → 0, the region directly after the stress overshoot becomes significantly flatter, and the Lissajous curves resemble more those of constitutive models derived for wormlike micellar systems.We also note for both Figures 10 and 9 that the scale of the  axis in each plot is not fixed.The parameters  and  do affect significantly the maximum stresses reached in the oscillation.In this sense, the stress overshoots that we are discussing here are relative.
In the standard form of the sPTT model ( = 0.5), the growth term for  12 due to the rate of deformation is equal to (/)cos() since  22 = 1.As  → 1, there is significant deviation of  22 from unity during the oscillation, and hence deviation of the growth term from a pure cosine wave.Since the time rate of change of  12 is governed by a balance between the growth and the elastic recoil, the deviation of  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  ( ) within the constitutive model (i.e. ( ) (A − I) versus ( ( )A − I)). Figure 11 shows  22 cos() during one oscillation for the toy sPTT model with  = 1/100 at  = 200 and  = 0.5 for various values of .Clear overshoots of  22 cos() are observed at roughly /  = 0.3 and /  = 0.8 as  → 1.
To highlight more clearly the role of  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 Note that these are just, respectively, the first and second terms on the right hand side of Equations (2.2)b and (2.3)b.Figure 12 shows   ,   , and  12 for the toy sPTT model with  = 1/100 when  = 0.5 and  = 1, respectivly, during a quarter of the oscillation period.During this quarter of the period,  is increasing from 0 to 1. Naturally, stress overshoots are observed in this region for the case when   >   , which is seemingly the case for the toy sPTT model when  → 1 and thus  ( ) = 1.Whilst it is difficult to assess explicitly the role that  22 plays in the evolution of   and   due to the strong coupling in the equations, it is clearly the case that the sharp decrease in   , which can only be caused by a change in  22 , occurs before any observed decrease in   , which allows for a region where   is significantly larger than   .In Figure 12(b), the point of intersection of   ( ) and   ( ) of course defines the exact position for the maxima of  12 (and hence  ,12 ) associated with the stress overshoot.As  → 1, there is a large region where   ≈   , which as mentioned might represent the system approaching steady-shear.
Figure 13 shows   ,   , and  12 for the FENE-P model (or toy FENE-P model with  = 1) with  2 = 100 at  = 1 and  = 100.Again, the decrease in   , caused by the time-dependence of  22 , is significant and is primarily responsible for the generation of the pronounced stress overshoots.  grows much more sharply for the FENE-P model than for the toy sPTT model due to the difference between  ( ) FP and  ( ) sPTT .The overshoots in the FENE-P response are likely exacerbated somewhat by this.We should note that if  ( ) sPTT is replaced by  ( ) FP in the original sPTT model, stress overshoots can still occur even when  22 = 1 solely due to the increased non-linearity of   .However, these overshoots are significantly smaller than those observed when   is non-linear.This replacement of  ( ) sPTT by  ( ) FP in the sPTT model essentially corresponds to a toy FENE-CR model, which we explore in the Appendix A.
As discussed in Section 3.1, the standard sPTT model has universal solutions both in steadyshear and in LAOS (for constant values of ) for constant values of  √ .However, the FENE-P model LAOS response only has universal solutions for constant values of /() for large enough values of  2 .In Section 3.1, we show that replacing  ( ) sPTT with  ( ) 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  ( ) used.Here, we use the toy models to also investigate the effect of the positioning of  ( ) on the scaling of the LAOS response (i.e. ( ) (A − I) versus ( ( )A − I)).In Figure 14, we show the LAOS response of the toy sPTT model with  = 1 (such that  ( ) = 1) for  = 0.5 and  √  = 20 for various values of .The solution seems to only be universal for small enough values of , which is strongly reminiscent of the FENE-P model response presented in Section 3.1.Note that  = 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 (2006) and Davoodi et al. (2022), the generic network model (Equation (1.15)) can be written for steady and homogeneous flows in stress tensor form as Note here that both the time and the velocity gradient are made dimensionless by , and so the entire upper-convected time derivative is multiplied by .For  ( ) ≠  ( ), the last term on the right hand side of Equation (3.12) is nonzero.For the toy sPTT model with  = 1,  (  ) =  (  ) sPTT and  (  ) = 1.In this case, the solution of Equation (3.12) in SSSF is given by the following system 1 For finite values of   (or  √ ), but in the limit that  → 0, the term − ( ,11 + ,22 + ,33 ) on the right hand side of the diagonal components approaches zero, in which case,  ,22 =  ,33 → 0 and Equation 3.13 reduces to which is the solution for the original sPTT model.As discussed in Section 3.1, the solution to Equation 3.14 only depends on the parameter  √  (something which is also seen by introducing the new variable  =    ,11 /(1 − ) into Equation (3.14)).It highlights that the breakdown of the universal scaling with

√
can be caused, at least in SSSF, by setting  ≠  even when  ( ) sPTT is used for the extensibility function.Therefore, considering also the scaling analysis for the FENE-CR model presented in Section 3.1, the breakdown of the universal scaling for low values of  2 in the FENE-P LAOS response is likely a combined effect of the functional form of  ( ) and its position in the conformation tensor form of the constitutive model (i.e. ( ) (A − I) versus ( ( )A − I)).We note briefly that, for the FENE-P model, the extra term on the right hand side of Equation (3.12) is multiplied by the Lagrangian derivative of 1/ (  ) due to the presence of  ( ) in the   -A relationship, which is why this term does not affect the SSSF solution for the FENE-P model.
For clarity, we will briefly summarise this subsection.When  ( ) ≠  ( ) in the network model framework, temporal changes in  22 cause both non-linear growth and non-linear recoil of  12 simultaneously.This causes a region where   is significantly larger than   which manifests as pronounced stress overshoots in the Lissajous curves.The presence of  ( ) in the   -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  ( ) used in the model and the fact that  ( ) is on the outside of the brackets in the recoil term (or  ( ) =  ( ) in the network model framework).The FENE-P response does not scale universally in LAOS due to the specific form of  ( ) used and the fact that  ( ) is on the inside of the brackets in the recoil term (or  ( ) ≠  ( ) in the network model framework).

Sequence of Physical Processes (SPP)
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 SPP analysis (Rogers et al. 2011), which will now be briefly outlined.
Whilst the moduli  In this approach, each point in the oscillation is given by a position vector in the strain-ratestress space () = ⟨(), (),  ,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  () = ()/| ()| points in the direction of the instantaneous trajectory (the overdot here denoting total differentiation with respect to time).The normal vector () =  ()/|  ()| points to the center of the curvature of the instantaneous trajectory.Finally, the bi-normal vector () =  ()×() 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  = 1 in our case due to the equations being solved in non-dimensional form) using the components of  as where   ,12 is the displacement stress.For a more detailed explanation of the SPP framework, the reader is referred to the works of Lee & Rogers (2017); Rogers (2012Rogers ( , 2017)); Rogers et al. (2011).
The time-dependent behaviour of   can be thought of as representing viscous backflow (Choi et al. 2019;Donley et al. 2022;Rogers 2017).In this subsection, we perform the SPP analysis for the sPTT (FENE-P) responses for a fixed  (/) = 0.5 with  2 = 1/ = 100.We also investigate the responses of the toy models outlines in Section 3.4.The SPP 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 3D Lissajous curves for the sPTT and FENE-P models as well as the Cole-Cole plots ( ′′  versus  ′  ).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 (Erturk et al. 2022;Park & Rogers 2020).As

√
is increased the deltoids increase in size, which physically might represent an increase in the degree of microstructural change during the oscillation (Park & Rogers 2020).For the FENE-P response, as  is increased (note  2 is fixed at 100 and we do not assume the solution truly scales with /()) the Cole-Cole plots resemble instead arrow-head shapes which are significantly larger in size than the deltoids observed in the sPTT response (see the scales of the axes).
Figure 16   is sharp enough that it becomes negative in this region before the stress overshoot.Then, similarly to the sPTT response, a decrease in  ′  with an increase in  ′′  is observed.This behaviour is, however, more extreme for the FENE-P response than for the sPTT response.Then, an increase in  ′  with a decrease in  ′′  is observed as the trajectory passes through the sharp stress overshoot.The decrease in  ′′  in this region for the FENE-P response is so large that  ′′  is negative after the stress overshoot.In the supplementary animations Movie 1 and Movie 2, we show the evolution of the Frenet-Serret frame along the Lissajous curves displayed in Figure 16, along with the projections of  in the  and  ,12 plane (highlighting the sign of  ′′  ) and the current position in the Cole-Cole plots.We now use the toy models outlined in Section 3.4 to further highlight the link between the nature of the constitutive model and its behaviour in terms of the SPP analysis.Figure 17   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  has a seemingly minor effect on the Cole-Cole plots.The general shape remains fairly constant, however, a region of negative  ′′  develops which corresponds to the developing stress overshoots (seen in Figure 9).The link between this negative region of  ′′  and the stress overshoot is seen clearly in the supplementary animation (Movie 1) for the FENE-P model response.Essentially, after the stress overshoot, the normal vector  points in the positive stress direction as the recoil fades, which points the bi-normal vector  towards the negative strain rate direction.For the toy FENE-P model, the value of  has a significant effect on the Cole-Cole plot.For  → 1, the extremities of  ′  and  ′′  become significantly large in magnitude.This indicates that the SPP analysis is highly sensitive to the presence of  ( ) in the   -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  ( ) in the   -A relationship appears to give rise to a region of backflow (i.e.negative  ′′  ) before the onset of the stress overshoot.This behaviour cannot be easily explained in terms of the evolution of A, as the region of backflow after the stress overshoot

1D modelling of LAOS
In the previous subsections, it is assumed that the shear rate is uniform in space.For constitutive models which have non-monotonic stress-shear rate relationships (such as Giesekus (Giesekus 1982) or Rolie-Poly (Likhtman & Graham 2003)), 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 or sPTT models 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 non-linear 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 which have monotonic underlying stress-shear rate constitutive curves (Adams & Olmsted 2009; Carter et al.Using the 1D 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 shearbanding.However, if shear-banding occurs, 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 et al. 2000;Olmsted 2008) and discontinuous strain rates in shear banding, we add a small diffusive term (∇ 2 A) to the right hand side of the constitutive models during the simulations.The value of  was fixed at 10 −9 .Note that  here is dimensionless and is given by  = / 2 where  is the diffusion coefficient in  2 /.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  2 1/ = 100 models, respectively, where the stress response has been computed with the 0D approximation (black dashed lines) and 1D simulation (solid lines).For both models, both the 0D and 1D methods give practically identical responses, indicating that the 0D approximation is adequate for properly describing the model responses and that no shear-banding is occurring in the 1D simulations.Since the methodology used for the 1D 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 non-linear in LAOS than it is in SSSF at least in the range of / and /() investigated.

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 non-linear 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 primarily due to  ,12 vs  the extensibility function being multiplied with A instead of (A − 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  22 and  33 during the oscillation which, due to the use of the upper-convected derivative, causes the elastic recoil of  12 to exceed the growth of  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 inhomogenous) flows, and so the differences between the model responses such as the stress overshoots will also occur in Eulerian steady flows which 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 previously been shown analytically for the FENE-P and sPTT models that the stress-strain rate curves scale with /() and  √ , respectively, for simple steady-shear, we have been able to show with our numerical results that the sPTT LAOS response also scales with

√
, but the FENE-P response in LAOS only scales with /() at large enough values of  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   result obtained from the SPP analysis is that the FENE models (both FENE-P and FENE-CR) exhibit backflow (i.e.negative  ′′  ) 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 1D 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 1D and 0D solutions were compared.with the non-linearity of   and thus the transient nature of  22 .This becomes particularly clear when Figures 23(c) and 10(f) are compared.Note when comparing these two figures that  ( ) is moving from the inside to the outside of the brackets in the recoil term, and  ( ) does not appear in the   -A relationship in either case.
Figure 24 shows the Cole-Cole plots for the toy FENE-CR model responses displayed in Figure 23.For  > 0, a region of backflow is observed before the small stress overshoot, which is also observed in the toy FENE-P model response for  > 0. This further highlights that the SPP framework can clearly identify the presence of  ( ) in the   − A relationship, and thus help to identify a suitable form of constitutive model from LAOS data.

Figure 1 :
Figure 1: Viscous Lissajous-Bowditch plots in / () − / () space for the FENE-P (sPTT) model.Black curves represent the UCM response.Black numbers in each plot represent the maximum value of  ,12 in the UCM response, since the -axis is scaled differently in each plot.

Figure 2 :
Figure 2: Elastic Lissajous-Bowditch plots in / () − / () space for the FENE-P (sPTT) model.Black curves represent the UCM response.Black numbers in each plot represent the maximum value of  ,12 in the UCM response, since the -axis is scaled differently in each plot.

Figure 3 :
Figure 3: Viscous Lissajous-Bowditch plots in /() ( √ ) − / () space for the FENE-P (sPTT) model.Blue numbers in each plot represent the maximum value of  ,12 in the sPTT response since the -axis is scaled differently in each plot.Plots in the black dashed box are shown at a larger scale below in Figure 5.

Figure 5 :
Figure 4: Elastic Lissajous-Bowditch plots in /() ( √ ) − / () space for the FENE-P (sPTT) model.Black curves represent the UCM response.Blue numbers in each plot represent the maximum value of  ,12 in the sPTT response, since the -axis is scaled differently in each plot. d

Figure 10 :
Figure 10: Viscous Lissajous curves for the toy FENE-P model ( 2 = 100) with various values of  between 1 and 0.  = 1 and  = 100.-axis limits are shown by the numbers adjacent to the ends of the axes. axes run from -1.04 to 1.04.

Figure 13 :
Figure 13:   ,   , and  12 vs  in one quarter of an oscillation for the FENE-P model ( 2 = 100) for  = 1 and  = 100.Dashed vertical line shows the cross-over point of   and   , and hence also the maxima of  12 .
independent during SAOS, and their values can be found from simple Fourier analysis, these moduli are transient during LAOS for a non-linear viscoelastic material or constitutive model.The SPP framework identifies, at each instant, these transient moduli, denoted as  ′  () and  ′′  (), by utilising the Frenet-Serret theorem.

′
and  ′′  during the (half) oscillation inform us of underlying physical phenomena occurring in the stress response.Increasing values of  ′  represents intra-cycle strain-hardening whilst decreasing values of  ′  represents intra-cycle strain-softening.Similarly, for the viscous modulus, increasing values of  ′′  represents intracycle shear-thickening whereas decreasing values of  ′′  represents intra-cycle shear-thinning.Negative instantaneous values of  ′  can be thought of as representing elastic recoil, and negative instantaneous values of  ′′
shows the Cole-Cole plots obtained from the SPP framework for the (a) toy FENE-P and (b) toy sPTT models.Note that  = 0.5 and  = 200 for the toy sPTT model, whilst  = 1 and  = 100

Figure 18 :
Figure 18: Cole-Cole plots comparing the toy sPTT model response with  = 1 and the toy FENE-P model response with  = 0.

Figure 19 :
Figure 19: Lissajous-Bowditch plots (viscous projection) for the sPTT model.Blue solid lines show the 0-D approximation solution from the previous section, black dashed lines show the results from the 1D simulations where the stress is computed at the top (moving) boundary.

′
and  ′′  .The Cole-Cole plots show that the FENE-P model exhibits significantly more complex rheological behaviour during the oscillation.A key  ,12 vs

Figure 20 :
Figure 20: Lissajous-Bowditch plots (viscous projection) for the FENE-P model ( 2 = 100).Red solid lines show the 0D approximation solution from the previous section, black dashed lines show the results from the 1D MOL simulations where the stress is computed at the top (moving) boundary.