An elastoviscoplastic constitutive equation is proposed to describe both the elastic and rate-dependent plastic deformation behavior of Carbopol^{®} dispersions, commonly used to study yield-stress fluids. The model, a variant of the nonlinear Maxwell model with stress-dependent relaxation time, eliminates the need for a separate Herschel–Bulkley yield stress. The stress dependence of the viscosity was determined experimentally by evaluating the steady-state flow stress at a constant applied shear rate and by measuring the steady-state creep rate at constant applied shear stress. Experimentally, the viscosity’s stress-dependence was confirmed to follow the Ree–Eyring model. Furthermore, it is shown that the Carbopol^{®} dispersions used here obey time-stress superposition, indicating that all relaxation times experience the same stress dependence. This was demonstrated by building a compliance mastercurve using horizontal shifting on a logarithmic time axis of creep curves measured at different stress levels and by constructing mastercurves of the storage- and loss-modulus curves determined independently by orthogonal superposition measurements at different applied constant shear stresses. Overall, the key feature of the proposed constitutive equation is its incorporation of a nonlinear stress-activated change in relaxation time, which enables a smooth transition from elastic to viscous behavior during start-up flow experiments. This approach bypasses the need for a distinct Herschel–Bulkley yield stress as a separate material characteristic. Additionally, the model successfully replicates the observed steady-state flow stress in transient-flow scenarios and the steady-state flow rate in creep experiments, underlining its effectiveness in capturing the material’s dynamic response. Finally, the one-dimensional description is readily extended to a full three-dimensional finite-strain elastoviscoplastic constitutive equation.

## I. INTRODUCTION

Simple or ideal yield-stress fluids is a designation often used for nonthixotropic viscoplastic soft materials that exhibit transition from solidlike behavior at low stress levels to flow above a critical stress [1–3]. Many real-world examples fall into this category, ranging from consumer products to food compounds [2], and various levels of description of increasing complexity are available to characterize the rheological behavior of these interesting materials.

The simplest model for describing yield stress fluid behavior is the two-parameter Bingham model. This model assumes rigid behavior (no deformation) below, and Newtonian fluid behavior by viscous flow above a constant yield stress, enabling analytical solutions for many relevant flow problems [4,5]. More complexity is allowed for by the three parameter Herschel–Bulkley [6] or the two parameter Casson [7] models, providing non-Newtonian viscous flow behavior above a limiting value of stress, named the “yield stress.” These models are designed to accurately depict flow stress in relation to the strain rate under *steady-state conditions*, rather than transient behavior. It is important to highlight the contrast in yield stress definitions between these models and those found in solid-state (polymer) mechanics textbooks. In solid-state mechanics, the yield stress is usually defined as the peak stress on a stress-strain curve, marking the shift from elastic to plastic deformation in a *transient experiment*. The yield stress of solid-state mechanics, defined as maximum stress in a stress-strain experiment, is typically strain-rate dependent, and, consequently, all flow-stress values in a Herschel–Bulkley plot would be regarded as strain-rate-dependent yield stresses in the solid-state mechanics community. To clarify this distinction and to avoid confusion, we will denote the limiting yield stress in a Herschel–Bulkley plot as the “Herschel–Bulkley yield stress.”

The assumption of the existence of a Herschel–Bulkley yield stress, below which there is only elastic deformation, has met with criticism [8] due to the emergence of a limiting zero-shear plateau value of the viscosity for typical yield-stress fluid systems such as Carbopol^{®} microgel dispersions of various volume fractions [9]. Experiments by other groups, however, showed that the plateau zero-shear viscosity in Carbopol^{®} microgels increases in time due to stress-induced aging [10,11], a phenomenon that is well known from glassy polymer systems [12,13].

The linear and nonlinear viscoelastic response of yield-stress fluids has received less attention, and if addressed, most research was directed toward creep below the yield stress in relation to the debate about the existence of a true yield stress. Caton and Baravian [14], Lidon *et al.* [15], and others observed Andrade-like creep in Carbopol^{®} microgel dispersions below the yield stress, where the creep rate depends on time through a power-law, the physical origin of which is, however, not fully understood and may arise from the interplay between thermal activation and elastic stress redistribution. Aime *et al.* [16] showed that the *linear viscoelastic behavior* of a fractal colloidal gel consisting of charged silica particles was accurately described by a fractional Maxwell model. Saramito [17] derived an *elastoviscoplastic* model, extending both the Bingham viscoplastic and the Oldroyd viscoelastic models, and later replaced the Bingham model by the Herschel–Bulkley viscoplastic model [18]. Also, Dimitriou *et al.* [19] focused on the *yielding behavior* of Carbopol^{®} dispersions, regarding this material as an elastoplastic material combined with a kinematic-hardening model, which was later extended to a full three-dimensional constitutive description of the yielding process more from the perspective of plasticity [20]. Fielding [21] presented a detailed *mesoscopic* elastoviscoplastic equation based on the soft glassy rheology model [22], which was shown to capture many aspects of the linear and nonlinear viscoelastic behavior of soft glassy materials, with a focus on rheological aging, and also pertains to a class of elastoplastic models [23].

