No yield stress required: Stress-activated flow in simple yield-stress fluids

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.

Here, $σ d$ is the deviatoric part of the Cauchy stress tensor, $D p$ is the plastic strain-rate tensor, and $γ ˙ eq$ is the equivalent shear rate, which depends on the second invariant of the plastic strain-rate tensor $D p$ as $γ ˙ eq= 2 D p : D p$ and which reduces to the norm of the shear rate $γ ˙$ 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” $η 0$ when the equivalent shear rate is smaller than a characteristic shear rate $γ ˙ 0$⁠.

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 $σ 0$⁠, depending on the power-law index, the power-law fluid departs from linear behavior $( a σ=1)$ at a steeper slope compared to the Eyring equation. However, especially at $n≈0.1$⁠, the shape of the two curves is quite similar, which means they can be made to approximately overlap by selecting the appropriate $σ 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, $σ eq≪ σ 0$⁠, as the stress-shift factor $a σ=1$⁠, the viscosity equals $η 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 $γ ˙ 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.

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.

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², 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: $∼0.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 °$ pitched blade turbine (diameter: 40 mm, power number: $∼0.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 μ 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 (⁠ $σ=$ 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×250×250 μ 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 (⁠ $γ ˙=1 s − 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.

Transient measurements performed at shear rates of $γ ˙= 5 × 10 − 4 s$ to 1 s⁻¹ 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 $σ$ during a start-up experiment at constant applied shear rate $γ ˙$⁠, a steady-state viscosity $η$ can be calculated as $η=σ/ γ ˙$ [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 $σ( γ ˙ e q)=η ( γ ˙ e q ) γ ˙ e q$ and $σ 0= η 0 γ ˙ 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 $σ 01=( 2.7 ± 0.2) Pa$⁠, $σ 02=10.6±0.7 Pa$⁠, $log⁡( γ ˙ 01/ s − 1)=−14.7±0.7$⁠, and $log⁡( γ ˙ 02/ s − 1)=−0.34±0.08$⁠, from which it can be calculated, using Eq. (11), that $η 01=1.59× 10 15$ and $η 02=3.0× 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 $®$⁠. This has been done for various applied shear rates, ranging from $γ ˙=5× 10 − 4$ to 1 s $− 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].

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 σ 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).

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 ′⊥$⁠) and the orthogonal loss moduli (G” $⊥$⁠). It should be noted that in the absence of stress, the Carbopol^®940 sample behaves elastically, meaning that $G ′=G≈300 Pa$ and $G ″≈0 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.

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 $γ ˙ 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 $γ ˙=0.1 s − 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 $γ ˙=0.1 s − 1$⁠.

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.

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 ′$ and $G ″$⁠) 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 $γ ˙=0.1 s − 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.

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).

The authors have no conflicts to disclose.

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

In the approximation of small elastic strain, $γ e$⁠, and small applied shear rate, $γ ˙$⁠, 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 $γ e 2$⁠: $N 1= σ 11− σ 22=G γ e 2$ and that the second normal stress coefficient $N 2$ equals zero and $N 2= σ 22− σ 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 $σ 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].

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 − 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.