# Laun's rule for predicting the first normal stress coefficient in complex fluids: A comprehensive investigation using fractional calculus

Mohua Das, Joshua David John Rathinaraj, Liviu Iulian Palade, Gareth H. Mckinley

► To cite this version:

Mohua Das, Joshua David John Rathinaraj, Liviu Iulian Palade, Gareth H. Mckinley. Laun's rule for predicting the first normal stress coefficient in complex fluids: A comprehensive investigation using fractional calculus. Physics of Fluids, 2024, 36 (1), 10.1063/5.0179709. hal-04662319

HAL Id: hal-04662319

https://hal.science/hal-04662319v1

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Laun's rule for predicting the first normal stress coefficient in complex fluids: A comprehensive investigation using fractional calculus 🧠

Special Collection: Tanner: 90 Years of Rheology

Mohua Das 🐦 ; Joshua David John Rathinaraj 🐦 ; Liviu Iulian Palade 🐦 ; Gareth H. McKinley FRS 🐦

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

Check for updates

Physics of Fluids 36, 013111 (2024)

https://doi.org/10.1063/5.0179709

<figure><img src="image_5.png" /></figure>

View
Online

<figure><img src="image_6.png" /></figure>

Export
Citation

## Articles You May Be Interested In

Why the Cox–Merz rule and Gleissle mirror relation work: A quantitative analysis using the Wagner integral framework with a fractional Maxwell kernel

*Physics of Fluids* (March 2022)

A constitutive model for entangled polydisperse linear flexible polymers with entanglement dynamics and a configuration dependent friction coefficient. Part II. Modeling “shear modification” following cessation of fast shear flows

J. Rheol. (January 2018)

Evaluating the exit pressure method for measurements of normal stress difference at high shear rates

J. Rheol. (May 2020)

<figure><img src="image_7.png" /></figure>

AIP
Publishing

Physics of Fluids

Special Topics Open
for Submissions

Learn More

<figure><img src="image_8.png" /></figure>

AIP
Publishing
# Laun's rule for predicting the first normal stress coefficient in complex fluids: A comprehensive investigation using fractional calculus

F

Cite as: Phys. Fluids **36**, 013111 (2023); doi: 10.1063/5.0179709

Submitted: 4 October 2023 · Accepted: 14 November 2023 ·

Published Online: 11 January 2024

<figure><img src="image_9.png" /></figure>

View Online

<figure><img src="image_10.png" /></figure>

Export Citation

<figure><img src="image_11.png" /></figure>

CrossMark

Mohua Das,¹,ª ID Joshua David John Rathinaraj,¹,ᵇ ID Liviu Iulian Palade,²,ᶜ ID and Gareth H. McKinley FRS¹,ᵈ ID

## AFFILIATIONS

¹Hatsopoulos Microfluids Laboratory, Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA

²Université de Lyon, CNRS, Institut Camille Jordan UMR 5208, INSA-Lyon, 21 Avenue Jean Capelle, 69621 Villeurbanne, France and Université de Lyon, CNRS, LaMCoS UMR 5259 and Pôle de Mathématiques, INSA-Lyon, 69621 Villeurbanne, France

Note: This paper is part of the special topic, Tanner: 90 Years of Rheology.
Dedicated to Roger I. Tanner FRS on his 90th Birthday.

a) Electronic mail: mohuadas@mit.edu

b) Electronic mail: joshuajr@mit.edu

c) Electronic mail: liviu-iulian palade@insa-lyon.fr

d) Author to whom correspondence should be addressed: gareth@mit.edu

## ABSTRACT

Laun's rule [H. M. Laun, "Prediction of elastic strains of polymer melts in shear and elongation," J. Rheol. **30**, 459–501 (1986).] is commonly used for evaluating the rate-dependent first normal stress coefficient from the frequency dependence of the complex modulus. We investigate the mathematical conditions underlying the validity of Laun's relationship by employing the time-strain-separable Wagner constitutive formulation to develop an integral expression for the first normal stress coefficient of a complex fluid in steady shear flow. We utilize the fractional Maxwell liquid model to describe the linear relaxation dynamics compactly and accurately and incorporate material nonlinearities using a generalized damping function of Soskey–Winter form. We evaluate this integral representation of the first normal stress coefficient numerically and compare the predictions with Laun's empirical expression. For materials with a broad relaxation spectrum and sufficiently strong strain softening, Laun's relationship enables measurements of linear viscoelastic data to predict the general functional form of the first normal stress coefficient but often with a noticeable quantitative offset. Its predictive power can be enhanced by augmenting the original expression with an adjustable power-law index that is based on the linear viscoelastic characteristics of the specific material being considered. We develop an analytical expression enabling the calculation of the optimal power-law index from the frequency dependence of the viscoelastic spectrum and the strain-softening characteristics of the material. To illustrate this new framework, we analyze published data for an entangled polymer melt and for a semiflexible polymer solution; in both cases our new approach shows significantly improved prediction of the experimentally measured first normal stress coefficient.

© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0179709

## I. INTRODUCTION

One key rheological phenomenon observed in a wide range of complex materials, including polymers, printing inks, cosmetics, food products, and other soft solids, is the presence of normal stress differences during simple shearing deformations.² Normal stress differences act in the transverse directions to the flow direction and can be considered as the fluidic counterpart of the “Poynting effect” observed in nonlinear elasticity.³,⁴ As the landmark text on *Engineering Rheology*

by Tanner⁴,⁵ explains, these nonlinear normal stress differences cause various rheological phenomena, such as extrudate swelling, the Weissenberg rod climbing effect, onset of purely elastic instabilities, etc.²⁻⁵,⁷ Thus, determining the magnitude of the normal stress differences in a complex fluid system holds significance not just in the design of various technological processes like extrusion and injection molding,⁴ but also in understanding a range of natural phenomena.⁸,⁹
To quantify the magnitude of the normal stress differences experienced by a complex material under simple shear flow with $\dot{\gamma}_{yx} = \partial v_x / \partial y$, two functions known as the first and second normal stress differences ($N_1$ and $N_2$) are typically measured. These quantities are defined as follows:⁴,¹⁰

$$N_1 = \sigma_{xx} - \sigma_{yy}, \quad \text{and} \quad N_2 = \sigma_{yy} - \sigma_{zz}. \qquad (1)$$

The magnitudes of $N_1$ and $N_2$ (compared to the corresponding shear stress, $\sigma_{yx}$) helps to distinguish between strongly viscoelastic and predominantly viscous (inelastic) behaviors.¹¹ Inelastic liquids exhibit immeasurably small normal stress differences²,¹⁰ and can be modeled much more simply as generalized Newtonian fluids.⁴

Specifically, in the context of viscoelastic fluids, the first normal stress difference, $N_1$, provides valuable information about the material's nonlinear elasticity during shear and can be viewed as generating a net tension along streamlines in a flowing material. In the 1970s, Roger Tanner and colleagues elegantly showed using ideas from nonlinear rubber elasticity how this idea of streamline tension could be used to make an *a priori* estimate of the extent of die swell to be expected in a viscoelastic material exiting an extruder nozzle.¹²

In the case of polymeric fluids, the experimentally determined magnitude of $N_2$ is noticeably smaller than that of $N_1$ (typically ranging between 10% and 30% of $N_1$ for steady, simple shear flow).⁴,¹⁰ Consequently, as Tanner has noted, accurately measuring $N_2$ presents experimental challenges as it is much smaller than $N_1$.¹³ One elegant way of effectively isolating and demonstrating the importance of second normal stresses was pioneered over 50 years ago by Tanner and coworkers and is known as the “tilted trough” experiment.¹⁴ In recent years, this apparatus has attracted renewed interest for helping in unambiguously distinguishing the signs and magnitudes of the second normal stress difference in dense suspensions under shear.¹⁵⁻¹⁷ For highly filled systems with viscoelastic matrices, it is often found that $|N_2| > |N_1|$.¹⁷,¹⁸ For a comprehensive discussion of the current state of understanding, see the review of Maklad and Poole.¹⁹ Thankfully, in a wide-range of unfilled polymeric systems, the effects of $N_2$ are less significant in many flows of practical importance,¹⁰,¹⁷ and in the rest of this work, we focus on the magnitude and functional form of the first normal stress difference.

Notably, for simple fluids with fading memory, at very low shear rates $N_1$ increases proportionally to the square of the shear rate ($N_1 \propto \dot{\gamma}^2$).⁴,¹⁰ Consequently, by scaling $N_1$ and $N_2$ with the imposed shear rate, two material functions known as the first (or primary) and second normal stress coefficients ($\Psi_1$ and $\Psi_2$) can be defined as follows:²,¹⁰,¹⁴,²⁰

$$\Psi_1 = \frac{N_1(\dot{\gamma})}{\dot{\gamma}^2}, \quad \text{and} \quad \Psi_2 = \frac{N_2(\dot{\gamma})}{\dot{\gamma}^2}. \qquad (2)$$

Measuring the nonlinear rheological behavior of highly viscoelastic materials using conventional rheometers (especially at the high shear rates encountered in some industrial scenarios²¹) presents several challenges, including edge fracture of the free surface at high shear rates²² and limitations in the sensitivities of torque and axial force transducers.¹⁴,²³ As a result of these difficulties, researchers have attempted to establish relationships between nonlinear viscoelastic properties, such as the rate-dependent shear viscosity and the first normal stress difference, and linear viscoelastic data, considering the

relative convenience of obtaining accurate linear viscoelastic data across a wide range of frequencies and time scales.

One of the earliest and most widely recognized relationships between linear viscoelastic data and the steady shear viscosity function was noted by Cox and Merz in 1958.²⁴ The Cox–Merz rule establishes an approximate equality between the nonlinear shear-rate-dependent steady shear viscosity, $\eta(\dot{\gamma})$, and the magnitude of the complex viscosity, $\eta^*(\omega)$, as follows:

$$\eta(\dot{\gamma})|_{\dot{\gamma}=\omega} = |\eta^*(\omega)| = \frac{G''(\omega)}{\omega} \left[ 1 + \left( \frac{G'(\omega)}{G''(\omega)} \right)^2 \right]^{0.5}, \qquad (3)$$

where $G'(\omega)$ and $G''(\omega)$ represent the frequency-dependent storage and loss moduli, respectively.

Later, in 1986, Laun proposed a similar empirical relation that connects the primary normal stress coefficient in a steady shear flow to the components of complex modulus.¹ This relationship is also discussed in Tanner's book⁴ and can be expressed as follows:

$$\Psi_1(\dot{\gamma})|_{\dot{\gamma}=\omega} = \frac{2G'(\omega)}{\omega^2} \left[ 1 + \left( \frac{G'(\omega)}{G''(\omega)} \right)^2 \right]^{0.7}. \qquad (4)$$

Alternatively, Eq. (4) can be written in the following form:

$$\Psi_1(\dot{\gamma})|_{\dot{\gamma}=\omega} = \frac{2G'(\omega)}{\omega^2} (\mathrm{cosec}(\delta))^{1.4}, \qquad (5)$$

where the phase angle, $\delta = \tan^{-1}(G''(\omega)/G'(\omega))$. Notably, at vanishingly low frequency and shear rate, Eq. (4) satisfies the second-order-fluid limit expected for simple fluids with fading memory,⁵ i.e.,

$$\Psi_{1,0} = \lim_{\omega \to 0} \left[ \frac{2G'(\omega)}{\omega^2} \right]. \qquad (6)$$

Laun's empirical rule has proven useful in predicting the first normal stress function of concentrated solutions and melts of various commercial polymers like polyethylene, polypropylene, polystyrene, and polybutadiene⁶ from measurements of the linear complex modulus. However, there has been ongoing discussion regarding the constant value of the power-law index selected by Laun in Eq. (4), i.e., the value 0.7. Several researchers have found different values for the power-law index when examining the consistency of experimental data, as discussed below.

To assess the suitability of Laun's empirical relationship for entangled polymer liquids, Wen et al.²⁵ conducted a study involving both monodisperse and a series of bidisperse solutions of polystyrene. These solutions were created by dissolving polystyrene with varying molecular weights in either diethyl phthalate or tricresyl phosphate. A significant amount of cosolvent, methylene chloride, was added to ensure complete dissolution and homogeneity. The results of this study reveal that the power-law index (0.7) proposed by Laun is a reliable approximation for polydisperse systems. However, for the monodisperse system, an exponent of 0.5 was found to be a more accurate choice than 0.7. A number of several other researchers, such as Baird,²¹ Teixeira et al.,²⁶ Kádár et al.,²⁷ and Tang et al.,⁹ have reported diverse power-law index values, ranging between 0 and 0.7, to achieve a satisfactory agreement between experimentally measured data and Laun's empirical expression.

The diverse findings of these previous systematic experimental studies inspired us to evaluate the accuracy and range of validity of
Laun's empirical relationship. To investigate the accuracy of Laun's prediction, we considered two distinct materials as representative and diverse samples: an entangled polymer melt (melt of low-density polyethylene, referred to as LDPE Melt I) and an entangled semiflexible polymer solution (a solution of 0.5 wt. % Xanthan gum in water). We referred to two primary sources for collecting the needed experimental data for the primary normal stress coefficients, storage and loss moduli of the selected samples. The data for LDPE Melt I was digitized from Laun's 1986 paper,¹ while the data for the 0.5 wt. % Xanthan gum solution was obtained from the study by Jaishankar and McKinley.²⁸ Laun's predicted values of the normal stress coefficient were calculated by utilizing the collected data for storage and loss moduli, employing Eq. (4).

The comparisons between the measured and predicted values of the primary normal stress coefficient for both samples are presented in Fig. 1. The results shown in Fig. 1(a) demonstrate that for an entangled and polydisperse melt (LDPE Melt I), the predictions using Laun's empirical formula exhibit good agreement with the experimental data. However, for the 0.5 wt. % Xanthan gum solution, a consistent offset is observed between the measured data and the predicted values, which becomes more prominent with increasing shear rate [Fig. 1(b)]. Consequently, it can be concluded that the power-law index chosen by Laun in his empirical formula [Eq. (4)] does not accurately predict the normal stress coefficient data for the semiflexible concentrated Xanthan gum solution.

To generalize Laun's relationship, in this study, we introduce a simple modification. Specifically, we replace the constant value of the power-law index (0.7) in Eq. (4) with a variable exponent denoted by a variable real number $l$. Equation (4) can then be expressed in the following generalized form:

$$ \Psi_1(\dot{\gamma})|_{\dot{\gamma}=\omega} = \frac{2G'(\omega)}{\omega^2} \left[ 1 + \left( \frac{G'(\omega)}{G''(\omega)} \right)^2 \right]^l \quad (7) $$

<figure><img src="image_12.png" /><figcaption>FIG. 1. Comparison of primary normal stress coefficient data (represented by filled purple circles) (a) for an entangled polyethylene melt as measured by reproduced with permission from H. M. Laun, J. Rheol. 30, 459–501 (1986). Copyright 1986 Society of Rheology¹ and (b) by Jaishankar and McKinley for a concentrated solution of a semiflexible polymer (reproduced with permission from Jaishankar and McKinley, J. Rheol. 58, 1751–1788 (2014). Copyright Society of Rheology.²⁸). The values of primary normal stress coefficients predicted by Laun's empirical relationship using dynamic moduli data from the aforementioned papers are represented by filled blue squares joined by a dotted line.</figcaption></figure>

In a recent study, Rathinaraj et al.²⁹ used continuum mechanics and fractional calculus to derive the underlying conditions for which empirical rules, such as the Cox–Merz rule and the Gleissle mirror relationship, are satisfied. Taking inspiration from their work, the present study aims to provide a detailed analysis of the conditions determining the validity of Laun's empirical relationship for complex fluids. Following the framework established by Rathinaraj et al.,²⁹ we employ a similar approach to assess the validity of Laun's rule and the range of the values of the exponent $l$.

To begin, we transform the storage and loss moduli terms in Eq. (7) into analytical expressions based on the parameters of the fractional Maxwell liquid (FML) constitutive model. Subsequently, we derive an integral expression for the primary normal stress coefficient using Wagner's integral formulation, which is a time-strain separable (TSS) form of the Kaye–Bernstein Kearsley Zapas (K-BKZ) constitutive equation.²⁹,³⁰–³² Within this framework, the relaxation kernel can be expressed in terms of the Mittag–Leffler function to predict the linear relaxation dynamics in arbitrary deformations. Additionally, we incorporate nonlinear strain-softening at large strain amplitudes using the Soskey–Winter form of the damping function,³³ which is a modified version of the sigmoidal form originally proposed by Papanastasiou et al.³⁴

By comparing the functional forms and numerical values of the primary normal stress coefficients obtained using Eq. (7) plus the FML model parameters and the Wagner integral formulation, we explore the material properties essential for more accurately determining the value of the power-law index ($l$) in Laun's relationship.

We perform asymptotic analyses at very high and very low shear rates to gain further insights. These analyses show that complex fluids with strong strain-softening at large strain amplitudes and a broad relaxation spectrum do generally follow a relationship of the form suggested by Laun. We specify limits based on fluid rheological properties to distinguish between fluids with sufficiently strong and with insufficient (or weak) damping. For complex fluids with sufficiently strong

<figure><img src="image_13.png" /></figure>
damping and a broad relaxation spectrum, we present a more comprehensive framework for selecting the optimal power-law index, $l_{opt}$, taking into account both the linear viscoelastic response and the nonlinear damping characteristics of a material. To validate the effectiveness of this framework, we replot the normal stress coefficient data for LDPE Melt I and the 0.5 wt. % Xanthan gum solution using the optimal value, $l_{opt}$, obtained from our proposed framework. By quantitatively comparing these results with the measured primary normal stress coefficient data, we demonstrate the efficacy of the proposed framework in accurately predicting the primary normal stress coefficient of complex fluids.

## II. MATHEMATICAL FORMULATION

The fractional Maxwell liquid (FML) model has been found to be highly effective in compactly and accurately describing the linear rheological behavior of numerous viscoelastic liquids.²⁸,³⁵ It is also thermodynamically consistent³⁶ and predicts a bounded viscosity $\eta_0$ in the limit of zero shear rate.³⁷ The corresponding objective constitutive equation³⁸ has also been proven to have suitable stability properties (Hadamard well-posedness, stable solutions).³⁹,⁴⁰ The three-parameter FML model is an important limiting case of the more general four-parameter fractional Maxwell model (FMM). It consists of a single fractional “spring-pot”²⁸ and a Newtonian dashpot connected in series.³⁷ The one-dimensional constitutive relationship for the FML model can be represented in the following form:²⁹

$$\sigma(t) + \frac{\eta_0}{\mathbb{G}} \frac{d^{1-\beta}}{dt^{1-\beta}} \sigma(t) = \eta_0 \frac{d\gamma(t)}{dt}. \quad (8)$$

In Eq. (8), $\eta_0$ represents the zero shear viscosity of the fluid, and the material parameter, $\mathbb{G}$, is best characterized as a quasi-property with dimensions $(\text{Pa} \cdot \text{s}^\beta)$.²⁸,²⁹ The value of the exponent $\beta$ is constrained in the range of $0 \le \beta < 1$.²⁸

Considering the physical dimensions of Eq. (8), we can identify the single characteristic relaxation time (with SI units of seconds) of the FML model to be²⁹

$$\tau_c = (\eta_0 / \mathbb{G})^{\frac{1}{1-\beta}}. \quad (9)$$

This represents the appropriate average value of a broad underlying continuous relaxation spectrum.²⁷

By taking the Fourier transform of Eq. (8), we obtain an expression for the complex modulus, $G^*(\omega)$. Separating the real and imaginary parts, we obtain expressions for the equations of the storage modulus and loss modulus, respectively,

$$\frac{G'(\omega)}{G_c} = \frac{(\text{De})^{2-\beta} \cos(\pi\beta/2)}{1 + 2(\text{De})^{1-\beta} \cos(\pi(1-\beta)/2) + (\text{De})^{2(1-\beta)}}, \quad (10)$$

$$\frac{G''(\omega)}{G_c} = \frac{\text{De} + (\text{De})^{2-\beta} \sin(\pi\beta/2)}{1 + 2(\text{De})^{1-\beta} \cos(\pi(1-\beta)/2) + (\text{De})^{2(1-\beta)}}, \quad (11)$$

where $G_c = \eta_0 \tau_c^{-1} = \mathbb{G} \tau_c^{-\beta}$ represents the characteristic elastic modulus of the FML model,²⁹ and $\text{De} = \tau_c \omega$ denotes the Deborah number.²⁹

Using the expressions from Eqs. (10) and (11), the generalized form of Laun's rule for the expected functional dependence of the primary normal stress coefficient [Eq. (7)] can be expressed in a non-

dimensional form in terms of the FML constitutive parameters $(\beta, \mathbb{G}, \eta_0)$ as

$$\frac{\Psi_1(\dot{\gamma})}{\eta_0 \tau_c} \equiv \hat{\Psi}_1(\text{De}) = \frac{2(\text{De})^{-\beta} \cos(\pi\beta/2)}{1 + 2(\text{De})^{1-\beta} \cos(\pi(1-\beta)/2) + (\text{De})^{2(1-\beta)}} \times \left[ 1 + \left\{ \frac{\cos(\pi\beta/2)}{(\text{De})^{\beta-1} + \sin(\pi\beta/2)} \right\}^2 \right]^t. \quad (12)$$

As noted earlier, we can derive an integral expression for the primary normal stress coefficient of a nonlinear viscoelastic material by utilizing the time-strain separable Wagner integral model.⁴¹ In this formulation, the tensorial stress-strain relationship is defined as follows:⁹,³¹,³²

$$\sigma(t) = \int_{-\infty}^{t} m(t-t')h(\gamma)(\mathbb{C}^{-1} - \mathbb{I})dt', \quad (13)$$

where $m(t-t')$ is a linear viscoelastic property, known as the hereditary or memory function of the material.³⁰,⁴² This function characterizes the time-dependent fading memory of the material, considering its deformation history. It can be obtained by taking the first derivative of the relaxation modulus, $G(t-t')$, with respect to the past time $t'$ as follows:⁹

$$m(t-t') = \frac{\partial G(t-t')}{\partial t'}. \quad (14)$$

Conversely, the linear viscoelastic relaxation modulus can be expressed as

$$G(t-t') = \int_{-\infty}^{t} m(t-t')dt'. \quad (15)$$

For a step-strain response, the analytical expression for the relaxation modulus, $G(t-t')$, for the FML model is given by Rathinaraj et al.,²⁹

$$G(t-t') = \mathbb{G}(t-t')^{-\beta} E_{1-\beta, 1-\beta} \left( -\frac{\mathbb{G}}{\eta_0} (t-t')^{1-\beta} \right). \quad (16)$$

where $E_{a, b}(z)$ is the two-parameter Mittag-Leffler function (MLF).²⁸,²⁹ Therefore, by combining Eqs. (14), (16), and also Eq. (9), we can derive the following expression for the memory function of the FML model:²⁹,⁴³

$$m(t-t') = -\mathbb{G}(t-t')^{-(1+\beta)} E_{1-\beta, -\beta} \left( -\left(\frac{t-t'}{\tau_c}\right)^{1-\beta} \right). \quad (17)$$

The function $h(\gamma)$ in Eq. (13) represents the damping function, which enables the description of the nonlinear response of viscoelastic materials under large strains.⁴⁴ It is defined as the ratio of the strain-dependent nonlinear relaxation modulus $G(t, \gamma)$ to the linear relaxation modulus $G(t)$

$$h(\gamma) = \frac{G(t, \gamma)}{G(t)}, \quad (18)$$

in fluids with fading memory. Values of $h(\gamma)$ begin at unity and ultimately decrease toward zero.⁴⁴ When $h(\gamma) = 1$, Eq. (13) reduces to Lodge's rubber-like liquid (RLL) constitutive equation.⁴⁴ Therefore, $h(\gamma) = 1$ denotes a quasi-linear response at large strains, and $h(\gamma) < 1$ represents softening of the material at large deformation, leading to a reduction in the value of the strain-dependent relaxation modulus. In this study, we consider a generalized version of the damping function proposed by Soskey and Winter³³
$$h(\gamma) = \frac{1}{1 + a\gamma^b}, \qquad (19)$$

where a and b are the fitting parameters. Rathinaraj and coauthors²⁹ provided a more physically motivated reformulation of Eq. (19) as follows:

$$h(\gamma) = \frac{1}{1 + (\gamma/\gamma^*)^m}. \qquad (20)$$

Equation (20) is a two-parameter equation with a positive damping exponent ($m > 0$) and critical strain denoted by $\gamma^*$ (see Ref. 29). When the imposed strain reaches the critical value, $\gamma \to \gamma^*$, the damping function reduces to 0.5, i.e., $h(\gamma) \to 0.5$ and the relaxation modulus is decreased by 50% from the corresponding linear viscoelastic value.

The term $(C^{-1} - I)$ in Eq. (13) defines the appropriate strain tensor for general large deformations. Here, $C^{-1}$ represents the Finger deformation tensor,⁹ and I denotes the identity tensor so that $\sigma(t)$ represents the additional or “extra” stress tensor resulting from deformation. For a simple shearing flow, the strain accumulated between time $t'$ and $t$ can be written as²,¹⁰

$$C^{-1}(t', t) = \begin{bmatrix} 1 + [\gamma(t) - \gamma(t')]^2 & [\gamma(t) - \gamma(t')] & 0 \\ [\gamma(t) - \gamma(t')] & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}. \qquad (21)$$

In Eq. (21), $[\gamma(t) - \gamma(t')]$ corresponds to the accumulated strain between the past time $t'$ and present time $t$. For steady simple shear, $\gamma(t) - \gamma(t') = \dot{\gamma}_0(t - t')$, where $\dot{\gamma}_0$ is the constant shear rate.²⁹

By substituting Eqs. (17), (20), and (21) into Eq. (13), we obtain an expression for the first normal stress difference of the fractional Wagner model in a steady shear flow:

$$\begin{aligned} N_1(\dot{\gamma}_0) &= \sigma_{11} - \sigma_{22} \\ &= -\mathcal{G} \int_{-\infty}^{t} (t-t')^{-(1+\beta)} E_{1-\beta, -\beta} \left( -\left(\frac{t-t'}{\tau_c}\right)^{1-\beta} \right) \\ & \times \frac{1}{1 + (\dot{\gamma}_0(t-t')/\gamma^*)^m} [\dot{\gamma}_0(t-t')]^2 dt', \end{aligned} \qquad (22)$$

From Eq. (25), we can see that the ratio $f$ depends on four constitutive model parameters: the FML exponent, $\beta$; the critical strain, $\gamma^*$; the damping exponent, $m$; and Laun’s power-law index, $l$.

The slope $\beta$ can be calculated using the frequency dependence of the dynamic moduli data at high frequencies as follows:

$$\beta = \lim_{\omega\tau_c \gg 1} \frac{d(\ln G'(\omega))}{d(\ln \omega)} = \lim_{\omega\tau_c \gg 1} \frac{d(\ln G''(\omega))}{d(\ln \omega)}. \qquad (26)$$

where the negative sign in front of the integral comes from differentiating the Mittag–Leffler function (MLF) that describes the fading memory in the material. Substituting $\mathcal{G} = \eta_0/\tau_c^{1-\beta}$, as well as defining a dimensionless time as $\gamma = (t - t')/\tau_c$, and a Weissenberg number, $Wi = \dot{\gamma}_0\tau_c$, Eq. (22) can be written more compactly in the following form:

$$N_1(\dot{\gamma}_0) = -\eta_0 \tau_c \int_0^\infty \frac{y^{1-\beta} E_{1-\beta, -\beta} (-y^{1-\beta})}{1 + (Wi y/\gamma^*)^m} \dot{\gamma}_0^2 dy. \qquad (23)$$

From the form of the terms in the denominator [Eq. (23)], it is clear that the importance of nonlinearities in the first normal stress difference depends on the relative magnitudes of the Weissenberg number, $Wi$, the critical strain, $\gamma^*$, and the damping exponent, $m$.

Scaling the normal stress with the product $\eta_0\tau_c$, we can also derive an expression for the non-dimensional first normal stress coefficient as a function of the Weissenberg number, $Wi$. The resulting expression is as follows:

$$\frac{\Psi_1(\dot{\gamma}_0)}{\eta_0 \tau_c} \equiv \dot{\Psi}_1(Wi) = - \int_0^\infty \frac{y^{1-\beta} E_{1-\beta, -\beta} (-y^{1-\beta})}{1 + (Wi y/\gamma^*)^m} dy. \qquad (24)$$

Thus, we have obtained two expressions for calculating the non-dimensional first normal stress coefficient. The first one exploits Laun’s rule in terms of the FML model representation of the complex modulus, $\dot{\Psi}_1(\text{De})$ [Eq. (12)]. The second formulation is obtained by using the Wagner integral equation, $\dot{\Psi}_1(\text{Wi})$ [Eq. (24)]. For this second formulation, we employed the time derivative of the two-parameter Mittag–Leffler form of the relaxation modulus [Eq. (17)] to represent linear relaxation dynamics and a generalized damping function of the Soskey–Winter form to allow us to describe the nonlinear deformation. If $l=0.7$ holds for all materials exhibiting time-strain separability, then the ratio of $\dot{\Psi}_1(\text{Wi})/\dot{\Psi}_1(\text{De})$ will be exactly equal to unity. For clarity, in our subsequent discussions, we denote this ratio as $f$ and it is given by the following expression:

$$f \equiv \frac{\dot{\Psi}_1(\text{Wi})}{\dot{\Psi}_1(\text{De})} = \frac{-\int_0^\infty y^{1-\beta} E_{1-\beta, -\beta} (-y^{1-\beta}) \frac{1}{1 + (Wi y/\gamma^*)^m} dy}{\frac{2(\text{De})^{-\beta} \cos(\pi\beta/2)}{1 + 2(\text{De})^{1-\beta} \cos(\pi(1-\beta)/2) + (\text{De})^{2(1-\beta)}} \times \left[ 1 + \left\{ \frac{\cos(\pi\beta/2)}{(\text{De})^{\beta-1} + \sin(\pi\beta/2)} \right\}^2 \right]^l} \qquad (25)$$

For example, in the Rouse model $\beta = 1/2$⁴⁵,⁴⁶ and for the Zimm model $\beta = 2/3$.⁴⁶ In our initial investigations, we thoroughly examined various scenarios by considering specific values for the model parameters. From our comprehensive analysis, it was established that Laun’s rule cannot be satisfied at all deformation rates for materials with a very narrow relaxation spectrum characterized by a single characteristic relaxation timescale, such as the Maxwell spectrum (where $\beta = 0$).
Specifically, in this particular case ($\beta \to 0$), the functions $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$ do converge at low values of Wi and De to constant values that are equal to each other. Consequently, the ratio between them becomes unity, and we define the limiting ratio here as

$$f_0 = \frac{\lim_{\text{Wi}\to 0} \hat{\Psi}_1(\text{Wi})}{\lim_{\text{De}\to 0} \hat{\Psi}_1(\text{De})} = 1. \qquad (27)$$

However, at higher deformation rates, when $\beta \to 0$, the frequency dependence of the denominator and the strain rate dependence of the numerator in Eq. (25) do not agree, and $f \neq 1$. While Laun's relationship can be employed to predict the first normal stress coefficient for low values of Wi and De, this is, in fact, just a restatement of second-order fluid theory, i.e.,

$$\Psi_1(\dot{\gamma}) = \frac{2G'(\omega)}{\omega^2} = 2\eta_0 \tau_c, \qquad (28)$$

which is not valid for higher values of De and Wi.

Furthermore, highly elastic or cross-linked materials without strain-softening such that $h(\gamma) = 1$ at all strains do not conform to Laun's relationship. Detailed discussion of these cases is provided in the Appendix.

To further investigate the behavior of $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$ for different values of the model parameters beyond those specified earlier, we generated plots based on their respective equations [Eqs. (12) and (24)]. To be specific in our discussions, here we report results for $\beta$ values corresponding to the Rouse model ($\beta = 1/2$)$^{45,46}$ and the Zimm model ($\beta = 2/3$).$^{45}$ We first employ constant values of the damping exponent, $m = 2.0$, and critical strain, $\gamma^* = \sqrt{15}/2$, which correspond closely to the damping function predicted by the Doi-Edwards theory with the independent alignment (IA) approximation.$^{44}$ The corresponding plots are shown in Figs. 2(a) and 2(b). We examine two different values of Laun's exponent, i.e., $l=0$ and $l=0.7$. These values represent the minimum and maximum values of the function $\Psi_1(\text{De})$

according to Laun's relationship [Eq. (4)].$^{30}$ In the insets, we show the evolution of the ratio $f$ [Eq. (25)].

As shown in Fig. 2(a), at very low values of De and Wi, all three curves overlap, and the ratio $f$ tends to converge toward the value of 0.92. However, the values of $f$ diverge at higher values of De and Wi. Notably, for the specified model parameters ($\beta, \gamma^*, m$) = (1/2, $\sqrt{15}/2$, 2.0), the agreement between $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$ is better for the $l=0$ value than for $l=0.7$. A similar trend can be observed in Fig. 2(b) as well, where the curves converge at very low De and Wi and diverge at very high De and Wi.

From these representative figures, we can draw some conclusions. First, at extremely low values of De and Wi, the ratio, $f$, is independent of $l$, indicating that it depends only on the values of $\beta, \gamma^*, \text{ and } m$. Second, at higher values of De and Wi, the curves exhibit similar power-law features but are slightly shifted by an offset factor. This suggests that adjusting the exponent $l$ to values other than 0 and 0.7 might make the ratio nearly equal to unity across a wide range of deformation rates, i.e., it may minimize discrepancies between $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$.

The behavior of polymer melts at low shear rates suggests that the first normal stress coefficient ($N_1$) becomes proportional to the square of the shear rate $\dot{\gamma}^2$ when $\dot{\gamma} \to 0$ (see, for example, the monograph by Dealy$^9$). Consequently, the normal stress coefficient should, in principle, approach a constant value at very low shear rates. However, it is evident from Fig. 2 that both $\lim_{\text{Wi}\to 0} \hat{\Psi}_1(\text{Wi})$ and $\lim_{\text{De}\to 0} \hat{\Psi}_1(\text{De})$ diverge. These observations further motivated us to explore the relationship among $l, \beta, \gamma^*, \text{ and } m$ in order to determine the optimal power-law index, $l_{\text{opt}}$. Additionally, we seek to comprehend the underlying reasons for the diverging nature of the first normal stress coefficient curves when $\dot{\gamma} \to 0$ and $\omega \to 0$.

To gain insights into the limiting behavior of these functions, we evaluate the asymptotes of Eqs. (12) and (24) at low and high De and Wi. We start by considering the case when $\text{Wi}/\gamma^* \ll 1$. From Eq. (24), we can rewrite the dimensionless timescale, $\gamma$, as follows:

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 2. Comparison of the agreement between $\hat{\Psi}_1(\text{Wi})$ [Eq. (24)] and $\hat{\Psi}_1(\text{De})$ [Eq. (12)] plotted against the Weissenberg number (Wi) and Deborah number (De). Here, we used the following model parameters in panel (a) $\beta = 1/2$, $m = 2.0$, and $\gamma^* = \sqrt{15}/2$, and in panel (b) $\beta = 2/3$, $m = 2.0$, and $\gamma^* = \sqrt{15}/2$. For both figures, the ratio, $f$ [Eq. (25)], approaches a constant value at very low values of De and Wi, while it diverges at higher values of De and Wi.</figcaption></figure>
$$y = \frac{t - t'}{\tau_c} = \frac{\dot{y}_0 (t - t')}{\dot{y}_0 \tau_c} = \frac{\gamma(t', t)}{\text{Wi}}.$$

Therefore, using the following recurrence property of the MLF:³⁷

$$E_{a,b}(z) = \frac{1}{z} \left[ E_{a,b-a}(z) - \frac{1}{\Gamma(b-a)} \right], \quad (29)$$

where $\Gamma(b-a)$ is the gamma function, and considering the following asymptotic behavior of the MLF function for large arguments:⁴⁷

$$E_{a,b}(z) = - \sum_{k=1}^{N} \frac{z^{-k}}{\Gamma(b-ak)} + O(z^{-(N+1)}); \quad |z| \gg 1, \quad (30)$$

we can simplify Eq. (24) to obtain a new expression for the non-dimensional normal stress coefficient, $\hat{\Psi}_1(\text{Wi})$,

$$\hat{\Psi}_1(\text{Wi}) = - \sum_{k=1}^{N} \frac{\int_0^{\infty} y^{-k(1-\beta)} \frac{1}{1 + (\text{Wi}y/\gamma^*)^m} dy}{(-1)^k \Gamma(-1 - k(1 - \beta))}. \quad (31)$$

By making a variable substitution, $x = (\text{Wi}y/\gamma^*)$ and using Euler's reflection formula,³⁷ i.e.,

$$\int_0^\infty \frac{x^a}{1+x^m} dx = \frac{\pi}{m \sin\left[\pi\left(\frac{1+a}{m}\right)\right]}, \quad \text{for } 0 < \frac{1+a}{m} < 1, \quad (32)$$

we can further simplify Eq. (31) as follows:

$$\hat{\Psi}_1(\text{Wi}) \approx - \sum_{k=1}^{N} \frac{(\text{Wi}/\gamma^*)^{k(1-\beta)-1}}{(-1)^k \Gamma(-1-k(1-\beta))} \times \frac{\pi}{m \sin\left[\pi\left(\frac{1+q}{m}\right)\right]}, \quad (33)$$

where $q = -k(1 - \beta)$. It is important to note that the asymptotic expression we derive for the non-dimensional normal stress coefficient [Eq. (33)] is valid only when the following condition holds:

$$0 < (1 + q)/m < 1. \quad (34)$$

By following a similar procedure, we also derive an asymptotic expression of the non-dimensional normal stress coefficient given by Eq. (24) for the case of $\text{Wi}/\gamma^* \gg 1$. To do so, we utilize the following two properties of the MLF:²⁸,³⁷

$$E_{a,b}(z) = \frac{b}{\Gamma(1+b)} + z E_{a,b+a}(z), \quad (35)$$

and

$$E_{a,b}(z) = \sum_{k=0}^{N} \frac{z^k}{\Gamma(ak+b)} + O(z^{N+1}); \quad |z| \ll 1. \quad (36)$$

Therefore, for very high values of Wi, Eq. (24) can be approximated as

$$\hat{\Psi}_1(\text{Wi}) \approx \frac{\beta}{\Gamma(1-\beta)} \left(\frac{\text{Wi}}{\gamma^*}\right)^{\beta-2} \frac{\pi}{m \sin\left[\pi\left(\frac{2-\beta}{m}\right)\right]} + \sum_{k=0}^{\infty} \frac{(-1)^k (\text{Wi}/\gamma^*)^{-(p+1)}}{\Gamma((1-\beta)k + (1-2\beta))} \frac{\pi}{m \sin\left[\pi\left(\frac{1+p}{m}\right)\right]}, \quad (37)$$

where $p = (2+k)(1-\beta)$. Equation (37) is valid only if the following two inequalities are satisfied:

$$0 < \frac{2-\beta}{m} < 1, \quad \text{and} \quad 0 < \frac{1+p}{m} < 1. \quad (38)$$

Thus, in general, Eq. (33) and Eq. (37) are both hold true under the condition where

$$m > (2 - \beta).$$

We refer to this limit as the strong damping limit because the large value of $m$ results in the rapid sigmoidal-like decay of the damping function, $h(\gamma)$, at large strains $\gamma > \gamma^*$ [cf. Eq. (20)]. We also determine the asymptotic limits of the linear viscoelastic prediction, $\hat{\Psi}_1(\text{De})$, from Eq. (12). To summarize the asymptotic behavior of both $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$, we report the corresponding limits in Table I. When considering the asymptotic forms of $\hat{\Psi}_1(\text{Wi})$, we have focused solely on the first terms of the series solutions. This approach is justified because these terms dominate the behavior of $\hat{\Psi}_1(\text{Wi})$ for both high and low values of Wi.

In one early study, Larson²⁰ derived an analytical expression for the first normal-stress difference considering a power-law form of the relaxation modulus, i.e., $G(t) = ct^{-d}$, where $c$ has units of Pa·s$^d$. The expression is as follows:

$$N_1 = \frac{2c}{(2-d)\Gamma(d)} \tau_c^{2-d} \dot{\gamma}^2, \quad (39)$$

where $\tau_c$ corresponds to the cutoff relaxation time, calculated as $\tau_c = 1.7/\dot{\gamma}$ (refer to Eq. 23 in Larson's 1985 paper²⁰). Substituting this expression for $\tau_c$ in Eq. (39) and using Eq. (2) to compute the exponent of the primary normal stress coefficient, we obtain a scaling, $\hat{\Psi}_1 \sim \dot{\gamma}^{d-2}$, which is similar to the exponent we obtain for $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$ for $\text{Wi}/\gamma^* \gg 1$ and $\text{De} \gg 1$, i.e., $(\beta - 2)$ (cf. Table I). Notably, the relaxation modulus, described by Eq. (16), exhibits power-law-like behavior, $G(t) = \mathcal{G}t^{-\beta}$, when $t \gg \tau_c$. Therefore, the power-law decay constant $d$ in Larson's derivation corresponds to $\beta$ in our analysis.

To validate the accuracy of these approximate asymptotic expressions, we compare the analytical approximations (Table I) and the full numerical solutions of $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$ in Fig. 3 using the same model parameters as those in Fig. 2. The analytical solutions closely align with the full numerical solutions, particularly at extremely high and low values of Wi and De, as depicted in Figs. 3(a) and 3(b).

<table><caption>TABLE I. The asymptotes of $\hat{\Psi}_1(\text{Wi})$ and $\hat{\Psi}_1(\text{De})$.</caption><thead><tr><th></th><th>Ψ<sub>1</sub>(Wi)</th><th>Ψ<sub>1</sub>(De)</th></tr></thead><tbody><tr><td>Wi/γ* ≪ 1</td><td>(Wi/γ*)<sup>-β</sup> / (β-2) π / m sin(πβ/m)</td><td>2(De)<sup>-β</sup> cos(πβ/2)</td></tr><tr><td>De ≪ 1</td><td></td><td></td></tr><tr><td>Wi/γ* ≫ 1</td><td>β / (β-1) (Wi/γ*)<sup>(β-2)</sup></td><td>2(De)<sup>(β-2)</sup> cos(πβ/2) ×</td></tr><tr><td>De ≫ 1</td><td>π / m sin[π(2-β)/m]</td><td>[1 + cot<sup>2</sup>(πβ/2)]<sup>l</sup></td></tr></tbody></table>
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 3. Verifying the validity of the asymptotic limits of $\dot{\Psi}_1(\text{Wi})$ and $\dot{\Psi}_1(\text{De})$. In (a), the prediction of $\dot{\Psi}_1$ using the K-BKZ formulation [Eq. (24)] is represented by the thick red solid line. The green solid and blue dashed lines depict the first term and the first two terms of the series expression for the asymptotic form of $\Psi_1(\text{Wi})$ [Eq. (37)]. In (b), the red full line corresponds to the prediction of $\Psi_1$ using Laun's empirical formula with the FML model [Eq. (12)]. The green lines represent the asymptotic limits of $\Psi_1(\text{De})$ at low and high values of De (Table I).</figcaption></figure>

Moreover, from Table I and considering $\text{Wi} = \text{De}$, we can also determine the lower and upper limits of the ratio, $f$. When $\text{Wi} = \text{De} \ll 1$, the ratio, denoted by $f_0$ for clarity, becomes

$$f_0 = \frac{\lim_{\text{Wi}\to 0} \dot{\Psi}_1(\text{Wi})}{\lim_{\text{De}\to 0} \dot{\Psi}_1(\text{De})} = \frac{\pi \gamma^{*\beta}}{2 m \Gamma(\beta - 2) \sin\left(\frac{\pi \beta}{m}\right) \cos\left(\frac{\pi \beta}{2}\right)} \quad (40)$$

Careful consideration of Eq. (40) shows that the ratio $f_0$ does not vary with $\text{Wi}$, $\text{De}$, nor even with Laun's exponent $l$. This ratio is, however, a function of $\beta$, $m$, and $\gamma^*$, suggesting that, at low deformations, $f_0$ approaches a constant that is not precisely unity as evident from the insets of Fig. 2. For $\text{Wi} = \text{De} \gg 1$, the ratio $f_\infty$ can be expressed as follows:

$$f_\infty = \frac{\lim_{\text{Wi}\to\infty} \dot{\Psi}_1(\text{Wi})}{\lim_{\text{De}\to\infty} \dot{\Psi}_1(\text{De})} = \frac{\beta \pi [1 + \cot^2(\pi \beta / 2)]^{-l}}{2 \Gamma(1 - \beta) \gamma^{*(\beta-2)} m \cos\left(\frac{\pi \beta}{2}\right) \sin\left[\pi \left(\frac{2 - \beta}{m}\right)\right]} \quad (41)$$

In contrast to $f_0$, the ratio $f_\infty$ is a function of $\beta$, $m$, and $\gamma^*$ as well as $l$.

As was evident from Fig. 2, both $\dot{\Psi}_1(\text{Wi})$ and $\dot{\Psi}_1(\text{De})$ diverge as $\text{Wi} = \text{De} \to 0$; i.e., our calculations show that there is no well-defined value of a limiting first normal stress coefficient $\Psi_{1,0}$. Now, from the asymptotic limits derived when $\text{Wi} = \text{De} \ll 1$ (see Table I), the underlying cause for this, initially surprising, conclusion becomes apparent. Considering fixed values of $\beta$, $m$, and $\gamma^*$, we can see that the limits $\lim_{\text{Wi}\to 0} \dot{\Psi}_1(\text{Wi})$ as well as $\lim_{\text{De}\to 0} \dot{\Psi}_1(\text{De})$ vary as $\dot{\gamma}^{-\beta}$ and $\omega^{-\beta}$, respectively. In other words, the power-law exponent characterizing the variation in the normal stress coefficient at low deformation rates is

$$\beta_0 = - \frac{d(\ln \dot{\Psi}_1(\text{Wi}))}{d \ln(\text{Wi})} = - \frac{d(\ln \dot{\Psi}_1(\text{De}))}{d \ln(\text{De})} \quad (42)$$

The key conclusion obtained from this discussion is that for materials described by the fractional Maxwell model, there is no unique defined value of the primary normal stress coefficient in the limit of low deformation rates. The first normal stress difference decays at low shear rates as $N_1 \sim \dot{\gamma}_0^{2-\beta}$, and the corresponding material property diverges as $\Psi(\dot{\gamma}) \sim \dot{\gamma}_0^{-\beta}$. To further support this assertion, we can separately derive an expression for $\Psi_{1,0}$ in terms of the underlying relaxation spectrum $H(\tau)$ captured by the fractional derivative. This is given by the following equation adapted from Chap. 3 of Dealy and Wissbrun:⁹

$$\Psi_{1,0} = 2 \int_{0}^{\infty} G(s) s \, ds, \quad (43)$$

where $s = t - t'$.

Alternatively, when expressed in terms of a continuous relaxation distribution, Eq. (43) can be written as⁹

$$\Psi_{1,0} = 2 \int_{0}^{\infty} H(\tau) \, \tau \, d\tau. \quad (44)$$

Moreover, Eq. (43) can also be written in terms of measurable linear viscoelastic properties as follows:⁹

$$\Psi_{1,0} = 2\eta_0^2 J_s^0. \quad (45)$$

Here, $\eta_0$ is the zero shear viscosity, and $J_s^0$ represents the steady-state compliance.⁹ The recent study by Rathinaraj and coauthors²⁹ mathematically proved that the fractional Maxwell liquid (FML) model does recover a bounded $\eta_0$ at very low shear rates in contrast to the full fractional Maxwell model (FMM). Consequently, the diverging behavior of $\Psi_{1,0}$ must be the same as the divergent nature of $J_s^0$.
The analytical expression for the creep compliance $J(t-t')$ of the four parameter FMM is given as follows:²⁸

$$J(t-t') = \frac{(t-t')^{\alpha}}{\mathbb{V}\Gamma(\alpha+1)} + \frac{(t-t')^{\beta}}{\mathbb{G}\Gamma(\beta+1)}. \quad (46)$$

By substituting $s = t - t'$, $\alpha = 1$ (for the FML model), and $\mathbb{V} = \eta_0$, we derive the following expression for the creep compliance of the FML model:

$$J(s) = \tau_c/\eta_0 \left[ (s/\tau_c) + (s/\tau_c)^\beta / \Gamma(1+\beta) \right]. \quad (47)$$

The steady-state creep compliance is defined by Dealy and Wissbrun⁹

$$J_s^0 = J(s) - s/\eta_0. \quad (48)$$

By combining Eq. (48) with Eq. (47) and making the substitution $s/\tau_c = \frac{2\eta_0}{\mathbb{W}i}$ we can show that

$$\hat{J}_s^0 = \frac{J_s^0 \eta_0}{\tau_c} \approx \left( \frac{\mathbb{W}i}{\dot{\gamma}_0 s} \right)^{-\beta} \frac{1}{\Gamma(1+\beta)}. \quad (49)$$

Thus, in the FML model, unless $\beta = 0$, the limiting value of $J_s^0$ (and thus also $\Psi_{1,0}$) diverges instead of approaching a constant value. It is worth noting that when we also substitute $\beta = 0$, which corresponds to the case of a single-mode Maxwell relaxation kernel, we obtain a constant value for $J_s^0$ and, consequently, a constant value for $\Psi_{1,0}$. This limit is consistent with Eq. (27) and $f_0 = 1$.

Furthermore, when considering the asymptotic limit of $\hat{\Psi}_1(\text{De})$ in the case of $\mathbb{W}i = \text{De} \ll 1$ (Table I), it is observed that the function varies as $(\text{De})^{-\beta}$. This observation is supported by the relationship obtained from the simple second-order fluid theory¹⁰,²⁴

$$\lim_{\dot{\gamma} \to 0} \Psi_1 = \lim_{\omega \to 0} \frac{2G'(\omega)}{\omega^2}. \quad (50)$$

Notably, Friedrich and Braun also stated in their 1992 paper⁴⁸ that “the coefficient of the first normal stress difference, $\Psi_{1,0}$ does not exist.”

### III. OPTIMIZATION OF THE POWER-LAW INDEX ($l_{opt}$)

In the present formulation, to satisfy Laun's rule, we wish the ratio $f$ defined in Eq. (25) to be unity across a wide range of deformation rates. Our fractional K-BKZ model contains five material constants ($\beta, \mathbb{G}, \eta_0$, and $\gamma^*$, $m$) resulting in a high dimensional optimization problem. We thus first consider various limiting cases. First, in order to achieve $f_0 = 1$ from Eq. (40), we have to find the optimal values of the model parameters ($\beta, m, \gamma^*$) characterizing a given complex fluid. Once these optimal values have been identified, we can proceed to compute an optimal value of Laun's exponent, $l_{opt}(\beta, m, \gamma^*)$, to ensure $f \approx 1$ across a wide range of deformation rates. The procedure for determining the optimal parameters is explained in detail in the subsequent section.

Recognizing that $\beta$ is set by the slope of the linear viscoelastic relaxation spectrum and that reptation theory for entangled polymer melts suggests $m \approx 2$, we initiate our analysis by calculating the needed value of the critical strain, denoted as $\gamma_{opt}^*$, in order to force the ratio $f_0$ to be exactly equal to unity while holding fixed values of $\beta$ and $m$. From Eq. (40), we can determine that for $\beta = 1/2$ and $m = 2.0$, the optimal critical strain should be $\gamma_{opt}^* = 2.264$ to achieve $f_0 = 1$. Similarly, for $\beta = 2/3$ and $m = 2.0$, the optimal critical strain was determined to be $\gamma_{opt}^* = 2.166$. These are quite representative values for a range of typical polymer melts.

To understand the consequences of these choices for $\gamma_{opt}^*$, we generate plots of the ratio, $f$, as a function of $\mathbb{W}i$ and $\text{De}$, as shown in Fig. 4. From Figs. 4(a) and 4(b), we can see that the curves corresponding to both $l=0$ and $l=0.7$ converge to unity at very low

<figure><img src="image_18.png" /><figcaption>FIG. 4. Ratiometric plot of $f$ using optimal values of the critical strain, $\gamma_{opt}^*$. Model parameters used in panel (a) $\beta = 1/2$, $m = 2.0$, and $\gamma_{opt}^* = 2.264$, and in panel (b) $\beta = 2/3$, $m = 2.0$, and $\gamma_{opt}^* = 2.166$. For both figures, the ratio converges to unity at sufficiently low values of $\mathbb{W}i$ and $\text{De}$. However, for both figures, the red lines, corresponding to $l=0.7$, are observed to be closer to the target value of unity indicated by the $f_{target} = 1$ line (represented by broken blue line). This shows that Laun's suggested value of $l=0.7$ is significantly more accurate than $l=0$.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 5. Contour plots of the ratio $f(Wi = De, l)$ for predicting the optimal power-law index to be used in Laun's rule. (a) Model parameters used: $\beta = 1/2$, $m = 2.0$, and $\gamma^* = 2.264$. The optimal power-law index $l_{opt}$, which results in a ratio of $f \approx 1$, falls within the range of 0.54-0.61. (b) Model parameters used: $\beta = 2/3$, $m = 2.0$, and $\gamma^* = 2.166$; the optimal power-law index $l_{opt}$ is found to be in the range of 0.71-0.92, which is greater than the value $l = 0.7$ originally proposed by Laun.</figcaption></figure>

values of Wi and De and this is thus for any choice of $l$ as expected from Eq. (40). It is also clear that Laun's proposed value of $l = 0.7$ is considerably more accurate in ensuring that the ratio $f = \Psi_1(Wi)/\Psi_1(De)$ is close to unity across a wide range of deformation rates. However, what is not clear from Fig. 4 is if this is truly an optimal choice for the exponent $l$.

In order to examine the variation of $f$ with respect to Wi and De for different $l$ values, we generate contour plots of the ratio $f$ [Eq. (25)] as a function of $l$, and Wi or De (Fig. 5). We also indicate specific contour values at which $f = 0.98$ and $f = 1.02$ by dashed and dashed lines, respectively (i.e., indicating that within these boundary lines, the error is within $\pm 2\%$). From Fig. 5(a), we can see that in order to achieve a ratio of $f \approx 1$, $l_{opt}$ should in fact be chosen within the range of 0.54-0.61 for model parameters ($\beta, m, \gamma^*_{opt}$) = (1/2, 2.0, 2.264). Similarly, for $\beta = 2/3, m = 2$, and $\gamma^*_{opt} = 2.166$, $l_{opt}$ should fall within the range of 0.71-0.92. The resulting accuracy of Laun's rule is within 2% for both cases. This is typically much higher accuracy than what is attainable with experimental measurements of $N_1(\gamma)$.

Finally, to calculate the exact value of $l_{opt}$, we used the following approach. For each value of $l$, we calculate the root-mean-square error of the ratio $f$ [Eq. (25)] and the desired value of the ratio using the following formula:

$$ \varepsilon_l = \sqrt{\frac{\sum_{i=1}^{N} (f_i - f_{\text{target}})^2}{N}}, \quad (51) $$

where $f_{target}$ denotes the desired value of the ratio, i.e., unity, and $N$ represents the number of data points spanning the range $Wi_{min}(i = 1) = 10^{-2}$ to $Wi_{max}(i = N) = 10^{+2}$. The value of the power-law index that yields the minimal error is considered the optimal value, $l_{opt}$. For ($\beta, m, \gamma^*_{opt}$) = (1/2, 2.0, 2.264), we find $l_{opt} \approx 0.59$, and for $\beta = 2/3, m = 2$, and $\gamma^*_{opt} = 2.166$, we obtain $l_{opt} \approx 0.82$.

Importantly, we can now conclude that the optimal value of the Laun exponent is not constant and also not constrained to the range $0 \le l \le 0.7$, as initially assumed in the study. As we show below, the value of $l$ can take a range of values depending on the material characteristics.

Up to this point, we have verified our fractional viscoelastic framework for obtaining $l_{opt}$ considering only two specific values of $\beta$

<figure><img src="image_20.png" /><figcaption>FIG. 6. Plot of the optimal power-law index, $l_{opt}$ for varying $\beta$ and $\gamma^*$ with a fixed damping exponent value $m = 2.0$. The symbols represent the optimal value obtained by computing the ratio $f$ [Eq. (25)] over a range of $l$ values and selecting the optimal value, $l_{opt}$, that minimizes the root-mean-square error [Eq. (51)]. Additionally, the dash-dot lines on the graph represent the predicted values of $l_{opt}$ using the asymptotic limits for high values of Wi and De, as given by Eq. (52). Finally, the solid line corresponds to Laun's suggested value $l = 0.7$, which is exactly satisfied only for specific combinations of $\beta, m$, and $\gamma^*$ (represented by the pink inverted triangles).</figcaption></figure>
consistent with the Rouse ($\beta = 1/2$) and Zimm ($\beta = 2/3$) relaxation spectra, respectively. The next objective is to explore how sensitive the exponent $l$ is to different values of $\beta$. To achieve this, we apply the aforementioned optimization technique to calculate $l_{opt}$ across a range of $\beta$ values spanning from $0.1 \le \beta \le 0.85$, which corresponds to the practical range commonly encountered in numerous complex fluids. The resulting curves are shown in Fig. 6.

From Fig. 6, we can reconfirm that Laun's proposed value of the power-law index, $l=0.7$, is only optimal for specific combinations of $\beta$, $m$, and $\gamma^*$. For example, for materials whose strain-softening effect is characterized by a critical strain, $\gamma^* = 2.50$, and damping exponent, $m=2.0$ (represented by the pink inverted triangle in Fig. 6), $l=0.7$ is only appropriate if $\beta = 0.45$. Another observation from Fig. 6 is that predicting $l_{opt}$ using a damping function consistent with the Doi-Edwards model parameters ($m=2.0, \gamma^* = \sqrt{15}/2$) consistently yields values lower than Laun's suggested power-law index, regardless of the value of $\beta$. Moreover, it is worth noting that it is also possible to obtain negative values of $l_{opt}$ based on the specific ranges of parameter values!

In our analysis, we optimized the values of $\gamma^*$ while maintaining constant $\beta$ and $m$ values in order to derive the optimal exponent $l$ that most accurately predicts the value and the scaling of the primary normal stress coefficient. It is important to recognize that these values ($\beta$, $m$, and $\gamma^*$) vary from material to material but they can be determined from the results of SAOS and step-strain tests that are typically easier to perform (and do not require a rheometer equipped with a normal force transducer) (refer to Sec. II for more details).

Finally, we attempt to establish a simple mathematical relationship between $l_{opt}$ and these three experimentally accessible constitutive parameters ($\beta$, $m$, and $\gamma^*$). By rearranging Eq. (41), and setting $f_{\infty} = 1$, we derive the following equation:

$$l_{\text{opt}} = \frac{\ln \left( \pi \beta (\gamma^*)^{(2-\beta)} \right) - \ln \left( 2m \Gamma(1-\beta) \sin \left( \frac{\pi(2-\beta)}{m} \right) \cos \left( \frac{\pi\beta}{2} \right) \right)}{\ln \left( 1 + \cot^2 \left( \frac{\pi\beta}{2} \right) \right)} \quad (52)$$

This constraint [Eq. (52)] is applicable for materials that demonstrate sufficient strain-softening in step strain experiments. We refer to the constraint $m > (2 - \beta)$ as the “strong damping” limit, whereas Rathinaraj et al.²⁹ found that when considering the Cox-Merz rule, the constraint $m > (1 - \beta)$ was sufficient to ensure strong damping in a general complex fluid. The new condition is clearly more restrictive, but the condition $m > (2 - \beta)$ is still satisfied for a wide range of complex fluids.⁴⁴ We find that Laun's relationship is not applicable (and was never intended to be) in cases of weak damping corresponding to, for example, gel networks or physically associated materials such as hydrogels, dough, etc. Evaluating Eq. (52) for fixed values of $\beta$ and $m$, it is found that $l_{opt}$ increases monotonically with increasing $\gamma^*$. On the contrary, for fixed values of $\beta$ and $\gamma^*$, $l_{opt}$ decreases with increasing values of $m$. To facilitate comparison with direct numerical simulations, we include the predictions of Eq. (52) as dashed lines in Fig. 6. It is clear that this analytic expression closely agrees with the values obtained through numerical simulations for $\beta$ values up to at least $\beta = 0.66$.

Finally, to showcase the broad applicability of our analysis, we now validate this framework by reconsidering the experimental data discussed in the introduction (Sec. I).

IV. REPRESENTATIVE EXPERIMENTAL EXAMPLE

To evaluate the accuracy of our approach for predicting $l_{opt}$, we employ the same representative samples initially used in Sec. I to illustrate the accuracy and limitations of Laun's correlation. Specifically, we consider LDPE Melt I and a 0.5 wt. % solution of Xanthan gum in water.

In Figs. 7(a) and 7(b), we present the linear viscoelastic relaxation modulus, $G(t)$, as well as the damping values $h(\gamma)$ reproduced from Laun's 1986 paper.¹ Subsequently, we apply Eq. (16) to fit the measured data and determine the best-fit values of the three parameters $\mathcal{G}$, $\beta$, and $\eta_0$. The fitted values are as follows: $\mathcal{G} = 9.60 \times 10^3 \text{ Pa} \cdot \text{s}^\beta$, $\beta = 0.47$, and $\eta_0 = 5.8 \times 10^4 \text{ Pa} \cdot \text{s}$. The parameter $\beta$ describes the short-time relaxation data, $G(t) \sim t^{-\beta}$ or, correspondingly, the high-frequency dependence of complex modulus data, $G^* \sim \omega^\beta$.

Furthermore, we fit the measured damping function values using Eq. (20) and obtain the fitted parameters $\gamma^* = 4.10$ and $m = 1.83$. To fit both measured relaxation modulus, $G(t)$, data and damping function, $h(\gamma)$, values, we employ a constrained nonlinear optimization toolbox “fmincon” within MATLAB. We chose the “interior-point” algorithm to obtain optimal values of the fitting parameters. It is worth noting that the material properties of LDPE Melt I satisfy the strong damping limit, i.e., $m > (2 - \beta)$, enabling us to apply our framework to predict $l_{opt}$. However, the critical strain is much larger than the Doi-Edwards value, $\gamma^* = \sqrt{15}/2$, most probably because of the polydisperse nature of the melt and the presence of long chain branches.

Using the obtained values of $\beta$, $m$, and $\gamma^*$, we calculate the optimal value of $l$ for LDPE Melt I, which is determined to be $l_{opt} = 2.23$. The resulting curves corresponding to $l=0, l=0.7$, and $l_{opt} = 2.23$ are presented in Fig. 7(c). It is evident that the prediction of the primary normal stress coefficient $\Psi_1$ from the linear viscoelastic data using the power-law index $l_{opt} = 2.23$ exhibits superior agreement with the measured data compared to the cases of $l=0$ and $l=0.7$.

For a more quantitative demonstration of this, we also present the non-dimensionalized mean square error in the inset of the same figure [Fig. 7(c)], which is calculated using the following equation:

$$\hat{e}_m = \sum_{i=1}^{N} \left( 1 - \frac{\Psi_{1,i}^{(p)}(\dot{\gamma}_i)}{\Psi_{1,i}^{(m)}(\dot{\gamma}_i)} \right)^2 \quad (53)$$

where $\Psi_{1,i}^{(p)}(\dot{\gamma}_i)$ and $\Psi_{1,i}^{(m)}(\dot{\gamma}_i)$ indicate the predicted and measured values of the primary normal stress coefficient, respectively, at each value of shear rate. The reduction in the mean square error value calculated using Eq. (53) for $l_{opt} = 2.23$ further confirms the improved prediction of the first normal stress coefficient that can be obtained using our framework.

To validate our framework for the 0.5 wt. % Xanthan gum solution, we employ a similar procedure. We reproduce in Fig. 8(a) the $G(t)$ data from the work of Jaishankar and McKinley.²⁷ By fitting the measured data using the FML model in Eq. (16), we determine the relaxation exponent to be $\beta = 0.23$. The values of $\gamma^* = 1.83$ and $m = 2.0$ are taken directly from Fig. 6 of the same paper.²⁸ These values again satisfy our strong damping constraint, i.e., $m > (2 - \beta)$. Utilizing our optimization procedure with these $\beta$, $m$, and $\gamma^*$ values, we find that because of the much lower value of the critical strain $\gamma^*$, the optimal value for the exponent $l$ is reduced to $l_{opt} = 0.134$. Subsequently, in Fig. 8(b), we present the reported experimental data for $\Psi_1(\dot{\gamma})$ alongside the corresponding curves obtained by employing Laun's rule for the power-law index values $l=0, l_{opt} = 0.134$, and
<figure><img src="image_21.png" /><figcaption>FIG. 7. (a) Linear viscoelastic relaxation modulus data for LDPE Melt I, reproduced with permission from H. M. Laun, J. Rheol. 30, 459–501 (1986). Copyright 1986 Society of Rheology.¹ The solid line on the plot represents the fit obtained using the three-parameter Fractional Maxwell Liquid (FML) model with the following values: $G = 9.60 \times 10^3 \text{ Pa} \cdot \text{s}^\beta$, $\beta = 0.47$, and $\eta_0 = 5.8 \times 10^4 \text{ Pa} \cdot \text{s}$. (b) Measured damping values are well described by the Soskey–Winter damping function with $\gamma^* = 4.10$ and $m = 1.83$. (c) Experimental values of the primary normal stress coefficient plotted against the shear rate. The original data points represented by large purple filled circles correspond to the experimental results (see Laun'). The hollow triangular green symbols correspond to predicted values using Laun's empirical relationship with $l=0$, based on the reported discrete $G'(\omega)$ data.¹ The blue rectangular symbols represent predictions using Laun's empirical relationship with $l=0.7$, utilizing both $G'(\omega)$ and $G''(\omega)$ data reported by Laun.¹ The red solid line on the plot represents the predicted variation in $\Psi_1(\dot{\gamma})$ obtained using the optimal power-law index value, $l_{\text{opt}}(\beta, m, \gamma^*) = 2.23$ and Eq. (12).</figcaption></figure>

$l=0.7$. Notably, again for this 0.5 wt. % Xanthan gum solution, there is excellent agreement between the measured normal stress data and the prediction for $\Psi_1(\dot{\gamma})$, utilizing our optimized value of $l_{\text{opt}}$ and Eq. (7).

V. CONCLUSION

We have assessed the accuracy of Laun's empirical relationship for predicting the first normal stress coefficient of complex materials. Although the empirical formula [Eq. (4)] successfully predicts the

general shape of the first normal stress coefficient curve, there is a noticeable discrepancy or shift between the predicted values and the actual observed data. It is important to note that this difference is small for the specific entangled and polydisperse polymer melt considered by Laun,¹ i.e., LDPE Melt I. However, it can be more significant for other complex fluids like the concentrated semiflexible polymer solution studied by Jaishankar and McKinley,²⁸ or the systems studied by Tang *et al.*⁹
<figure><img src="image_22.png" /><figcaption>FIG. 8. (a) Linear viscoelastic relaxation modulus data for a 0.5 wt. % Xanthan gum solution, reproduced with permission from Jaishankar and McKinley, J. Rheol. 58, 1751–1788 (2014). Copyright Society of Rheology.²⁸ The solid line on the plot represents the fit using the Fractional Maxwell (FML) model with the following parameters: $G = 6.188 \text{ Pa} \cdot s^{\beta}$, $\beta = 0.225$, and $\eta_0 = 6.92 \times 10^2 \text{ Pa} \cdot \text{s}$. (b) Experimental values of the primary normal stress coefficient plotted against the shear rate. The data points represented by the purple filled circles correspond to the experimental results (see Jaishankar and McKinley²⁸). The hollow triangular green symbols correspond to predicted values using Laun’s empirical relationship with $l=0$, based on the discrete $G'(\omega)$ data.²⁸ The blue rectangular symbols represent predictions using Laun’s empirical relationship with $l=0.7$, utilizing both $G'(\omega)$ and $G''(\omega)$ data.²⁸ The red solid line represents predicted variation in $\Psi_1(\dot{\gamma})$ data obtained using the optimal power-law index value, $l_{opt}(\beta, m, \gamma^*) = 0.134$ and Eq. (12).</figcaption></figure>

To gain a mathematical understanding of the specific conditions under which Laun’s relationship holds true for a complex fluid, we considered the fractional Maxwell liquid (FML) model of linear visco-elasticity to express Laun’s empirical relationship compactly in terms of three linear viscoelastic model parameters $\eta_0$, $G$, and $\beta$ [Eq. (12)]. Next, we also derived an integral mathematical formulation for the primary normal stress coefficient, which is given by Eq. (24). This formulation is based on the time-strain separable Wagner model in conjunction with the fractional relaxation modulus predicted by the FML model, which is of the Mittag–Leffler form. We also incorporated a generalized Soskey–Winter form of the damping function $h(\gamma)$ to describe the material’s nonlinear response in large step strains in order to evaluate the resulting integral. For clarity, we denoted the two competing expressions for the non-dimensional primary normal stress coefficient as $\dot{\Psi}_1(\text{De})$ and $\dot{\Psi}_1(\text{Wi})$, respectively.

By investigating the limiting cases of $\dot{\Psi}_1(\text{De})$ and $\dot{\Psi}_1(\text{Wi})$, we have shown that when the exponent $\beta$ is set to zero (corresponding to the single-mode Maxwell model), Laun’s relationship does not hold. Furthermore, materials with a narrow relaxation spectrum ($\beta \to 0$) and/or no strain-softening effect, represented by a critical strain of $\gamma^* \to \infty$, do not conform to Laun’s relationship.

However, for materials exhibiting broad relaxation spectra and significant strain softening ($m > 2 - \beta$), upon comparing Eqs. (12) and (24), we determined that the power-law index proposed by Laun, denoted as $l=0.7$ (pink line in Fig. 6), is not universally constant. It is clear from Fig. 4 that Laun’s rule provides a marked improvement from the simplest approach of setting $l=0$ in Eq. (7), which is strictly valid only in the limit of low frequencies and gives

$$N_1(\dot{\gamma}) \approx 2G'(\omega)|_{\dot{\gamma}=\omega}.$$

Detailed numerical computations show that the optimal choice of exponent $l$ is a function of the FML model parameter $\beta$ (which characterizes the power-law scaling of the viscoelastic relaxation spectrum at high frequencies and short times) as well as the nonlinear strain-damping parameters $m$ and $\gamma^*$. Consequently, it is possible to improve the predictive agreement between Laun’s empirical formula and the observed data by appropriately selecting an optimal value of the exponent $l$.

We have presented a framework for calculating the appropriate optimal value of the exponent, denoted as $l_{opt}$. Based on our analysis, we have established a constraint that distinguishes materials with weak damping from those with strong damping. Our analysis reveals that if the damping parameter $m > (2 - \beta)$, the material can be considered to exhibit strong strain softening, making Laun’s rule applicable to such materials. Conversely, when $m$ is less than or equal to $(2 - \beta)$, Laun’s rule is not applicable.

Furthermore, our asymptotic analysis of $\dot{\Psi}_1(\text{De})$ reaffirms the low shear rate behavior of these polydisperse complex fluids. Specifically, we find that when $\omega \ll \tau_c^{-1}$, $\dot{\Psi}_1(\text{De})$ remains unaffected by the power-law index $l$ and diverges as $\text{De}^{-\beta}$. However, as the frequency $\omega$ and shear rate $\dot{\gamma}$ increase beyond $\tau_c^{-1}$, we notice that the shape of the curves for the primary normal stress coefficient becomes dependent upon the specific value of $l$.

In order to select the appropriate optimal value for $l$, it is necessary to know the values of the parameters $\beta$, $m$, and $\gamma^*$. The exponent $\beta$ can be calculated by measuring the slope of the $G'(\omega)$ or $G''(\omega)$ curves at very high frequencies or the power-law decay of the relaxation modulus at short times. In addition to SAOS data, another important component that needs to be determined within this framework is the form of the strain-dependent damping function. The exponent $\beta$
can also be calculated by fitting the linear viscoelastic relaxation modulus data using Eq. (16). Therefore, when step-strain test data are available over a range of strain amplitudes, it becomes possible to calculate all three required model parameters. Once the values of $\beta, m$, and $\gamma^*$ are obtained, we have outlined the following two approaches to calculate $l_{opt}$:

* The first approach involves calculating the ratio ($f$) between $\Psi_1(\text{Wi})$ and $\Psi_1(\text{De})$ [cf. Eq. (25)] for varying values of $l$. Subsequently, the root-mean-square error given by Eq. (51) is computed. The optimal value $l_{opt}$ then corresponds to the one that yields minimal error.

* The second approach is based on the asymptotic limits for $\Psi_1(\text{Wi})$ and $\Psi_1(\text{De})$ given in Table I at high values of Wi and De. We derive an explicit formula, given in Eq. (52), to calculate $l_{opt}$ based on the values of $\beta, m$, and $\gamma^*$.

The comparison of these two approaches (cf. Fig. 6) reveals that up to values of $\beta \le 0.66$, the approximation for $l_{opt}$ given by Eq. (52) closely aligns with the values obtained from the first approach, which involves numerical calculations.

Performing step-strain tests over a range of strain amplitudes can often be challenging or tedious. However, if one has access to small amplitude oscillatory shear (SAOS) experimental data, one can easily calculate the value of $\beta$ from the slope of the frequency-dependent complex modulus, $G^*(\omega)$. For the parameters $m$ and $\gamma^*$, it is possible to use as an initial approximation the values consistent with the Doi-Edwards theory, where $m=2$ and $\gamma^* = \sqrt{15}/2$. Subsequently, an estimate for $l_{opt}$ can be calculated using Fig. 6 and/or Eq. (52) as illustrated in Sec. III. It is important to note that this alternative framework will only yield accurate predictions for $l_{opt}$ if the material exhibits time-strain separability and the critical strain ($\gamma^*$) is reasonably close to the Doi-Edwards limit. Given the difficulties in making accurate measurements of $N_1(\dot{\gamma})$, even this approximate approach can yield substantially improved a priori estimates for the primary normal stress difference in complex fluids. Once this functional form is known, it then can be used in many of the engineering applications outlined in Tanner's groundbreaking book.⁴

Finally, we note that our proposed framework for determining $N_1(\dot{\gamma})$ from the available linear viscoelastic data is restricted to materials that exhibit time-strain separability. Thus, a remaining area for future research involves developing methods to predict the optimal power-law index of the functional form of the first normal stress difference $N_1(\dot{\gamma})$ for materials that do not exhibit time-strain separability.

## ACKNOWLEDGMENTS

The authors would like to acknowledge Dr. Bavand Keshavarz for a number of helpful discussions during the preparation of this work.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Mohua Das:** Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); & Writing – review & editing (equal). **Joshua David John**

**Rathinaraj:** Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); & Writing – review & editing (equal). **Liviu Iulian Palade:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); & Writing – review & editing (equal). **Gareth H. McKinley:** Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); & Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## APPENDIX: LAUN'S RULE FOR SINGLE-MODE RELAXATION MATERIALS AND MATERIALS WITH WEAK DAMPING

When $\beta=0$, the constitutive equation of the FML model, represented by Eq. (8), simplifies to the following equation for the single-mode Maxwell model:

$$ \sigma(t) + \frac{\eta_0}{G_0} \frac{d\sigma(t)}{dt} = \eta_0 \frac{d\gamma}{dt}. \quad (A1) $$

Here, $G_0$ represents the elastic modulus of the material and $\tau_0$ represents the single characteristic relaxation time, $\tau_0 = \eta_0/G_0$. The equations for the storage and loss moduli of the single-mode Maxwell model are given by

$$ G'(\omega) = G_0 \frac{(\omega\tau_0)^2}{1 + (\omega\tau_0)^2}, \quad (A2) $$

$$ G''(\omega) = G_0 \frac{\omega\tau_0}{1 + (\omega\tau_0)^2}. \quad (A3) $$

By substituting Eqs. (A2) and (A3) in Eq. (7), we can write the expression for the non-dimensionalized primary normal stress coefficient as follows:

$$ \hat{\Psi}_1(\text{De}) = \frac{\Psi_1(\dot{\gamma})}{2\eta_0\tau_0} = \frac{[1 + (\text{De})^2]^l}{[1 + (\text{De})^2]^l}. \quad (A4) $$

The limiting behaviors of Eq. (A4) can be described as follows:

* For $\text{De} \ll 1$, $\hat{\Psi}_1(\text{De}) \approx 1$.

* For $\text{De} \gg 1$, $\hat{\Psi}_1(\text{De}) \approx (\text{De})^{2(l-1)}$.

Next, we need to consider the normal stress differences predicted by the single-mode model (i.e., the left hand side of Laun's rule). We thus consider the general expression given by Eq. (22). Notably, $\beta=0$ and $\gamma^* \to \infty$ (i.e., $h(\gamma) = 1, \forall \text{Wi}$) is the case of the upper-convected Maxwell (UCM) model. When $\beta=0$ in Eq. (16), the Mittag-Leffler relaxation kernel reduces to the Maxwell-Debye form, $\exp(-t/\tau)$. Consequently, for the UCM model, we can represent the non-dimensionalized primary normal stress coefficient as

$$ \hat{\Psi}_1(\text{Wi}) = \frac{\Psi_1(\dot{\gamma})}{2\eta_0\tau_0} = \frac{1}{2} \int_0^\infty y^2 \exp(-y) dy = 1, \quad (A5) $$
$$ \therefore \Psi_1 = 2\eta_0\tau_0 = 2G_0\tau_0^2. \qquad (A6) $$

Equation (A6) represents the same analytical result expected for a simple fluid with fading memory [Eq. (28)]. Comparing Eqs. (A4) and (A5), we can say that, for the case of the UCM model, Laun's rule holds true when $\dot{\gamma} = \omega \ll (1/\tau_c)$ because in this limit both $\dot{\Psi}_1(\text{De})$ and $\dot{\Psi}_1(\text{Wi})$ approach unity. However, for high deformation rates and frequencies, Laun's rule is not valid except for the special case of $l=1$.

The limit $\beta=0$ and a finite value of the critical strain, $\gamma^*$ represent the Wagner model with a Maxwell kernel. For this, we can express the equation for $\Psi_1(\text{Wi})$ by using the Wagner integral formulation with a single mode exponential relaxation kernel and the Soskey-Winter damping function as follows:

$$ \dot{\Psi}_1(\text{Wi}) = \frac{\Psi_1(\dot{\gamma})}{2\eta_0\tau_0} = \int_0^\infty \frac{\gamma^2 \exp(-y)}{2[1 + (\text{Wi}y/\gamma^*)^m]} dy. \quad (A7) $$

The asymptotic limits of Eq. (A7) for very low and very high values of Wi are given as follows:

* When $\text{Wi} \ll 1$, $\dot{\Psi}_1(\text{Wi})$ approaches unity.

* When $\text{Wi} \gg 1$, $\frac{1}{1+(\text{Wi}y/\gamma^*)^m} \approx (\text{Wi}y/\gamma^*)^{-m}$

$$ \therefore \dot{\Psi}_1(\text{Wi}) \approx (1/2)(\text{Wi}/\gamma^*)^{-m} \Gamma(3-m). \quad (A8) $$

From these asymptotic analyses, we can conclude that for a material with a very narrow relaxation spectrum corresponding to a single relaxation time, $\tau_0$, the primary normal stress coefficient converges to a constant value ($= 2\eta_0\tau_0$) for very low values of De and Wi, i.e., when $\dot{\gamma}(=\omega) \to 0$. This finding aligns with the conclusion presented in Sec. II, which suggests that the single-mode Maxwell model with strain-damping can accurately describe $\Psi_{1,0}$. However, when comparing the asymptotic behavior of $\dot{\Psi}_1(\text{De})$ and $\dot{\Psi}_1(\text{Wi})$ for high values of De and Wi, it becomes apparent that the slopes of the two expressions in general differ. Thus, Laun's rule is not valid when $\text{De} \gg 1$ and $\text{Wi} \gg 1$ except for one special case when the following two relationships are satisfied:

$$ l = 1 - \frac{m}{2}, \quad \text{and} \quad \gamma^* = \left[ \frac{2}{\Gamma(3-m)} \right]^{1/m}. \qquad (A9) $$

Finally, we note that, in the case of weak damping (for which $m \le (2-\beta)$), the integral in Eq. (24) diverges. As a result, Laun's relationship is not applicable for physically crosslinked materials such as metal-crosslinked hydrogel networks or highly viscoelastic gluten gels and bread doughs with weak damping (see Rathinaraj et al.<fnref n="28" /> for additional details).

## REFERENCES

<footnote n="1">¹H. M. Laun, "Prediction of elastic strains of polymer melts in shear and elongation," J. Rheol. 30, 459-501 (1986).</footnote>

<footnote n="2">²R. B. Bird, R. C. Armstrong, and O. Hassager, *Dynamics of Polymeric Liquids*, Vol. 1: *Fluid Mechanics* (Wiley, 1987).</footnote>

<footnote n="3">³M. M. Denn, *Polymer Melt Processing: Foundations in Fluid Mechanics and Heat Transfer* (Cambridge University Press, 2008).</footnote>

<footnote n="4">⁴R. I. Tanner, *Engineering Rheology* (OUP Oxford, 2000).</footnote>

<footnote n="5">⁵Dedicated to Prof. Roger I. Tanner FRS on his 90th birthday.</footnote>

<footnote n="6">⁶C. A. García-Franco, "A note on the elasticity of polymer melts described by primary normal stress difference ($N_1$)," in *Macromolecular Symposia* (Wiley Online Library, 2013), Vol. 325, pp. 184-193.</footnote>

<footnote n="7">⁷V. Sharma and G. H. McKinley, "An intriguing empirical rule for computing the first normal stress difference from steady shear viscosity data for concentrated polymer solutions and melts," *Rheol. Acta* 51, 487-495 (2012).</footnote>

<footnote n="8">⁸C. Friedrich and L. Heymann, "Primary normal-stress coefficient prediction at high shear rates," *Rheol. Acta* 27, 567-574 (1988).</footnote>

<footnote n="9">⁹D. Tang, F. H. Marchesini, L. Cardon, and D. R. D'hooge, "Evaluating the exit pressure method for measurements of normal stress difference at high shear rates," *J. Rheol.* 64, 739-750 (2020).</footnote>

<footnote n="10">¹⁰J. M. Dealy and K. F. Wissbrun, *Melt Rheology and Its Role in Plastics Processing: Theory and Applications* (Springer Science & Business Media, 2012).</footnote>

<footnote n="11">¹¹R. Steller, "Determination of the first normal stress difference from viscometric data for shear flows of polymer liquids," *Rheol. Acta* 55, 649-656 (2016).</footnote>

<footnote n="12">¹²R. I. Tanner, "A theory of die-swell," *J. Polym. Sci. Part A-2* 8, 2067-2078 (1970).</footnote>

<footnote n="13">¹³R. I. Tanner and M. Keentok, "Shear fracture in cone-plate rheometry," *J. Rheol.* 27, 47-57 (1983).</footnote>

<footnote n="14">¹⁴R. I. Tanner, "Some methods for estimating the normal stress functions in viscometric flows," *Trans. Soc. Rheol.* 14, 483-507 (1970).</footnote>

<footnote n="15">¹⁵É. Couturier, F. Boyer, O. Pouliquen, and É. Guazzelli, "Suspensions in a tilted trough: Second normal stress difference," *J. Fluid Mech.* 686, 26-39 (2011).</footnote>

<footnote n="16">¹⁶S.-C. Dai, E. Bertevas, F. Qi, and R. I. Tanner, "Viscometric functions for non-colloidal sphere suspensions with Newtonian matrices," *J. Rheol.* 57, 493-510 (2013).</footnote>

<footnote n="17">¹⁷R. I. Tanner, "Rheology of noncolloidal suspensions with non-Newtonian matrices," *J. Rheol.* 63, 705-717 (2019).</footnote>

<footnote n="18">¹⁸S. E. Mall-Gleissle, W. Gleissle, G. H. McKinley, and H. Buggisch, "The normal stress behaviour of suspensions with viscoelastic matrix fluids," *Rheol. Acta* 41, 61-76 (2002).</footnote>

<footnote n="19">¹⁹O. Maklad and R. J. Poole, "A review of the second normal-stress difference; its importance in various flows, measurement techniques, results for various complex fluids and theoretical predictions," *J. Non-Newtonian Fluid Mech.* 292, 104522 (2021).</footnote>

<footnote n="20">²⁰R. G. Larson, "Constitutive relationships for polymeric materials with power-law distributions of relaxation times," *Rheol. Acta* 24, 327-334 (1985).</footnote>

<footnote n="21">²¹D. G. Baird, "First normal stress difference measurements for polymer melts at high shear rates in a slit-die using hole and exit pressure data," *J. Non-Newtonian Fluid Mech.* 148, 13-23 (2008).</footnote>

<footnote n="22">²²M. Keentok, A. G. Georgescu, A. A. Sherwood, and R. I. Tanner, "The measurement of the second normal stress difference for some polymer solutions," *J. Non-Newtonian Fluid Mech.* 6, 303-324 (1980).</footnote>

<footnote n="23">²³R. V. More, R. Patterson, E. Pashkovski, and G. H. McKinley, "Rod-climbing rheometry revisited," *Soft Matter* 19, 4073-4087 (2023).</footnote>

<footnote n="24">²⁴W. P. Cox and E. H. Merz, "Correlation of dynamic and steady flow viscosities," *J. Polym. Sci.* 28, 619-622 (1958).</footnote>

<footnote n="25">²⁵Y. H. Wen, H. C. Lin, C. H. Li, and C. C. Hua, "An experimental appraisal of the Cox-Merz rule and Laun's rule based on bidisperse entangled polystyrene solutions," *Polymer* 45, 8551-8559 (2004).</footnote>

<footnote n="26">²⁶P. F. Teixeira, L. Hilliou, J. A. Covas, and J. M. Maia, "Assessing the practical utility of the hole-pressure method for the in-line rheological characterization of polymer melts," *Rheol. Acta* 52, 661-672 (2013).</footnote>

<footnote n="27">²⁷R. Kádár, I. F. C. Naue, and M. Wilhelm, "First normal stress difference and in-situ spectral dynamics in a high sensitivity extrusion die for capillary rheometry via the 'hole effect'," *Polymer* 104, 193-203 (2016).</footnote>

<footnote n="28">²⁸A. Jaishankar and G. H. McKinley, "A fractional K-BKZ constitutive formulation for describing the nonlinear rheology of multiscale complex fluids," *J. Rheol.* 58, 1751-1788 (2014).</footnote>

<footnote n="29">²⁹J. D. J. Rathinaraj, B. Keshavarz, and G. H. McKinley, "Why the Cox-Merz rule and Gleissle mirror relation work: A quantitative analysis using the Wagner integral framework with a fractional Maxwell kernel," *Phys. Fluids* 34, 033106 (2022).</footnote>

<footnote n="30">³⁰R. I. Tanner, "From A to (BK) Z in constitutive relations," *J. Rheol.* 32, 673-702 (1988).</footnote>

<footnote n="31">³¹R. G. Larson, *Constitutive Equations for Polymer Melts and Solutions* (Butterworth-Heinemann, 2013).</footnote>

<footnote n="32">³²J. M. Dealy and R. G. Larson, *Structure and Rheology of Molten Polymers: From Polymerization to Processability via Rheology* (Hanser Publishers, 2006).</footnote>
<footnote n="33">³³P. R. Soskey and H. H. Winter, "Large step shear strain experiments with parallel-disk rotational rheometers," *J. Rheol.* **28**, 625–645 (1984).</footnote>

<footnote n="34">³⁴A. C. Papanastasiou, L. E. Scriven, and C. W. Macosko, "An integral constitutive equation for mixed flows: Viscoelastic characterization," *J. Rheol.* **27**, 387–410 (1983).</footnote>

<footnote n="35">³⁵L. I. Palade, V. Verney, and P. Attané, "A modified fractional model to describe the entire viscoelastic behavior of polybutadienes from flow to glassy regime," *Rheol. Acta* **35**, 265–273 (1996).</footnote>

<footnote n="36">³⁶A. Lion, "On the thermodynamics of fractional damping elements," *Continuum Mech. Thermodyn.* **9**, 83–96 (1997).</footnote>

<footnote n="37">³⁷J. D. J. Rathinaraj, G. H. McKinley, and B. Keshavarz, "Incorporating rheological nonlinearity into fractional calculus descriptions of fractal matter and multi-scale complex fluids," *Fractal Fractional* **5**, 174 (2021).</footnote>

<footnote n="38">³⁸L. I. Palade, P. Attané, R. R. Huilgol, and B. Mena, "Anomalous stability behavior of a properly invariant constitutive equation which generalises fractional derivative models," *Int. J. Eng. Sci.* **37**, 315–329 (1999).</footnote>

<footnote n="39">³⁹A. Heibig and L. I. Palade, "On the rest state stability of an objective fractional derivative viscoelastic fluid model," *J. Math. Phys.* **49**, 043101 (2008).</footnote>

<footnote n="40">⁴⁰A. Heibig and L. I. Palade, "Well posedness of a linearized fractional derivative fluid model," *J. Math. Anal. Appl.* **380**, 188–203 (2011).</footnote>

<footnote n="41">⁴¹M. H. Wagner and J. Meissner, "Network disentanglement and time-dependent flow behaviour of polymer melts," *Die Makromol. Chem.* **181**, 1533–1550 (1980).</footnote>

<footnote n="42">⁴²N. W. Tschoegl, *The Phenomenological Theory of Linear Viscoelastic Behavior: An Introduction* (Springer-Verlag, Berlin, Heidelberg, 1989).</footnote>

<footnote n="43">⁴³I. Podlubny, *Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications* (Elsevier, 1998).</footnote>

<footnote n="44">⁴⁴V. H. Rolón-Garrido and M. H. Wagner, "The damping function in rheology," *Rheol. Acta* **48**, 245–284 (2009).</footnote>

<footnote n="45">⁴⁵R. L. Bagley and P. J. Torvik, "A theoretical basis for the application of fractional calculus to viscoelasticity," *J. Rheol.* **27**, 201–210 (1983).</footnote>

<footnote n="46">⁴⁶A. W. Wharmby and R. L. Bagley, "Generalization of a theoretical basis for the application of fractional calculus to viscoelasticity," *J. Rheol.* **57**, 1429–1440 (2013).</footnote>

<footnote n="47">⁴⁷A. Jaishankar and G. H. McKinley, "Power-law rheology in the bulk and at the interface: Quasi-properties and fractional constitutive equations," *Proc. R. Soc. A* **469**, 20120284 (2013).</footnote>

<footnote n="48">⁴⁸C. Friedrich and H. Braun, "Generalized Cole-Cole behavior and its rheological relevance," *Rheol. Acta* **31**, 309–322 (1992).</footnote>