The aim of this study is to create a macroscopic phenomenological constitutive model that captures both linear and nonlinear viscoelastic behaviors, as well as the transition to shear-rate dependent yielding in Carbopol^{®}, a prototypical simple yield-stress fluid. This model sidesteps the need for an explicit Herschel–Bulkley yield stress. Instead, it employs concepts from the plastic flow of amorphous materials to explain stress-activated flow in dense suspensions, though it does not cover the entire spectrum of elastoplastic material behavior. The conjecture is that the specific Carbopol^{®} dispersion used in this study can be regarded as an example of an ideal elastoviscoplastic fluid, in which stress-activated events relax the elastic stress. Our approach involves examining stress activation by scaling creep experiments and conducting independent assessments of how stress influences the relaxation during flow, utilizing orthogonal superposition rheometry. The core of our proposed constitutive equation is a simple nonlinear viscoelastic model with a stress-dependent relaxation time. This key feature effectively accounts for stress impacts on plastic flow. We will demonstrate how this model can be expanded into a comprehensive three-dimensional nonlinear viscoelastic constitutive equation.

## II. THEORY

Here, $ \sigma d$ is the deviatoric part of the Cauchy stress tensor, $ D p$ is the plastic strain-rate tensor, and $ \gamma \u02d9 eq$ is the equivalent shear rate, which depends on the second invariant of the plastic strain-rate tensor $ D p$ as $ \gamma \u02d9 eq= 2 D p : D p$ and which reduces to the norm of the shear rate $ \gamma \u02d9$ for the case of simple shear flow. It should be noted that $ D p$ is defined as the symmetric part of the plastic velocity gradient and $ D p= D p d$ as plastic flow is assumed to proceed at constant volume [27]. The Bird–Carreau equation predicts non-Newtonian flow behavior with a “zero-shear viscosity” $ \eta 0$ when the equivalent shear rate is smaller than a characteristic shear rate $ \gamma \u02d9 0$.

*stress*-dependent power-law viscosity $ \eta PL( \sigma eq)$ as

It should be noted that Eq. (11) is a single viscosity expression, in which both structural relaxation processes experience the same stress and strain rate. As such, it is fundamentally different from the original Ree–Eyring model, Eq. (10), where each process experiences a different part of the total stress. Nevertheless, as the total flow stress as a function of strain rate for Eq. (11) is almost indistinguishable from the original Ree–Eyring model, Eq. (10), we will also refer to Eq. (11) as the “Ree–Eyring model.”

The Eyring and Ree–Eyring models are based on the idea that shear flow is a stress-biased thermally activated process [29]. The colloidal suspensions studied here are expected to be athermal, yet they are also shear activated. The viscosity of soft particle suspensions has been observed to scale exponentially with an excess entropy, which characterizes the microstructure thermodynamically [40]. Hence, for a system that is shear activated, but with negligible thermal forces, the functional form of the Eyring model, usually limited to thermal systems, is a rational and thermodynamics-based choice, also for these predominantly athermal systems.

A comparison of the stress-dependent shift factor for the power-law model, the Eyring and the Ree–Eyring equations are shown in Fig. 1. Here, it can be seen that at equal $ \sigma 0$, depending on the power-law index, the power-law fluid departs from linear behavior $( a \sigma =1)$ at a steeper slope compared to the Eyring equation. However, especially at $n\u22480.1$, the shape of the two curves is quite similar, which means they can be made to approximately overlap by selecting the appropriate $ \sigma 0$ value.

Figure 1 also shows that using the Ree–Eyring function it is possible to introduce additional plateaus in the viscosity as a function of stress. There is another viscosity model that allows for this feature, which is the SMD-model from Souza Mendes and Dutra [41], which was used in a strain-dependent elastoviscoplastic model [42]. So, rather than invoking a separate material function (the Herschel–Bulkley yield stress), simple elastoplastic models use highly nonlinear functional forms of stress-dependent shift factors which reflect the stress-activated nature of plastic flow, which acts to unload the elastic stress.

In contrast to the Herschel–Buckley model, this nonlinear viscoelastic approach allows for simultaneous elastic and stress-activated plastic deformation. During start-up flow with a constant applied shear rate, the accumulation of elastic strain reduces as the plastic shear rate increases. In the solid-state polymer plasticity literature, the *yield stress* in case of ideal plasticity (no strain hardening or strain softening) is then defined as the stress at which the elastic shear rate becomes zero and the applied shear rate equals the plastic shear rate [44]. During plastic flow, the stress is then per definition equal to the yield stress, which might be pressure, temperature, and strain-rate dependent, see, for example, Bauwens-Crowet *et al.* [45]. Interestingly, this elastoviscoplastic definition of the yield stress is different from how it is defined in the rheological literature, where, especially in relation to the Herschel–Bulkley model, the term “yield stress” typically is reserved for the limiting value of the flow stress at low shear rates, below which it is assumed that there is no plastic deformation [1]. To avoid confusion, we have chosen to use the term “Herschel–Bulkley yield stress” to refer to this rheologically defined yield stress, obtained by extrapolation of the flow stress to the zero shear rate.

The evolution equations, Eqs. (12a)–(12d), depict a standard nonlinear viscoelastic formulation of the rate-dependent yielding of a solid material without the need for a yield criterion. At low stress levels, $ \sigma eq\u226a \sigma 0$, as the stress-shift factor $ a \sigma =1$, the viscosity equals $ \eta 0$, which is typically high, resulting in predominantly elastic behavior. At higher stress levels, the stress-dependent viscosity rapidly decreases nonlinearly toward lower values, leading to an increasing plastic shear rate $ \gamma \u02d9 p$ and, finally, to a fully developed plastic flow, defined by an elastic shear rate that becomes zero during the experimental time-frame. The yielding condition of zero elastic shear rate can then be used to experimentally determine the stress-dependent viscosity during plastic flow. It should be noted, that in the limit of low stress this model reduces to a linear viscoelastic Maxwell model, implying there is always plastic flow as long as the stress is nonzero.

*stress-induced mobility*, the use of which is sometimes also referred to as “stress-clock materials” [46].

In what follows, this model will be applied to a dense microgel suspension (a particular grade of Carbopol^{®} microgels) as an example of a yield-stress fluid, and some implications of the proposed model will be discussed, including a full three-dimensional elastoviscoplastic constitutive formulation, valid for both finite (large) elastic and plastic deformations.

## III. MATERIALS AND METHODS

A batch of 0.5 wt. % Carbopol^{®}940 in water was prepared following a protocol adapted from Varges *et al.* [47] Briefly, dry powder (4 g) of Carbopol^{®}940 (Acros Organic) was sieved using a stainless steel coarse sieve (mesh size: 1 mm^{2}, Ikea), and added to 790 g of MilliQ water while continuously stirred at 800 rpm. The stirring was conducted in a Heidolph RZR2100 mixer motor with a Rushton six-bladed turbine (diameter: 53 mm, power number: $\u223c0.12$, and pumping number: $0.72$). This mixture was stirred for 20 min at 800 rpm after which it was let to rest for 30 min. 6 g of triethanolamine (TEA) from Sigma-Aldric ( $99%$ purity) was added to the mixture while gently stirring at 240 rpm. Once TEA was added, the mixture was kept stirring at 240 rpm for three days to ensure homogeneity of the batch. After one day of mixing, two $ 45 \xb0$ pitched blade turbine (diameter: 40 mm, power number: $\u223c0.53$, and pumping number: $0.79$) were added to the setup. Evaporation was limited by wrapping the vessel with Parafilm and placing two gaskets on the impeller shaft on opposite sides of the bottle cap. The applied protocol lead to samples that were not thixotropic on experimental accessible time scales [48].

Rheological measurements were performed with a modular compact rheometer (MCR) 702 from Anton Paar. For the orthogonal superposition (OSP) measurements, a prototype configuration for the MCR702 (Anton Paar, Austria) using a linear motor mounted in the lower flange was employed to perform small amplitude oscillation orthogonal to a steady flow generated by the standard rotational motor. The prototype works with a double-wall Couette geometry designed for orthogonal superposition measurements (length: 40 mm, average radius: 12.8 mm, average radii ratio: 1.03 , average gap: 0.36 mm). The sample gap is open at the bottom and connected to an inner reservoir permitting the flow of the material, thus reducing pumping effects [49]. The Couette geometry was sandblasted to reduce wall slip, with asperity heights between 4 and $7 \mu m$. The same geometry was employed to perform creep and start-up measurements. In this case only the rotational motor performed the experiments. The limit of the instrument and the data refinement for OSP measurements is pointed out in Appendix B.

Slip effects at the geometry surfaces are reported in the literature for viscoplastic materials [50], and they are expected in our sample. These are characterized by a deviation in the flow curve at low shear rate, with an apparent material response more compliant than the expected one. Rough geometries should reduce this deviation and suppress wall-slip, with the exception of small gaps where confinement effects are significant. Thus, the effect of wall-slip was checked by performing measurements at different gap distances in a serrated plate-plate configuration for the lowest stress applied during experiments ( $\sigma =$ 20 Pa), Fig. 2, and by comparing those with the sandblasted Couette used in this paper and a smooth version of the same Couette. The serrated plate-plate configuration (diameter: 40 mm) used for the verification experiments has a custom-made mesh with cubic pillars ( $250\xd7250\xd7250 \mu m$) spaced 190 μm in both planar Cartesian directions. The focus on a single stress value during our analysis force us to use a single-point correction of the response obtained with plate-plate geometry, which is automatically implemented by the instrument, for all measurements, by using an effective radius of $2/3$ of the maximum radius for all calculations. For rough geometries, the curves show geometry independent response until a considerable small gap was considered (H=400 μm), and even in the latter case, the material response does not show the same intensity of slip as in the smooth Couette case. The geometry independent response and the use of the lower extreme of the experimental stress range validates all data presented in this study against the occurrence of wall-slip.

Since the material is also sensitive to evaporation, which limits the quality of the data at long measurement times, great care was taken into account to minimize it using solvent traps or immersion. The reliability of each measurement was confirmed through a set of two fast constant shear rate experiments ( $ \gamma \u02d9=1 s \u2212 1$ for 30 s), before and after the measurement, to check that no changes in the experimental curves were present. In addition, a small waiting time of 30 s was introduce between the first control experiment and the proper measurement. For start-up flows, the material was covered with low viscous paraffin oil (technical grade, VWR Chemicals) and experiments were conducted in one session without noticeable evaporation. An additional step was performed for creep measurements in order to avoid excessive impact of inerto-elastic effects. The sample was conditioned with a stress ramp to the desired stress from 0 to 0.1 s, then the data until 10 s were not considered to have a conservative cut-off of possible inerto-elastic effects.

## IV. RESULTS AND DISCUSSION

### A. Start-up flow

Transient measurements performed at shear rates of $ \gamma \u02d9= 5 \xd7 10 \u2212 4 s$ to 1 s^{−1} are shown in Fig. 3. From this figure it follows that, for the strain rates applied in this study, the specific Carbopol^{®} dispersions behave to a good approximation as simple yield-stress fluids, depicting (visco)elastic deformation followed by plastic deformation without overshoot (no strain softening or strain hardening).

According to the elastoviscoplastic model discussed in the Theory section, at the yield point, the elastic shear rate becomes zero and the plastic shear rate equals the total applied shear rate. Therefore, after reaching a steady-state flow stress $\sigma $ during a start-up experiment at constant applied shear rate $ \gamma \u02d9$, a steady-state viscosity $\eta $ can be calculated as $\eta =\sigma / \gamma \u02d9$ [44]. In addition, flow-stress values were determined from steady-state flow derived from the plateau-creep rate during creep experiments (see Sec. IV B). These flow stresses, determined from steady-state flow during transient and creep experiments, can be represented in three equivalent ways. First, by depicting the steady-state flow stress as a function of the applied shear rate (Fig. 4); second, by plotting the steady-state viscosity as a function of applied shear rate [Fig. 5(a)]; and third, by plotting the same viscosity as a function of the measured steady-state flow stress [Fig. 5(b)]. From Fig. 4, it is immediately clear that the Eyring model will not be able to describe these data as it is well known that the Eyring model appears as a straight line on a semi-logarithmic plot of the flow stress versus strain rate [44].

From Fig. 5(a), it can be seen that the viscosities, to a good approximation, appear to depict a power-law dependence on the shear-rate over the experimental range of shear-rates used. Therefore, using this representation would only allow an accurate determination of the power-law index $n$. Better fitting opportunities are provided by both the flow stress as a function of the shear rate (Fig. 4) and the viscosity as a function of stress [Fig. 5(b)]. In what follows, nonlinear fitting was performed using Mathematica® on the flow stress as a function of shear rate (Fig. 4), using Eqs. (1), (3), and (10), for the Ree–Eyring, Herschel–Bulkley, and power-law model, respectively. The material parameters determined were then used to calculate both the viscosity as a function of shear rate and shear stress, respectively [solid lines in Figs. 5(a) and 5(b)], by applying Eqs. (7) and (11), using the relations $\sigma ( \gamma \u02d9 e q)=\eta ( \gamma \u02d9 e q ) \gamma \u02d9 e q$ and $ \sigma 0= \eta 0 \gamma \u02d9 0$. Figure 4 shows that the Herschel–Bulkley model gives a good description of the flow-stress as a function of shear rate, but will not be considered here, as the viscosity of the Herschel–Bulkley model will approach infinity in the limit of zero shear rate. Comparing the power-law and Ree–Eyring model, it can be seen in Fig. 4 that the Ree–Eyring function provides the best fit to the experimental data and will thus be used in what follows. The Ree–Eyring parameters determined from the fit of Eq. (10) in Fig. 4 are $ \sigma 01=( 2.7 \xb1 0.2) Pa$, $ \sigma 02=10.6\xb10.7 Pa$, $log\u2061( \gamma \u02d9 01/ s \u2212 1)=\u221214.7\xb10.7$, and $log\u2061( \gamma \u02d9 02/ s \u2212 1)=\u22120.34\xb10.08$, from which it can be calculated, using Eq. (11), that $ \eta 01=1.59\xd7 10 15$ and $ \eta 02=3.0\xd7 10 4 Pa s$.

It should be noted that the available viscosity data as depicted in Figs. 5(a) and 5(b), do not allow for an accurate determination of the zero-shear viscosity, as neither the viscosity as a function of stress nor as a function of shear rate, clearly suggest an emerging plateau at low shear stress or shear rate. It should also be noted that the change in slope as depicted in Fig. 5(b), has been observed by others, for example, Erwin *et al.* [51]. Finally, the initial elastic shear modulus $G$ was chosen to fit the initial slope of the transient measurements at the highest shear rate and was found to be approximately $G=300 Pa$.

Having determined the necessary material parameters, predictions for shear-rate dependent stress-strain experiments according to the elastoviscoplastic model, Eq. (12), can be obtained by numerical integration of the evolution Eqs. (12a)–(12d) using Mathematica $ \xae$. This has been done for various applied shear rates, ranging from $ \gamma \u02d9=5\xd7 10 \u2212 4$ to 1 s $ \u2212 1$, and these are shown in Fig. 3. As can be seen, the model describes a sharp but continuous transition from elastic to plastic deformation, without the need to invoke an explicit Herschel–Bulkley yield stress. There is a good agreement between the experimental and predicted values of the steady-state flow-stress levels. This is not surprising, but simply reflects the good fit of the steady-state flow stresses with the Ree–Eyring equation depicted in Fig. 4. For the single relaxation time model [Eqs. (12)], the transition from elastic to plastic flow behavior shown in Fig. 3 is rather abrupt compared to the experimental data. However, using a spectrum of relaxation times, all endowed with the same stress dependence, is expected to result in a more smooth strain-rate dependent transition from elastic to plastic deformation and will also correctly describe the linear and nonlinear viscoelastic behavior (see Sec. IV D) [30].

### B. Creep experiments

When using a spectrum of relaxation times, the nonlinear behavior of an elastoviscoplastic material can then be described by assuming that all relaxation times obey the same stress dependence [30,52,53]. The similar stress dependence of the relaxation times can be critically investigated by performing constant-stress experiments (creep testing) at different stress levels and look for time-stress superposition in much the same way as the well-known time-temperature superposition (where it is assumed that all relaxation times have the same temperature dependence). If in the present case a time-stress superposition holds, it should be possible to build a creep mastercurve using creep tests performed at different stress levels (see Fig. 6) by horizontal shifting along a logarithmic time axis. This is shown for the Carbopol^{®} formulation used in this study in Fig. 7. In this figure, the reference stress for the master curve was chosen as 20 Pa. As can be seen, creep curves measured at higher stresses than 20 Pa, shift to longer times (to the right), while creep curves measured at lower stresses than 20 Pa would shift to shorter times (to the left). All curves overlap relatively smoothly to form a creep compliance mastercurve that spans close to 14 (!) decades. Shifting was performed according to the closed form shifting (CFS) algorithm from Gergesova *et al.* [54]. This algorithm avoids the ambiguity generated by other procedures, e.g., manual shifting or preliminary curve fitting of the experimental data, through a minimization of the area between two data sets at different applied stress. The authors of the paper showed the reliability of their algorithm estimating an error of few percent when applied to experimental data with up to $15 %$ noise.

It should be noted that the creep compliance mastercurve in Fig. 7 is a virtual curve that does not take into account potential molecular relaxation and inertia effects at very short times and possible effects of physical aging at long times [55]. In reality, especially physical aging during a true *experimental* creep test over long time scales, possibly accelerated by the application of small stress, could make an experimental creep compliance curve deviate considerably from the virtual creep compliance mastercurve depicted in Fig. 7 [10–12]. It should also be noted that the true virtual *linear* creep compliance mastercurve can now be obtained by shifting the mastercurve at 20 Pa in Fig. 7 to the right by an amount of $1/ a \sigma RE(20 Pa)$.

The shift factors necessary to build the creep mastercurve at 20 Pa are shown in Fig. 8, together with a line, calculated using the Ree–Eyring model, Eq. (14), using parameters determined from the fit of the measured viscosities from the transient and creep experiments to the Ree–Eyring model (see Fig. 4). As can be seen, also for the shift factors, the Ree–Eyring model gives a good description.

If steady-state flow is reached during a creep test, the elastic shear rate becomes zero and the plastic shear rate equals the total (measured) shear rate, also called the “plateau-creep rate.” The plateau-creep rate is, therefore, experimentally determined as the moment that the creep rate becomes constant, which can be conveniently determined experimentally by plotting the creep rate as a function of creep strain. In such a so-called “Sherby–Dorn” plot [56], the constant plateau-creep rate can easily be identified as a minimum (or a horizontal line).

^{-3}s

^{−1}, giving rise to a viscosity of

### C. Orthogonal superposition rheometry

Another type of constant stress measurements that can be used to directly interrogate the stress activation of the different relaxation times during plastic deformation are orthogonal superposition experiments. Orthogonal superposition (OSP) rheometry consists of the simultaneous application of a constant shear rate (resulting in a constant applied stress at the steady state) and a small amplitude oscillation shear in an orthogonal direction [57]. The frequency sweep retains the linear viscoelastic interpretation [58] because, in a first approximation (assuming an isotropic structure), the two flows are decoupled, which permits to study the relaxation spectra during flow [59–62].

The experimental results, displayed in Fig. 10, show how varying stress levels affect the orthogonal storage ( $ G \u2032\u22a5$) and the orthogonal loss moduli (G” $ \u22a5$). It should be noted that in the absence of stress, the Carbopol^{®}940 sample behaves elastically, meaning that $ G \u2032=G\u2248300 Pa$ and $ G \u2033\u22480 Pa$. These graphs make it evident that increasing stress significantly alters the curves from their linear viscoelastic response at zero stress. Notably, a crossover point between the moduli emerges at stresses of 117 Pa. This crossover indicates the material’s relaxation frequency under the applied stress. Furthermore, as this crossover shifts to higher frequencies with increasing stress, it confirms the acceleration of material relaxation and supports the notion that flow in these materials is stress-activated. Assuming isotropic behavior and neglecting the minute stresses due to the linear viscoelastic measurement in the orthogonal direction, the equivalent stress [see Eq. (4)] during an orthogonal superposition measurement is equal to the applied shear stress. In addition, the viscosity at that stress is simply the measured viscosity during the OSP experiment. These viscosities can be compared to the experimental viscosities presented in Figs. 5(a), 5(b), and 12. Selecting a reference stress and using the measured viscosity at that stress, the shift factor with respect to other stresses can now be calculated according to Eq. (15). In this manner, the shift factors obtained from creep measurements were used to shift the moduli, as depicted in Fig. 11, retrieving a single mastercurve for the orthogonal moduli.

Similar shifting procedures are found in the literature for concentrated colloidal dispersions [63,64], showing a correlation between the relaxation frequency and the applied shear rate. Dispersions with an effective volume fraction above the glass transition appear only to shift in the time-frequency domain. The absence of shifting in the moduli intensity is explained by the constrained structure that limits the available free-volume. These considerations can be used to get insight into the structure of Carbopol^{®}940, associating the microscopic dynamic to the one of glassy materials. Furthermore, Sung *et al.* [65] presented the shifting of dynamic data for colloidal gels with an open gel structure (rather than the dense suspension microstructure in the Carbopol^{®}940), where they compare the obtained shift factors with the viscosity flow-curve. Similarly to what is done in the present study, they show a strong agreement between the two curves and a smooth shifting of the OSP moduli. The similarity of behaviors in materials different from a microstructural point of view shows the underlying importance of stress activated plastic events in the plastic flow of complex fluids.

The estimated viscosities from the orthogonal superposition measurements are included in Fig. 12. The agreement of the extrapolated viscosities with the measured viscosities and the creation of a single master-curve give additional credit to the assumption of stress-induced mobility as the main driving force in plastic flow of Carbopol. The solid line in Fig. 12 represents the best fit of Ree–Eyring model, Eq. (11), to the experimental viscosity values obtained from the steady-state flow stresses during start-up flow and from the steady-state creep rates (see Fig. 4).

It should be noted that in the original Ree–Eyring formulation, Eq. (10), strain is assumed to be equal for all molecular processes, whereas the stresses due to the different processes are assumed to be additive. Therefore, in case where the original formulation of the Ree–Eyring function applies, time-stress superposition would be expected to fail in the stress region where the plastic strain rate is determined by both molecular processes [34]. This was not observed in this research, possibly while the stress range where this failure would occur is small compared to the experimental stresses used in this study. The successful construction of the creep compliance and orthogonal moduli mastercurves, whereby the resulting shift factors follow a physically “realistic” function (in this case the Ree–Eyring function), are a strong indication that time-stress superposition applies [66].

The Herschel–Bulkley description with a single rate-independent Herschel–Bulkley yield stress will not be able to capture the experimental plastic creep deformation behavior below that Herschel–Bulkley yield stress. In contrast, the present nonlinear viscoelastic constitutive equation, Eq. (12), accurately describes both the elastic response, the rate-dependent plastic deformation, and the nonlinear steady-state creep rate, all without the need for a Herschel–Bulkley yield stress. In the current model, the Carbopol^{®} dispersion are regarded as a nonlinear viscoelastic *liquid*. Even in the limit of very low stress, the material will still behave like a viscoelastic liquid, albeit with an extremely high viscosity, resulting in virtually pure elastic deformation at small strains (low stresses) on realistic (accessible) time scales.

### D. Multimode model

*linear*shear-relaxation modulus $G(t)$, we start from the creep-compliance mastercurve at 20 Pa. From this curve, the linear compliance curve was obtained by shifting the creep-compliance curve at 20 to 0 Pa, as shown in Fig. 13. Subsequently, the linear compliance curve was fitted to the fractional Maxwell-liquid model [67,68], consisting of a dashpot and a “spring-pot” in series, for which the creep compliance is given by

A full nonlinear viscoelastic constitutive description of the Carbopol^{®} sample used in this study is now obtained by combining Eq. (16) with the experimentally determined relaxation-time spectrum Eq. (18). For a start-up flow with a constant applied shear rate $ \gamma \u02d9 0$, Eqs. (16) and (18) were solved for the shear stress by iteration for each time step using Mathematica®. The comparison of a single-mode and a multimode response is given in Fig. 15. Here, it can be seen that the multimode description results in a more gradual transition of elastic to plastic deformation and accurately follows the experimental curves up to applied shear rates in the order of $ \gamma \u02d9=0.1 s \u2212 1$. For higher shear rates, the experimental curves are more steep. The multimode model still predicts the correct steady-state flow stress but depicts a too slow transition from elastic to plastic flow. This could be due to the transition to thixotropic flow behavior of the Carbopol^{®} formulation used in this study at these higher shear rates, indicated by the development of a slight overshoot during the start-up flow. In other words, this would mean that this Carbopol^{®} formulation only acts as a true simple yield-stress fluid up to shear rates of $ \gamma \u02d9=0.1 s \u2212 1$.

### E. Three-dimensional constitutive equation

Full three-dimensional elastoviscoplastic constitutive modeling of the deformation behavior of yield-stress and thixotropic materials has been presented before [20,69]. Dimitriou and McKinley [20] formulated a standard elastoviscoplastic equation for Carbopol, including kinematic hardening and de Souza Mendes and Thompson [42] used the SMD-viscosity model to describe a strain-dependent elastoviscoplastic constitutive equation.

*d*” and “

*T*” denote the deviatoric part and the transpose, respectively. Finally, the equivalent stress, $ \sigma eq$ is given by Eq. (4) and $ D t= \u2202 t+\u2207\u22c5 v$ is the material derivative. According to Eq. (20), the hydrostatic deformation is assumed to be elastic, whereas the accumulation of deviatoric elastic strain $ B ~ e$ is retarded by a stress-dependent plastic-strain rate $ L p$. It should be noted that Eq. (20) is thermodynamically admissible, as was shown by Hütter and Tervoort [71], who derived this class of elastoviscoplastic equations using the GENERIC formalism [73]. It should also be noted that for plastic deformation of isotropic materials, the plastic spin is zero, leaving the plastic velocity gradient symmetric and equal to the plastic strain-rate tensor $ L p= D p$ [27,70,74].

^{®}dispersions, not only the plastic but also the elastic distortions proceed to a good approximation at constant volume, so that the “incompressiblity” constraint can be applied. According to this constraint, the hydrostatic stress is solely determined by the boundary conditions and can be excluded from the constitutive relation. In this case, $ tr( D)=0$, $J=1$, $ B ~ e= B e$ with $det( B e)=1$ and $ L= L d$ so that the constitutive relation simplifies to

*elastic*strains, $ \gamma e$, and small applied shear rate, $ \gamma \u02d9$, the evolution equations, Eq. (21), reduce to the one-dimensional constitutive model, Eq. (12) (see Appendix A). In addition, it can be shown that in this case the first normal stress coefficient, $ N 1$, is predicted to be proportional to $ \gamma e 2$: $ N 1= \sigma 11\u2212 \sigma 22=G \gamma e 2$, that the second normal stress coefficient, $ N 2$ equals zero: $ N 2= \sigma 22\u2212 \sigma 33=0$ and that both normal stress coefficients are expected to be small at small elastic strains. Finally, Eq. (21) predicts that, again at small elastic strains in simple shear, the normal stress coefficients have little effect on the equivalent stress $ \sigma eq$ (which is proportional to the von Mises stress) in pure shear. Except for the second normal stress coefficient, all this agrees with some recent experimental observations [75], although other results suggest more complex chaotic behavior [76]. Experimental verification of normal stresses and the validity of the von Mises criterion for the particular grade of Carbopol

^{®}used in this work will be presented in a future publication.

## V. CONCLUSIONS

An elastoviscoplastic constitutive equation is proposed to describe the rate-dependent plastic deformation behavior of Carbopol^{®} dispersions, formulated to be an example of an ideal simple yield-stress fluid. Notably, this model sidesteps the traditional use of a discrete Herschel–Bulkley yield stress. Instead, it adopts a nonlinear Maxwell framework, combining a linear spring with a strongly nonlinear stress-dependent viscosity, implying that the relaxation behavior is governed by a stress-dependent relaxation time. Two key experiments were conducted: measuring the steady-state response following a start-up flow with a constant applied shear rate and a creep test with a constant applied shear stress. It was found that the stress-dependent flow behavior was accurately represented by the well-known Ree–Eyring viscosity model. The constitutive equation was found to describe a continuous transition from elastic to viscous flow during start-up flow experiments, and accurately reproduces the experimentally determined steady-state flow stresses. Moreover, the stress-activated plastic flow is consistent with thermodynamic arguments for the excess entropy [40].

In addition, it was shown for Carbopol^{®} dispersions that it is possible to construct a master creep compliance curve at a given reference shear stress from creep compliance curves measured at various shear stress levels, by mere horizontal shifting using a logarithmic time axis. This implies that to a good approximation, all relaxation times in Carbopol^{®} are the same function of the applied shear stress, leading to *time-stress superposition* analogous to the well-known time-temperature superposition, where it is assumed that all relaxation times are the same function of temperature. The experimental shift factors obtained from building the master creep compliance curve were shown to follow the Ree–Eyring function and agreed well with shift factors obtained from the steady-state stresses, measured during start-up and creep flow experiments.

To critically evaluate the applicability of time-stress superposition, orthogonal-superposition experiments were performed measuring the linear viscoelastic response of Carbopol^{®} dispersions during steady-state shear flow at a given shear stress. Also, for these experiments, it was shown to be possible to build master curves for the orthogonal linear viscoelastic storage and loss modulus ( $ G \u2032$ and $ G \u2033$) measured at different applied shear stresses, using a logarithmic frequency axis, applying the shift factors obtained from the fitting flow stresses in Fig. 4 to the Ree–Eyring model. The applicability of stress-time superposition was then successfully used to perform full nonlinear viscoelastic modeling of the start-up flows up to shear rates in the order of $ \gamma \u02d9=0.1 s \u2212 1$, using the complete relaxation-time spectrum, where all relaxation times shift in the same manner with the total stress. At higher shear rates, the predicted transition to fully developed plastic flow is slower than what is observed experimentally. It is postulated that at these strain rates, the Carbopol^{®} dispersion used in this study is not a simple yield-stress fluid anymore but starts to behave like a thixotropic material.

Finally, a full three-dimensional elastoviscoplastic for finite elastic and total deformations is proposed, which reduces to the one-dimensional Maxwell model used in the first part of this research and which assumes von-Mises-like flow behavior.

## ACKNOWLEDGMENTS

The authors thank Dr. Kirill SOJ Feldman (ETH Zurich) for experimental support, and Professor R. Bonnecaze (UT Austin) and Professor G. McKinley (MIT) for stimulating discussions. Anton Paar (Dr. J. Lauger) is gratefully acknowledged for the use of the OSP prototype. This work was supported by ETH zurich and the Swiss National Science Foundation, SNSF Grant No. 200020-192336 and the International Fine Particle Research Institute (IFPRI).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: RELATION BETWEEN THE 3D AND 1D CONSTITUTIVE MODEL

In the approximation of small elastic strain, $ \gamma e$, and small applied shear rate, $ \gamma \u02d9$, the (1,2) components of Eqs. (A2) and (A3) reduce to the one-dimensional constitutive equation, Eq. (12).

From Eq. (A3a) it follows that the first normal stress coefficient, $ N 1$, is expected to be proportional to $ \gamma e 2$: $ N 1= \sigma 11\u2212 \sigma 22=G \gamma e 2$ and that the second normal stress coefficient $ N 2$ equals zero and $ N 2= \sigma 22\u2212 \sigma 33=0$ and that both normal stress coefficients are expected to be small at small *elastic* strains. From Eq. (A3b) it follows that, again at small elastic strains, the normal stress coefficients have little effect on the equivalent stress $ \sigma eq$ (which is proportional to the von Mises stress) in pure shear. Except for the second normal stress coefficient, all this agrees with experimental observations [75].

### APPENDIX B: OSP LIMITS

*et al.*[49] were used to perform the transformation and from those the corrections that affect the final moduli are defined, to be precise the inertia correction (In: $ ( m + \beta A \rho ) \omega 2$), motor friction correction (Fr: $\zeta \omega $), and motor elasticity correction (El: $K$). In Eq. (B1), s indicate the displacement and F the raw force, while A and $\beta $ are geometrical parameter for which the exact formulas can be found in the cited paper,

The instrument limits the acquisition of the data between 0.94 and 126 rad/s. However, the most impactful limitations are given by the inertia correction and the minimum force detectable (that corresponds $ G min= 22.5 Pa$).

It is clear from the previous plot how the inertia correction affects the data even before the nominal inertia line is met. To avoid this deviation, the data are cut when 10% of the material response is given by the inertial response, Fig. 17.

The curves at 141 and at $164 Pa$ show a clear crossover frequency and its shift at higher frequency as the applied stress increase. However, these stresses correspond to a shear rate of $4.84$ and $12.1 s \u2212 1$, respectively, that are found to be out of the focus range of our model. Therefore, these are excluded from the body of the paper.

## REFERENCES

*Curr. Opin. Colloid Interface Sci.*

*Handbook of Numerical Analysis*

*Rheology of Disperse Systems*, edited by C. Mill (Pergamon, London, 1959), pp. 84–104.

*Dynamics of Polymeric Liquids*

*Polymers: from Structure to Properties*

*Mechanical Properties of Solid Polymers*

^{®}dispersions in water and in water/glycerol solutions

*Elastic and Creep Properties of Filamentous Materials and Other High Polymers*

*Physical Aging of Amorphous Polymers and Other Materials*

*The Mechanics and Thermodynamics of Continua*